# A Hybrid Algorithm Based on Optimal Quadratic Spline Collocation and Parareal Deferred Correction for Parabolic PDEs

Authors: Jun Liu; Yan Wang; Rongjian Li
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Parareal is a kind of time parallel numerical methods for time-dependent systems. In this paper, we consider a general linear parabolic PDE, use optimal quadratic spline collocation (QSC) method for the space discretization, and proceed with the parareal technique on the time domain. Meanwhile, deferred correction technique is also used to improve the accuracy during the iterations. In fact, the optimal QSC method is a correction of general QSC method. Along the temporal direction we embed the iterations of deferred correction into parareal to construct a hybrid method, parareal deferred correction (PDC) method. The error estimation is presented and the stability is analyzed. To save computational cost, we find out a simple way to balance the two kinds of iterations as much as possible. We also argue that the hybrid algorithm has better system efficiency and costs less running time. Numerical experiments by multicore computers are attached to exhibit the effectiveness of the hybrid algorithm.


---

## Body

## 1. Introduction

The parareal algorithm was firstly proposed by Lions et al. [1] in 2001 as a parallel-in-time approach for solving time-dependent differential equations, mainly for ordinary differential equations (ODEs) and parabolic equations [2–4]. The significant advantage of the parareal algorithm is that it does not require the consecutive solving of subproblems along the temporal direction. As a purely parallel algorithm, the parareal algorithm has been applied to many practical problems, such as hyperbolic problems [5], fluid mechanics [6], quantum control [7, 8], and optimized control problem [9, 10]. Recently, some innovative variants of the parareal were proposed by integrating with other efficient methods, such as spectral deferred corrections [11], domain decomposition technique [12], the Krylov subspace method [13], and adaptive technique [14]. We have also proposed a kind of parareal waveform relaxation methods for ODEs [15] and parabolic PDEs [16], as well as parareal deferred correction methods for PDEs [17].The parareal algorithm was initially derived from multiple shooting, time multigrid, and algebraic deferred correction [18, 19] through borrowing the idea of domain decomposition on the time domain. In particular, it employs two operators on two time decompositions of different step sizes with the name of coarse and fine grids to propagate the values along the time direction. The initial guesses of subproblems defined on the coarse grid are given by the coarse propagator for the first iteration. Then, the values at the interface of two adjacent coarse-grid subproblems are updated using fine propagators in iterations, which is quite similar to the Schwarz methods for elliptic systems. Spectral deferred correction method, a kind of iterative methods for solving ODEs, is originally proposed in [20] and has been combined with parareal in [11]. It makes use of the previously computed solution on the same time subdomain. High order approximations can be achieved while the computation cost is reduced sharply.In this paper, we apply the idea of correction to the problem of parabolic PDEs. In detail, we employ the optimal QSC method for the spatial variables and PDC method for the temporal variables, which will result in high accuracy approximations with low computational cost. The new method can also be implemented in parallel, just like the same parallel mechanism as the classical parareal method.The outline of our paper is as follows. In Section2, we give a brief review about parareal algorithm and deferred correction. In Section 3, we employ the optimal QSC method for the space discretization of parabolic PDE, which results in an ODEs system. Then we present the PDC method and apply the method to the ODEs system. In Section 4, we give a proper error estimation for the PDC method. In Section 5, we analyze the stability of the algorithm and show the relationship between the stability and the propagators employed. In Section 6, we present a simple way to choose the time step of the fine propagator. The speedup and system efficiency of the new algorithm are analyzed in Section 7. Numerical experiments are attached in Section 8.


## 2. Materials and Methods

### 2.1. Parareal Algorithm

We briefly review the parareal algorithm on the following nonlinear ODEs:(1)
d
u
t
d
t
=
f
u
t
,
t
,
t
∈
0
,
T
,
u
0
=
u
0
,where the nonlinear function f
:
R
n
×
R
→
R
n is Lipschitz on R
n and the function u
:
R
→
R
n is to be computed.The parareal algorithm requires a decomposition of the time domain[
0
,
T
] into N time subdomains [
T
n
,
T
n
+
1
], n
=
0,1
,
…
,
N
-
1, with 0
=
T
0
<
T
1
<
⋯
<
T
N
-
1
<
T
N
=
T. For simplicity, we suppose that all time subdomains have the same length Δ
T
≔
T
n
+
1
-
T
n. The approximations at these points {
T
n
∣
n
=
1,2
,
…
,
N
} can be obtained using the following scheme:(2)
U
0
k
+
1
=
u
0
,
U
n
+
1
k
+
1
=
G
T
n
+
1
,
T
n
,
U
n
k
+
1
+
F
T
n
+
1
,
T
n
,
U
n
k
-
G
T
n
+
1
,
T
n
,
U
n
k
,where G
(
T
n
+
1
,
T
n
,
U
n
) and F
(
T
n
+
1
,
T
n
,
U
n
) denote, respectively, the coarse and fine approximation at the point T
n
+
1 with initial value U
n at the point T
n. The coarse propagator G is cheap and lowly accurate, and the fine propagator F is of high precision but also more expensive. The initialization of the iterative algorithm (2) is usually realized by the propagation of the coarse propagator along {
T
n
∣
n
=
1,2
,
…
,
N
} consecutively; that is, U
n
+
1
0
=
G
(
T
n
+
1
,
T
n
,
U
n
0
).The parareal algorithm has several well-known properties in the following: (a) the initial values for the fine propagator on each time subdomain are known in every iteration, so the fine propagators on different time subdomains can be implemented by different processors. (b) The final result of parareal algorithm will achieve the accuracy of serial implement of the fine propagator on the whole time domain afterN iterations, where N is the number of time subdomains. Note that the number of iterations is usually set to be much less than the number of time subdomains in practice to avoid necessary computations, while achieving a competitive accuracy [13].


### 2.2. Deferred Correction Method

Deferred correction method is an iterative method for ODEs. We employ the scheme designed in [21] as one of the basic elements for our proposed method, in spite of many other variants of the original deferred corrections. We adopt the scheme in [21] because of its simplicity for implementation with high accuracy and also robustness for different problems.We also take system (1) as our example to describe the execution details of deferred correction method. We firstly define a grid on [
0
,
T
] with (
M
+
1
) equispaced nodes, where 0
=
t
0
<
t
1
<
⋯
<
t
M
=
T, with time step h
=
T
/
M. An rth order method on this grid could bring a numerical solution {
U
0
,
U
1
,
…
,
U
M
} with error U
i
=
u
(
t
i
)
+
O
(
h
r
). In addition, we can also generate a unique Mth-degree polynomial u
~
(
t
) on [
0
,
T
] with this grid using Lagrange interpolation. Then we can define an error function (3)
δ
t
=
u
t
-
u
~
t
,
t
∈
0
,
T
,which satisfies the differential equation(4)
d
δ
t
d
t
=
f
δ
t
+
u
~
t
,
t
-
d
u
~
t
d
t
,
t
∈
0
,
T
,
δ
0
=
0
.This equation is usually regarded as the neighbor problem of original system (1) in the regime of deferred correction methods. Meanwhile, if we use the same rth numerical method to solve this neighbor problem, we could get another group of numerical values {
Δ
0
,
Δ
1
,
…
,
Δ
M
} on the same grid with error Δ
i
=
δ
(
t
i
)
+
O
(
h
r
). The numerical solution {
U
i
+
Δ
i
,
i
=
0,1
,
…
,
M
} is a better approximation of the original system after one correction since we have U
i
+
Δ
i
=
u
(
t
i
)
+
O
(
h
2
r
) [22–24]. If we iterate such correction J times, the order of accuracy of the iterated deferred correction is O
(
h
m
i
n
[
(
J
+
1
)
r
,
M
]
). In practice, the degree M can not be very large.


## 2.1. Parareal Algorithm

We briefly review the parareal algorithm on the following nonlinear ODEs:(1)
d
u
t
d
t
=
f
u
t
,
t
,
t
∈
0
,
T
,
u
0
=
u
0
,where the nonlinear function f
:
R
n
×
R
→
R
n is Lipschitz on R
n and the function u
:
R
→
R
n is to be computed.The parareal algorithm requires a decomposition of the time domain[
0
,
T
] into N time subdomains [
T
n
,
T
n
+
1
], n
=
0,1
,
…
,
N
-
1, with 0
=
T
0
<
T
1
<
⋯
<
T
N
-
1
<
T
N
=
T. For simplicity, we suppose that all time subdomains have the same length Δ
T
≔
T
n
+
1
-
T
n. The approximations at these points {
T
n
∣
n
=
1,2
,
…
,
N
} can be obtained using the following scheme:(2)
U
0
k
+
1
=
u
0
,
U
n
+
1
k
+
1
=
G
T
n
+
1
,
T
n
,
U
n
k
+
1
+
F
T
n
+
1
,
T
n
,
U
n
k
-
G
T
n
+
1
,
T
n
,
U
n
k
,where G
(
T
n
+
1
,
T
n
,
U
n
) and F
(
T
n
+
1
,
T
n
,
U
n
) denote, respectively, the coarse and fine approximation at the point T
n
+
1 with initial value U
n at the point T
n. The coarse propagator G is cheap and lowly accurate, and the fine propagator F is of high precision but also more expensive. The initialization of the iterative algorithm (2) is usually realized by the propagation of the coarse propagator along {
T
n
∣
n
=
1,2
,
…
,
N
} consecutively; that is, U
n
+
1
0
=
G
(
T
n
+
1
,
T
n
,
U
n
0
).The parareal algorithm has several well-known properties in the following: (a) the initial values for the fine propagator on each time subdomain are known in every iteration, so the fine propagators on different time subdomains can be implemented by different processors. (b) The final result of parareal algorithm will achieve the accuracy of serial implement of the fine propagator on the whole time domain afterN iterations, where N is the number of time subdomains. Note that the number of iterations is usually set to be much less than the number of time subdomains in practice to avoid necessary computations, while achieving a competitive accuracy [13].


## 2.2. Deferred Correction Method

Deferred correction method is an iterative method for ODEs. We employ the scheme designed in [21] as one of the basic elements for our proposed method, in spite of many other variants of the original deferred corrections. We adopt the scheme in [21] because of its simplicity for implementation with high accuracy and also robustness for different problems.We also take system (1) as our example to describe the execution details of deferred correction method. We firstly define a grid on [
0
,
T
] with (
M
+
1
) equispaced nodes, where 0
=
t
0
<
t
1
<
⋯
<
t
M
=
T, with time step h
=
T
/
M. An rth order method on this grid could bring a numerical solution {
U
0
,
U
1
,
…
,
U
M
} with error U
i
=
u
(
t
i
)
+
O
(
h
r
). In addition, we can also generate a unique Mth-degree polynomial u
~
(
t
) on [
0
,
T
] with this grid using Lagrange interpolation. Then we can define an error function (3)
δ
t
=
u
t
-
u
~
t
,
t
∈
0
,
T
,which satisfies the differential equation(4)
d
δ
t
d
t
=
f
δ
t
+
u
~
t
,
t
-
d
u
~
t
d
t
,
t
∈
0
,
T
,
δ
0
=
0
.This equation is usually regarded as the neighbor problem of original system (1) in the regime of deferred correction methods. Meanwhile, if we use the same rth numerical method to solve this neighbor problem, we could get another group of numerical values {
Δ
0
,
Δ
1
,
…
,
Δ
M
} on the same grid with error Δ
i
=
δ
(
t
i
)
+
O
(
h
r
). The numerical solution {
U
i
+
Δ
i
,
i
=
0,1
,
…
,
M
} is a better approximation of the original system after one correction since we have U
i
+
Δ
i
=
u
(
t
i
)
+
O
(
h
2
r
) [22–24]. If we iterate such correction J times, the order of accuracy of the iterated deferred correction is O
(
h
m
i
n
[
(
J
+
1
)
r
,
M
]
). In practice, the degree M can not be very large.


## 3. QSC-PDC Methods for Parabolic PDEs

In this study, we consider the following linear PDE as our basic model to demonstrate our proposed method:(5)
∂
∂
t
u
x
,
t
=
∂
2
∂
x
2
u
x
,
t
+
f
x
,
t
,
x
∈
Ω
,
t
∈
0
,
T
,
u
x
,
t
=
0
,
x
∈
∂
Ω
,
t
∈
0
,
T
,
u
x
,
0
=
g
x
,
x
∈
Ω
,where Ω
≡
[
a
,
b
] is the spatial domain, [
0
,
T
] is the time domain with T
>
0, and u
(
x
,
t
) is the unknown function. We take the time domain decomposition for system (5) and design the generic iteration scheme:(6)
U
0
k
+
1
x
=
g
x
,
U
n
+
1
k
+
1
x
=
G
T
n
+
1
,
T
n
,
U
n
k
+
1
x
+
F
T
n
+
1
,
T
n
,
U
n
k
x
-
G
T
n
+
1
,
T
n
,
U
n
k
x
.The essential efforts of this paper are concentrated on the design of techniques, which could bring execution cost and high accuracy for scheme (6).


### 3.1. QSC Methods for Spatial Variables

We firstly consider reducing the computational cost in spatial domain. Note that, for every iterationk, we have to use the propagators G and F to solve a linear homogeneous PDE defined on a fixed spatial domain the same as that of the original system for every time subdomain. The discretization of the spatial and temporal variables will result in algebraic equations for subproblems. The approximate order of the discretization procedure is crucial to the scale of the algebraic equation. In this paper, we use the optimal QSC method for the spatial variable, which was originally used for two-point boundary value problem [25], then elliptic PDEs [26], and parabolic PDEs [27, 28].Similarly, we consider a uniform partition of the spatial domain,{
a
=
x
0
<
x
1
<
⋯
<
x
I
=
b
} with mesh size Δ
x
=
1
/
I. We denote P
2
,
Δ as the space of piecewise quadratic polynomials and S
2
,
Δ
=
P
2
,
Δ
∩
C
1
(
[
a
,
b
]
) as the space of quadratic splines. Remember that the dimension of the space S
2
,
Δ is I
+
2. We let {
ϕ
0
,
ϕ
1
,
…
,
ϕ
I
+
1
} be a set of piecewise polynomial basis functions of the space S
2
,
Δ.For a fixedt, system (5) is a two-point boundary value problem, whose quadratic spline approximate solution can be written as (7)
p
x
=
∑
i
=
0
I
+
1
c
i
ϕ
i
x
.Therefore, we can write the approximation of system (5) as (8)
u
Δ
x
,
t
=
∑
i
=
0
I
+
1
c
i
t
ϕ
i
x
,where c
i
(
t
), i
=
0,1
,
…
,
I, are degrees of freedom (DOFs). Taking the set {τ
0
=
a, τ
i
=
(
x
i
-
1
+
x
i
)
/
2; 1
≤
i
≤
I, τ
I
+
1
=
b} as the collocation points, we can obtain ODEs systems with respect to the DOFs.Starting from the quadratic polynomial(9)
ϕ
x
=
1
2
x
2
,
0
≤
x
≤
1
,
-
2
x
-
1
2
+
2
x
-
1
+
1
,
1
≤
x
≤
2
,
3
-
x
2
,
2
≤
x
≤
3
,
0
,
e
l
s
e
w
h
e
r
e
,we choose a set of quadratic B-splines (10)
ϕ
i
x
=
ϕ
x
-
a
Δ
x
-
i
+
2
,
i
=
0,1
,
…
,
I
+
1as the basis functions of the space S
2
,
Δ. Following [27], we denote S
~
2
,
Δ as the space of quadratic splines satisfying homogeneous Dirichlet boundary conditions. The dimension of the space S
~
2
,
Δ is I, and a set of basis functions of such space is {ϕ
~
1
=
ϕ
1
-
ϕ
0, ϕ
~
i
=
ϕ
i, i
=
2
,
…
,
I
-
1, ϕ
~
I
=
ϕ
I
-
ϕ
I
+
1}. Then the quadratic spline approximate solution of system (5) can be written as (11)
u
Δ
x
,
t
=
∑
i
=
1
I
c
i
t
ϕ
~
i
x
.With the values of the basis functions on collocation points {
τ
1
,
…
,
τ
I
}, we can obtain the following ODEs systems on the DOFs c
i
(
t
) of the optimal one-step QSC methods:(12)
Q
0
d
c
t
d
t
=
1
Δ
x
2
Q
2
+
1
24
Q
x
x
c
t
+
f
→
t
,
t
∈
0
,
T
,
c
0
=
c
0
,where c
(
t
)
=
[
c
1
(
t
)
,
…
,
c
I
(
t
)
]
T, the coefficient matrices are(13)
Q
0
=
1
8
5
1
0
⋯
0
1
6
1
⋯
0
⋱
⋱
⋱
0
⋯
1
6
1
0
⋯
0
1
5
,
Q
2
=
-
3
1
0
⋯
0
1
-
2
1
⋯
0
⋱
⋱
⋱
0
⋯
1
-
2
1
0
⋯
0
1
-
3
,
Q
x
x
=
-
11
16
-
14
6
-
1
0
⋯
0
-
5
6
-
4
1
0
0
⋯
0
1
-
4
6
-
4
1
0
⋯
0
⋱
⋱
⋱
⋱
⋱
0
⋯
0
1
-
4
6
-
4
1
0
⋯
0
0
1
-
4
6
-
5
0
⋯
0
-
1
6
-
14
16
-
11
,and f
→
(
t
) is the interpolating vector of f
(
x
,
t
) at all the collocation points. The vector c
0
∈
R
I satisfies Q
0
c
0
=
g
0, with g
0 the interpolation of g
(
x
) at collocation points.The discretization error of the spatial variable in the uniform norm isO
(
Δ
x
3
) globally and O
(
Δ
x
4
) at the nodes of the uniform partition and the collocation points [25]. In contrast, the classical finite difference discretization could only generate approximations with error O
(
Δ
x
2
) locally, which demonstrates the great benefits of the optimal QSC.


### 3.2. PDC Method for the Temporal Variables

The parareal algorithm has a very important property. The approximation at all the endpoints of time subdomains will achieve the accuracy of theF-propagator after N iterations [8, 18], where N is the number of time subdomains. In fact, for any n
<
N, the approximation at the time point T
n will be updated n times and reach the final approximation which is exactly the result by the F-propagator at this time point.In order to carry out the PDC method on ODEs system (12) with respect to DOFs, we discretize each time interval I
n
=
[
T
n
,
T
n
+
1
] into smaller equispaced grids T
n
=
t
0
<
t
1
<
⋯
<
t
m
<
⋯
<
t
M
=
T
n
+
1 and denote C
n
k the PDC approximation at the point T
n in kth parareal iteration. For simplicity, we only arrange one correction in each parareal iteration. The details of the hybrid algorithm are described as follows step by step:(i)
Use the coarse propagatorG to obtain rough numerical approximations at all the time interval points T
n, that is, C
n
0, n
=
0
,
…
,
N.
(ii)
Use the prediction step of the new propagatorF
D
C on each time interval I
n, usually an (
r
0
)th-order numerical method, to obtain fine approximations at the grid points on I
n with the initial value C
n
0 at T
n. We suppose that the numerical solution is η
→
n
[
1
]
=
(
η
n
,
0
1
,
η
n
,
1
1
,
…
,
η
n
,
m
[
1
]
,
…
,
η
n
,
M
[
1
]
), and F
D
C
(
T
n
+
1
,
T
n
,
C
n
0
)
=
η
n
,
M
[
1
], where η
n
,
0
[
1
]
=
C
n
0.
(iii)
Use the inductive formulae of parareal algorithm to obtain improved approximationsC
n
1 at T
n, n
=
0
,
…
,
N; that is,(14)
C
0
1
=
c
0
,
C
n
+
1
1
=
G
T
n
+
1
,
T
n
,
C
n
1
+
F
D
C
T
n
+
1
,
T
n
,
C
n
0
-
G
T
n
+
1
,
T
n
,
C
n
0
.
(iv)
Fork
=
2
,
…
,
k
0, (k
0 is the number of parareal iterations)
(1)
Solve the following neighbor problem of the original system:(15)
Q
0
d
e
n
k
-
1
t
d
t
=
Q
0
d
c
n
t
d
t
-
Q
0
d
η
n
k
-
1
t
d
t
=
1
Δ
x
2
Q
2
+
1
24
Q
x
x
e
n
k
-
1
t
+
η
n
k
-
1
t
+
f
→
t
-
Q
0
d
η
n
k
-
1
t
d
t
,
t
∈
I
n
,
e
n
k
-
1
T
n
=
C
0
k
-
C
0
k
-
1
,where e
n
(
k
-
1
)
(
t
)
=
c
n
(
t
)
-
η
n
(
k
-
1
)
(
t
) is the error function on each time interval I
n from the previous step, c
n
(
t
) is the true solution of (12) restricted on I
n, and η
n
(
k
-
1
)
(
t
) is an Mth-order polynomial interpolating η
→
n
[
k
-
1
]. We denote the numerical solution of (15) as δ
→
n
[
k
-
1
]
=
(
δ
n
,
0
k
-
1
,
δ
n
,
1
k
-
1
,
…
,
δ
n
,
m
[
k
-
1
]
,
…
,
δ
n
,
M
[
k
-
1
]
).
(2)
Update the numerical solutionη
→
n
[
k
]
=
η
→
n
[
k
-
1
]
+
δ
→
n
[
k
-
1
], and let F
D
C
(
T
n
+
1
,
T
n
,
C
n
k
-
1
)
=
η
n
,
M
[
k
].
(3)
Perform the following iterative formulae to get improved approximations:(16)
C
n
+
1
k
=
G
T
n
+
1
,
T
n
,
C
n
k
+
F
D
C
T
n
+
1
,
T
n
,
C
n
k
-
1
-
G
T
n
+
1
,
T
n
,
C
n
k
-
1
,where n
=
0,1
,
…
,
N
-
1.With these final DOFsC
n
k
0, we can obtain the approximations (17)
U
0
k
0
,
U
1
k
0
,
…
,
U
N
k
0
=
Q
0
C
0
k
0
,
C
1
k
0
,
…
,
C
N
k
0of original system (5) on the grid consisting in spatial collocation points and temporal points.Remark 1.
When solving neighbor problem (15), we can use some simple numerical integrators, such as the backward Euler method. An alternative formulation is to compute the approximate solution in two steps, and the order of accuracy O
(
Δ
x
4
) is also preserved [25, 27, 28]. The realization process of such strategy based on backward Euler method is as follows: (i)
first compute a second order approximationξ
n
,
0
k
-
1
,
…
,
ξ
n
,
m
k
-
1
,
…
,
ξ
n
,
M
+
1
k
-
1 by (18)
Q
0
-
Δ
t
Δ
x
2
Q
2
ξ
n
,
m
+
1
k
-
1
=
Q
0
ξ
n
,
m
k
-
1
+
Δ
t
Δ
x
2
Q
2
η
n
,
m
+
1
k
-
1
+
Δ
t
f
→
t
n
,
m
+
1
-
Δ
t
Q
0
d
d
t
η
n
k
-
1
t
n
,
m
+
1
,
(ii)
then compute an optimal order approximationδ
n
,
0
k
-
1
,
…
,
δ
n
,
m
k
-
1
,
…
,
δ
n
,
M
+
1
k
-
1 by (19)
Q
0
-
Δ
t
Δ
x
2
Q
2
δ
n
,
m
+
1
k
-
1
=
Q
0
δ
n
,
m
k
-
1
+
Δ
t
Δ
x
2
1
24
Q
x
x
ξ
n
,
m
+
1
k
-
1
+
Δ
t
Δ
x
2
Q
2
+
1
24
Q
x
x
η
n
,
m
+
1
k
-
1
+
Δ
t
f
→
t
n
,
m
+
1
-
Δ
t
Q
0
d
d
t
η
n
k
-
1
t
n
,
m
+
1
.The two-step version involves two linear systems over every time step, and both the two coefficient matrices of the above two systems are tridiagonal. So the total computational cost will be smaller than that of one-step integrator (15), in which the bandwidth of the coefficient matrix is 4, because of the matrix Q
x
x.


## 3.1. QSC Methods for Spatial Variables

We firstly consider reducing the computational cost in spatial domain. Note that, for every iterationk, we have to use the propagators G and F to solve a linear homogeneous PDE defined on a fixed spatial domain the same as that of the original system for every time subdomain. The discretization of the spatial and temporal variables will result in algebraic equations for subproblems. The approximate order of the discretization procedure is crucial to the scale of the algebraic equation. In this paper, we use the optimal QSC method for the spatial variable, which was originally used for two-point boundary value problem [25], then elliptic PDEs [26], and parabolic PDEs [27, 28].Similarly, we consider a uniform partition of the spatial domain,{
a
=
x
0
<
x
1
<
⋯
<
x
I
=
b
} with mesh size Δ
x
=
1
/
I. We denote P
2
,
Δ as the space of piecewise quadratic polynomials and S
2
,
Δ
=
P
2
,
Δ
∩
C
1
(
[
a
,
b
]
) as the space of quadratic splines. Remember that the dimension of the space S
2
,
Δ is I
+
2. We let {
ϕ
0
,
ϕ
1
,
…
,
ϕ
I
+
1
} be a set of piecewise polynomial basis functions of the space S
2
,
Δ.For a fixedt, system (5) is a two-point boundary value problem, whose quadratic spline approximate solution can be written as (7)
p
x
=
∑
i
=
0
I
+
1
c
i
ϕ
i
x
.Therefore, we can write the approximation of system (5) as (8)
u
Δ
x
,
t
=
∑
i
=
0
I
+
1
c
i
t
ϕ
i
x
,where c
i
(
t
), i
=
0,1
,
…
,
I, are degrees of freedom (DOFs). Taking the set {τ
0
=
a, τ
i
=
(
x
i
-
1
+
x
i
)
/
2; 1
≤
i
≤
I, τ
I
+
1
=
b} as the collocation points, we can obtain ODEs systems with respect to the DOFs.Starting from the quadratic polynomial(9)
ϕ
x
=
1
2
x
2
,
0
≤
x
≤
1
,
-
2
x
-
1
2
+
2
x
-
1
+
1
,
1
≤
x
≤
2
,
3
-
x
2
,
2
≤
x
≤
3
,
0
,
e
l
s
e
w
h
e
r
e
,we choose a set of quadratic B-splines (10)
ϕ
i
x
=
ϕ
x
-
a
Δ
x
-
i
+
2
,
i
=
0,1
,
…
,
I
+
1as the basis functions of the space S
2
,
Δ. Following [27], we denote S
~
2
,
Δ as the space of quadratic splines satisfying homogeneous Dirichlet boundary conditions. The dimension of the space S
~
2
,
Δ is I, and a set of basis functions of such space is {ϕ
~
1
=
ϕ
1
-
ϕ
0, ϕ
~
i
=
ϕ
i, i
=
2
,
…
,
I
-
1, ϕ
~
I
=
ϕ
I
-
ϕ
I
+
1}. Then the quadratic spline approximate solution of system (5) can be written as (11)
u
Δ
x
,
t
=
∑
i
=
1
I
c
i
t
ϕ
~
i
x
.With the values of the basis functions on collocation points {
τ
1
,
…
,
τ
I
}, we can obtain the following ODEs systems on the DOFs c
i
(
t
) of the optimal one-step QSC methods:(12)
Q
0
d
c
t
d
t
=
1
Δ
x
2
Q
2
+
1
24
Q
x
x
c
t
+
f
→
t
,
t
∈
0
,
T
,
c
0
=
c
0
,where c
(
t
)
=
[
c
1
(
t
)
,
…
,
c
I
(
t
)
]
T, the coefficient matrices are(13)
Q
0
=
1
8
5
1
0
⋯
0
1
6
1
⋯
0
⋱
⋱
⋱
0
⋯
1
6
1
0
⋯
0
1
5
,
Q
2
=
-
3
1
0
⋯
0
1
-
2
1
⋯
0
⋱
⋱
⋱
0
⋯
1
-
2
1
0
⋯
0
1
-
3
,
Q
x
x
=
-
11
16
-
14
6
-
1
0
⋯
0
-
5
6
-
4
1
0
0
⋯
0
1
-
4
6
-
4
1
0
⋯
0
⋱
⋱
⋱
⋱
⋱
0
⋯
0
1
-
4
6
-
4
1
0
⋯
0
0
1
-
4
6
-
5
0
⋯
0
-
1
6
-
14
16
-
11
,and f
→
(
t
) is the interpolating vector of f
(
x
,
t
) at all the collocation points. The vector c
0
∈
R
I satisfies Q
0
c
0
=
g
0, with g
0 the interpolation of g
(
x
) at collocation points.The discretization error of the spatial variable in the uniform norm isO
(
Δ
x
3
) globally and O
(
Δ
x
4
) at the nodes of the uniform partition and the collocation points [25]. In contrast, the classical finite difference discretization could only generate approximations with error O
(
Δ
x
2
) locally, which demonstrates the great benefits of the optimal QSC.


## 3.2. PDC Method for the Temporal Variables

The parareal algorithm has a very important property. The approximation at all the endpoints of time subdomains will achieve the accuracy of theF-propagator after N iterations [8, 18], where N is the number of time subdomains. In fact, for any n
<
N, the approximation at the time point T
n will be updated n times and reach the final approximation which is exactly the result by the F-propagator at this time point.In order to carry out the PDC method on ODEs system (12) with respect to DOFs, we discretize each time interval I
n
=
[
T
n
,
T
n
+
1
] into smaller equispaced grids T
n
=
t
0
<
t
1
<
⋯
<
t
m
<
⋯
<
t
M
=
T
n
+
1 and denote C
n
k the PDC approximation at the point T
n in kth parareal iteration. For simplicity, we only arrange one correction in each parareal iteration. The details of the hybrid algorithm are described as follows step by step:(i)
Use the coarse propagatorG to obtain rough numerical approximations at all the time interval points T
n, that is, C
n
0, n
=
0
,
…
,
N.
(ii)
Use the prediction step of the new propagatorF
D
C on each time interval I
n, usually an (
r
0
)th-order numerical method, to obtain fine approximations at the grid points on I
n with the initial value C
n
0 at T
n. We suppose that the numerical solution is η
→
n
[
1
]
=
(
η
n
,
0
1
,
η
n
,
1
1
,
…
,
η
n
,
m
[
1
]
,
…
,
η
n
,
M
[
1
]
), and F
D
C
(
T
n
+
1
,
T
n
,
C
n
0
)
=
η
n
,
M
[
1
], where η
n
,
0
[
1
]
=
C
n
0.
(iii)
Use the inductive formulae of parareal algorithm to obtain improved approximationsC
n
1 at T
n, n
=
0
,
…
,
N; that is,(14)
C
0
1
=
c
0
,
C
n
+
1
1
=
G
T
n
+
1
,
T
n
,
C
n
1
+
F
D
C
T
n
+
1
,
T
n
,
C
n
0
-
G
T
n
+
1
,
T
n
,
C
n
0
.
(iv)
Fork
=
2
,
…
,
k
0, (k
0 is the number of parareal iterations)
(1)
Solve the following neighbor problem of the original system:(15)
Q
0
d
e
n
k
-
1
t
d
t
=
Q
0
d
c
n
t
d
t
-
Q
0
d
η
n
k
-
1
t
d
t
=
1
Δ
x
2
Q
2
+
1
24
Q
x
x
e
n
k
-
1
t
+
η
n
k
-
1
t
+
f
→
t
-
Q
0
d
η
n
k
-
1
t
d
t
,
t
∈
I
n
,
e
n
k
-
1
T
n
=
C
0
k
-
C
0
k
-
1
,where e
n
(
k
-
1
)
(
t
)
=
c
n
(
t
)
-
η
n
(
k
-
1
)
(
t
) is the error function on each time interval I
n from the previous step, c
n
(
t
) is the true solution of (12) restricted on I
n, and η
n
(
k
-
1
)
(
t
) is an Mth-order polynomial interpolating η
→
n
[
k
-
1
]. We denote the numerical solution of (15) as δ
→
n
[
k
-
1
]
=
(
δ
n
,
0
k
-
1
,
δ
n
,
1
k
-
1
,
…
,
δ
n
,
m
[
k
-
1
]
,
…
,
δ
n
,
M
[
k
-
1
]
).
(2)
Update the numerical solutionη
→
n
[
k
]
=
η
→
n
[
k
-
1
]
+
δ
→
n
[
k
-
1
], and let F
D
C
(
T
n
+
1
,
T
n
,
C
n
k
-
1
)
=
η
n
,
M
[
k
].
(3)
Perform the following iterative formulae to get improved approximations:(16)
C
n
+
1
k
=
G
T
n
+
1
,
T
n
,
C
n
k
+
F
D
C
T
n
+
1
,
T
n
,
C
n
k
-
1
-
G
T
n
+
1
,
T
n
,
C
n
k
-
1
,where n
=
0,1
,
…
,
N
-
1.With these final DOFsC
n
k
0, we can obtain the approximations (17)
U
0
k
0
,
U
1
k
0
,
…
,
U
N
k
0
=
Q
0
C
0
k
0
,
C
1
k
0
,
…
,
C
N
k
0of original system (5) on the grid consisting in spatial collocation points and temporal points.Remark 1.
When solving neighbor problem (15), we can use some simple numerical integrators, such as the backward Euler method. An alternative formulation is to compute the approximate solution in two steps, and the order of accuracy O
(
Δ
x
4
) is also preserved [25, 27, 28]. The realization process of such strategy based on backward Euler method is as follows: (i)
first compute a second order approximationξ
n
,
0
k
-
1
,
…
,
ξ
n
,
m
k
-
1
,
…
,
ξ
n
,
M
+
1
k
-
1 by (18)
Q
0
-
Δ
t
Δ
x
2
Q
2
ξ
n
,
m
+
1
k
-
1
=
Q
0
ξ
n
,
m
k
-
1
+
Δ
t
Δ
x
2
Q
2
η
n
,
m
+
1
k
-
1
+
Δ
t
f
→
t
n
,
m
+
1
-
Δ
t
Q
0
d
d
t
η
n
k
-
1
t
n
,
m
+
1
,
(ii)
then compute an optimal order approximationδ
n
,
0
k
-
1
,
…
,
δ
n
,
m
k
-
1
,
…
,
δ
n
,
M
+
1
k
-
1 by (19)
Q
0
-
Δ
t
Δ
x
2
Q
2
δ
n
,
m
+
1
k
-
1
=
Q
0
δ
n
,
m
k
-
1
+
Δ
t
Δ
x
2
1
24
Q
x
x
ξ
n
,
m
+
1
k
-
1
+
Δ
t
Δ
x
2
Q
2
+
1
24
Q
x
x
η
n
,
m
+
1
k
-
1
+
Δ
t
f
→
t
n
,
m
+
1
-
Δ
t
Q
0
d
d
t
η
n
k
-
1
t
n
,
m
+
1
.The two-step version involves two linear systems over every time step, and both the two coefficient matrices of the above two systems are tridiagonal. So the total computational cost will be smaller than that of one-step integrator (15), in which the bandwidth of the coefficient matrix is 4, because of the matrix Q
x
x.


## 4. The Convergence of the PDC Method

The convergence analysis of the optimal QSC method can be found in the literature [25], and we just need to present the convergence of the PDC algorithm in the temporal domain. Bal [2], Gander, and Hairer [29] gave a convergence analysis on the classical parareal algorithm. Based on their work, we present the convergence analysis on the PDC method.We also assume that the coarse propagatorG satisfies error expansion(20)
S
T
n
+
1
,
T
n
,
U
-
G
T
n
+
1
,
T
n
,
U
=
c
p
+
1
U
Δ
T
p
+
1
+
c
p
+
2
U
Δ
T
p
+
2
+
⋯
,for Δ
T small and the Lipschitz condition(21)
G
t
+
Δ
T
,
t
,
U
-
G
t
+
Δ
T
,
t
,
V
≤
1
+
C
2
Δ
T
U
-
V
,where S
(
T
n
+
1
,
T
n
,
U
) means the true solution at point T
n
+
1 of the original system with the value U at the time point T
n.Within every time subinterval[
T
n
,
T
n
+
1
], the difference between the true solution and the approximation by the PDC method can be written as (22)
u
T
n
+
1
-
U
n
+
1
k
+
1
=
S
T
n
+
1
,
T
n
,
u
T
n
-
G
T
n
+
1
,
T
n
,
U
n
k
+
1
-
F
D
C
T
n
+
1
,
T
n
,
U
n
k
+
G
T
n
+
1
,
T
n
,
U
n
k
=
S
T
n
+
1
,
T
n
,
u
T
n
-
G
T
n
+
1
,
T
n
,
u
T
n
+
G
T
n
+
1
,
T
n
,
U
n
k
-
S
T
n
+
1
,
T
n
,
U
n
k
+
S
T
n
+
1
,
T
n
,
U
n
k
-
F
D
C
T
n
+
1
,
T
n
,
U
n
k
+
G
T
n
+
1
,
T
n
,
u
T
n
-
G
T
n
+
1
,
T
n
,
U
n
k
+
1
.Based on this relationship and also the conditions (20) and (21), we could get(23)
u
T
n
+
1
-
U
n
+
1
k
+
1
≤
C
1
Δ
T
p
+
1
u
T
n
-
U
n
k
+
1
+
C
2
Δ
T
u
T
n
-
U
n
k
+
1
+
S
T
n
+
1
,
T
n
,
U
n
k
-
F
D
C
T
n
+
1
,
T
n
,
U
n
k
.We first look at the last termS
(
T
n
+
1
,
T
n
,
U
n
k
)
-
F
D
C
(
T
n
+
1
,
T
n
,
U
n
k
), which is an upper bound of the error caused by the deferred correction. For simplicity, we assume the solution of the original system to be smooth enough and solve the original equation in the prediction steps and the neighbor equations in the correction steps by the forward Euler integrators.The prediction step in deferred correction includes a numerical integration with initial valueU
n
0 on every time interval I
n. Let η
→
n
[
1
]
=
(
η
n
,
0
[
1
]
,
η
n
,
1
[
1
]
,
…
,
η
n
,
M
[
1
]
) be the numerical solution from the forward Euler integrators on (
M
+
1
) uniformly distributed nodes in I
n with time step h. We know from [30] that the error satisfies u
n
(
M
h
)
-
η
n
,
M
[
1
]
∞
~
O
(
h
), where u
n
(
M
h
) is the true solution at T
n
+
1 with the same initial value U
n
0 at T
n. In fact, the error at T
n
+
1 is related to the length of I
n, so it is reasonable to hypothesize u
n
(
M
h
)
-
η
n
,
M
[
1
]
∞
≤
C
Δ
T
h, where C is a constant.The correction steps in the hybrid algorithm are different from these steps of the pure deferred correction. The initial values of neighbor problems are always0 for pure deferred correction. Meanwhile, for the correction steps which are embedded in parareal process, the initial values of the neighbor problems are kept being updated in parareal iterations. In other words, the initial value of problem (15) on time subinterval I
n is U
0
k
-
U
0
k
-
1. However, this formula of the initial value does not affect the iterative improvement property of the deferred correction. For the kth parareal iteration, we have(24)
u
t
n
,
m
-
η
n
,
m
k
=
u
t
n
,
m
-
η
n
,
m
k
-
1
-
δ
n
,
m
k
-
1
=
e
n
k
-
1
t
n
,
m
-
δ
n
,
m
k
-
1
,where t
n
,
m
=
(
n
Δ
T
+
m
h
). We can see from (24) that the error of deferred correction method equals the error of the corresponding neighbor problem. Similar to the inductive convergence proof of spectral deferred correction in [30], we can obtain that S
(
T
n
+
1
,
T
n
,
U
n
k
)
-
F
D
C
(
T
n
+
1
,
T
n
,
U
n
k
)
~
O
(
h
k
+
1
), considering the solution of the original system to be smooth enough.Based on these results, we have the convergence of the PDC method.Theorem 2.
Let the solution of the original system be smooth enough. If the coarse propagatorG satisfies condition (20) and condition (21), where c
j, j
=
p
+
1
,
p
+
2
,
…, are continuously differentiable, then the error at time point T
n after k PDC iterations is bounded:(25)
u
T
n
-
U
n
k
≤
C
3
C
1
C
1
T
n
Δ
T
p
k
+
1
k
+
1
!
e
C
2
T
n
-
T
k
+
1
+
∑
i
=
0
k
-
1
C
1
i
Δ
T
p
i
T
n
i
i
+
1
!
T
n
e
C
2
T
n
-
T
i
+
1
C
T
h
k
-
i
.Obviously, when the time step h approaches 0, such error bound will go to 0 with enough iterations.Proof.
We suppose that a constantC
T exists such that S
(
T
n
+
1
,
T
n
,
U
n
k
)
-
F
D
C
(
T
n
+
1
,
T
n
,
U
n
k
)
≤
C
T
Δ
T
h
k
+
1 holds for any n
=
0,1
,
…
,
N
-
1. Inequality (23) can be bounded by (26)
u
T
n
+
1
-
U
n
+
1
k
+
1
≤
C
1
Δ
T
p
+
1
u
T
n
-
U
n
k
+
1
+
C
2
Δ
T
u
T
n
-
U
n
k
+
1
+
C
T
Δ
T
h
k
+
1
.
Considering such inequality, we study the recurrence relation(27)
e
n
+
1
k
+
1
=
α
e
n
k
+
β
e
n
k
+
1
+
C
T
Δ
T
h
k
+
1
,
e
n
+
1
0
=
γ
+
β
e
n
0
,where α
=
C
1
Δ
T
p
+
1, β
=
1
+
C
2
Δ
T, and γ
=
C
3
Δ
T
p
+
1. It is easy to see that e
n
+
1
k
+
1 is an upper bound on u
(
T
n
+
1
)
-
U
n
+
1
k
+
1, and we can choose e
0
k
=
0, k
=
0,1
,
…. Multiplying (27) by ζ
n
+
1 and summing over n, we obtain (28)
∑
n
≥
0
e
n
+
1
k
+
1
ζ
n
+
1
=
∑
n
≥
0
α
e
n
k
ζ
n
+
1
+
∑
n
≥
0
β
e
n
k
+
1
ζ
n
+
1
+
∑
n
≥
0
C
T
Δ
T
h
k
+
1
ζ
n
+
1
,
∑
n
≥
0
e
n
+
1
0
ζ
n
+
1
=
∑
n
≥
0
γ
ζ
n
+
1
+
∑
n
≥
0
β
e
n
0
ζ
n
+
1
,where -
1
<
ζ
<
1
/
β
<
1.
We denoteρ
k
(
ζ
)
≔
∑
n
≥
0
e
n
+
1
k
ζ
n
+
1. This function satisfies the following recurrence relation: (29)
ρ
k
+
1
ζ
=
α
ζ
ρ
k
ζ
+
β
ζ
ρ
k
+
1
ζ
+
C
T
Δ
T
h
k
+
1
ζ
1
-
ζ
,
ρ
0
ζ
=
γ
ζ
1
-
ζ
+
β
ζ
ρ
0
ζ
.Then, we have (30)
ρ
k
ζ
=
α
ζ
1
-
β
ζ
ρ
k
-
1
ζ
+
C
T
Δ
T
h
k
ζ
1
-
ζ
1
-
β
ζ
,
ρ
0
ζ
=
γ
ζ
1
-
ζ
1
-
β
ζ
.Solving for ρ
k
(
ζ
), we obtain after induction (31)
ρ
k
ζ
=
α
ζ
1
-
β
ζ
k
ρ
0
ζ
+
C
T
Δ
T
ζ
1
-
ζ
1
-
β
ζ
h
k
+
a
h
k
-
1
+
⋯
+
a
k
-
1
h
=
γ
α
k
ζ
k
+
1
1
-
β
ζ
k
+
1
1
-
ζ
+
C
T
Δ
T
ζ
1
-
ζ
1
-
β
ζ
h
k
+
a
h
k
-
1
+
⋯
+
a
k
-
1
h
,where a
=
α
ζ
/
(
1
-
β
ζ
). We know that replacing the factor 1
-
ζ in the denominator by 1
-
β
ζ only will increase the coefficients in the power series of ρ
k
(
ζ
). It means the coefficients in the power series of ρ
k
(
ζ
) are smaller than the corresponding coefficients in the power series of ρ
~
k
(
ζ
), where(32)
ρ
~
k
ζ
=
γ
α
k
ζ
k
+
1
1
-
β
ζ
k
+
2
+
C
T
Δ
T
ζ
1
-
β
ζ
2
h
k
+
a
h
k
-
1
+
⋯
+
a
k
-
1
h
.
Using the binomial series expansion(33)
1
1
-
β
ζ
k
+
2
=
∑
j
≥
0
k
+
1
+
j
j
β
j
ζ
j
,it is easy to find that the coefficient of the term ζ
n in the expression ρ
~
k
(
ζ
) is (34)
c
n
k
=
γ
α
k
n
k
+
1
β
n
-
k
-
1
+
∑
i
=
0
k
-
1
n
i
+
1
α
i
β
n
-
i
-
1
C
T
Δ
T
h
k
-
i
=
C
3
Δ
T
p
+
1
C
1
k
Δ
T
p
+
1
k
n
!
n
-
k
-
1
!
k
+
1
!
1
+
C
2
Δ
T
n
-
k
-
1
+
∑
i
=
0
k
-
1
n
!
n
-
i
-
1
!
i
+
1
!
C
1
i
Δ
T
p
+
1
i
1
+
C
2
Δ
T
n
-
i
-
1
C
T
Δ
T
h
k
-
i
.For any x
>
0, we have 1
+
x
≤
e
x. So, it follows (35)
c
n
k
≤
C
3
C
1
C
1
Δ
T
p
+
1
k
+
1
k
+
1
!
1
+
C
2
Δ
T
n
-
k
-
1
∏
j
=
0
k
n
-
j
+
∑
i
=
0
k
-
1
C
1
i
Δ
T
p
+
1
i
i
+
1
!
1
+
C
2
Δ
T
n
-
i
-
1
C
T
Δ
T
h
k
-
i
∏
j
=
0
i
n
-
j
≤
C
3
C
1
C
1
T
n
Δ
T
p
k
+
1
k
+
1
!
e
C
2
T
n
-
T
k
+
1
+
∑
i
=
0
k
-
1
C
1
i
Δ
T
p
i
T
n
i
i
+
1
!
T
n
e
C
2
T
n
-
T
i
+
1
C
T
h
k
-
i
.Remark 3.
In practice, we usually takeΔ
T very small and h much smaller. We notice that the first term in (25) decreases to 0 superlinearly as k increases, and such decrease is very quick with the help of small Δ
T. For another hand, the increase of k also takes the summation in (25) much smaller despite more terms. Therefore, a couple of iterations are enough to pull the error down below some tolerance.


## 5. Stability Analysis of the PDC Method

To describe the stability of the PDC method, we apply such numerical method to the equation(36)
d
φ
t
d
t
=
λ
φ
t
,
t
∈
0,1
,
φ
0
=
1
.We follow the definition of stability in [20] and define the amplification factor A
m
(
λ
Δ
T
) for λ
∈
C as (37)
Am
⁡
λ
Δ
T
=
φ
~
1
,where φ
~
(
1
) denotes the numerical approximation at the point t
=
1. The numerical method is stable for some given value of λ, if (38)
Am
⁡
λ
Δ
T
≤
1
.For simplicity, we adopt forward and backward Euler methods in the analysis. The time domain[
0,1
] is divided into N subdomains. The coarse propagator is forced to take large time steps, and implicit Euler is better than explicit Euler for the coarse propagator [31]. Therefore, we choose one-step backward Euler method to be the coarse propagator on every subdomain, while the fine propagator is chosen as s steps of Euler method on every subdomain for the original equations and corresponding neighbor problems in the process of deferred correction.We notice that the approximationφ
~
(
1
) is related to the iteration number k of the PDC method. We denote by Φ
n
k the approximation of system (36) at nth subdomain point and on the kth iteration, where n
=
0,1
,
…
,
N and k
=
0,1
,
…. The approximations satisfy (39)
Φ
n
+
1
k
+
1
=
G
T
n
+
1
,
T
n
,
Φ
n
k
+
1
+
F
T
n
+
1
,
T
n
,
Φ
n
k
-
G
T
n
+
1
,
T
n
,
Φ
n
k
.Because of the one-step backward Euler method for the coarse propagator, we have (40)
G
T
n
+
1
,
T
n
,
Φ
n
k
=
1
1
-
λ
Δ
T
Φ
n
k
.We denote R
=
1
/
(
1
-
λ
Δ
T
) as the amplification factor of the coarse propagator and let r be the amplification factor of the fine propagator, which is from the deferred correction method. In fact, the first iteration of the PDC method is exactly the classical deferred correction, but with smaller numerical steps on each time subdomain. Therefore, r
=
(
1
/
(
1
-
λ
Δ
T
/
s
)
)
s for s steps of backward Euler integrators, and r
=
(
1
+
λ
Δ
T
/
s
)
s for s steps of forward Euler integrators. However, the expressions of r will be much more complicated when k
≥
1.For the initial approximation, we have(41)
Φ
n
0
=
R
n
φ
0
=
R
n
,
n
=
0,1
,
…
,
N
.For the first iteration (k
=
0), we have (42)
Φ
n
+
1
1
=
R
Φ
n
1
+
r
Φ
n
0
-
R
Φ
n
0
=
R
Φ
n
1
+
r
-
R
Φ
n
0
,
n
=
0,1
,
…
,
N
-
1
,and, after inductions, we get (43)
Φ
N
1
=
R
N
-
1
N
r
-
N
-
1
R
,which is also the amplification factor after one iteration of the PDC method. We will continue to analyze the stability in two different cases, respectively.Case 1.
For thes steps of implicit Euler integrators as the fine propagators(44)
A
m
I
m
I
m
λ
Δ
T
,
k
=
0
=
1
1
-
λ
Δ
T
N
-
1
N
1
-
λ
Δ
T
/
s
s
-
N
-
1
1
-
λ
Δ
T
.We denote z
=
λ
Δ
T
∈
C, and the factor A
m
I
m
I
m
(
λ
Δ
T
) is always smaller than 1 except two areas around its two poles z
=
1 and z
=
s. If we choose N
=
20 and s
=
20, the stability regions are shown in Figure 1.
Whenk
≥
1, we can not get the explicit expressions of the amplification factors. The stability regions are estimated numerically in Figure 2. We see that the size of the instable region expands as the iteration number is increased. Nevertheless, these instable regions always lie in the right half plane, so all the PDC methods with implicit fine propagators at these tested iterations are A-stable.Figure 1
Instability region (blank area) using implicit Euler integrators for both the coarse and fine propagators. Thex-axis is R
e
(
λ
Δ
T
) and the y-axis is I
m
(
λ
Δ
T
), where N
=
20 and s
=
20.Figure 2
Instability regions (blank area) of the PDC method using implicit Euler integrators for both the coarse and fine propagators after 2 iterations (a), 3 iterations (b), 4 iterations (c), and 6 iterations (d). Thex-axis is R
e
(
λ
Δ
T
) and the y-axis is I
m
(
λ
Δ
T
), where N
=
20 and s
=
20.
(a)
(b)
(c)
(d)Additionally, we notice that(45)
lim
Re
z
→
∞
⁡
A
m
I
m
I
m
z
,
0
=
0
,and the PDC method with both implicit propagators is L-stable at the first iteration. We take several values for λ along x-axis, and the corresponding amplification factors are in Table 1.Table 1
Amplification factors at different values ofλ.
λ
-
10
-
50
-
100
-
200
-
400
-
800
-
1600
k
=
1
1.68
e
-
04
9.41
e
-
10
4.00
e
-
14
2.52
e
-
19
6.15
e
-
25
9.49
e
-
31
1.16
e
-
36
k
=
2
1.26
e
-
05
3.39
e
-
09
2.09
e
-
13
1.42
e
-
18
3.48
e
-
24
5.38
e
-
30
6.56
e
-
36
k
=
3
4.97
e
-
05
8.46
e
-
09
7.62
e
-
13
5.66
e
-
18
1.39
e
-
23
2.15
e
-
29
2.62
e
-
35
k
=
4
4.55
e
-
05
1.55
e
-
08
2.07
e
-
12
1.69
e
-
17
4.18
e
-
23
6.45
e
-
29
7.87
e
-
35From Table1, we can conjecture that the PDC methods at different iterations could also be L-stable.Case 2.
For thes steps of explicit Euler integrators as the fine propagators, the amplification factor is (46)
A
m
I
m
E
x
λ
Δ
T
,
k
=
0
=
1
1
-
λ
Δ
T
N
-
1
N
1
+
λ
Δ
T
s
s
-
N
-
1
1
-
λ
Δ
T
.The regions of z
=
λ
Δ
T to make A
m
I
m
E
x
(
λ
Δ
T
,
k
=
0
)
≤
1 are related to the relationship between s and N.
Ifs
>
N
-
1, the factor A
m
I
m
E
x
(
z
,
0
) has two poles, that is, z
=
1 and infinity. We choose N
=
20 and s
=
20, and the stability regions of the PDC method at different iterations are shown in Figure 3, in which the stable areas are shrunken as the iteration proceeds.
Ifs
<
N
-
1, the factor A
m
I
m
E
x
(
z
,
0
) has one pole, that is, z
=
1. We choose N
=
20 and s
=
16, the stability regions of the PDC method at different iterations are shown in Figure 4. In fact, there are very small unstable areas located at the left half plane in Figures 4(a) and 4(b), so the PDC methods are not A-stable in such situations either.
Ifs
=
N
-
1, the factor A
m
I
m
E
x
(
z
|
z
|
=
+
∞
,
0
) is always a positive constant much smaller than 1 because the number of time steps s is usually a moderate integer. Therefore, we are reasonable in merging this case into the previous one.
In Figure4, we observe that there is a significant change about the pattern of the size of the stable region at the third iteration. The stability region before the third iteration is quite similar to the stability region of the backward Euler method, while the stability region after the third iteration has a similar shape of the stability region of the forward Euler method. The change means that the dominant control is transferred between the components of the PDC method. Figure 5 illustrates the relationship between this critical iteration and the number of interpolation points s.
We can see a common trend that the stability regions are decreased with more iterations according to Figures2–4 in this section. In fact, there are several factors having impact on the size of the stability region, including the coarse propagators, the fine propagators, the way of interpolation at those equispaced nodes, and the numerical differentiation formula of the original ODEs [32, 33]. We believe that the change of the stable regions is the result of all these factors working together accumulatively.Figure 3
Stability region (doted area) using implicit Euler integrators for the coarse propagators and explicit Euler integrators for the fine propagators after 2 iterations (a), 3 iterations (b), 4 iterations (c), and 6 iterations (d). In fact, there exists a relatively very small unstable area around the pointz
=
1 in (a). The x-axis is R
e
(
λ
Δ
T
) and the y-axis is I
m
(
λ
Δ
T
), where N
=
20 and s
=
20.
(a)
(b)
(c)
(d)Figure 4
Stability region (doted area) using implicit Euler integrators for the coarse propagators and explicit Euler integrators for the fine propagators after 1 iteration (a), 3 iterations (b), 4 iterations (c), and 5 iterations (d). Thex-axis is R
e
(
λ
Δ
T
) and the y-axis is I
m
(
λ
Δ
T
), where N
=
20, and s
=
16.
(a)
(b)
(c)
(d)Figure 5
The relationship between the critical iteration ands
(
<
N
), where N
=
20.


## 6. The Choice of the Fine Time Step

The PDC method includes two iterative processes, that is, parareal and deferred correction. Their individually different iteration speeds will alternatively influence the performance of the proposed hybrid algorithm. To this end, we study the balance strategy to avoid the waste of resources as well as a slowed-down speed. It is intuitively to tune the parameters of deferred correction since they are embedded in the framework parareal for our proposed PDC method. In fact, there are several ways to improve the accuracy of the final approximations for pure deferred corrections, such as smaller time step, higher order numerical integrators, or more iterations. In this study, we choose to tune the time step of the fine propagators of the PDC method.Here we take the following system as an example to show how to find an optimized time step for the PDC method:(47)
d
u
t
d
t
=
A
u
t
,
t
∈
0
,
T
,
u
0
=
U
0
.We assume the G propagator is one-step backward Euler integrator on I
n with G
(
T
n
+
1
,
T
n
,
U
n
k
)
=
(
I
-
A
Δ
T
)
-
1
U
n
k, and the fine propagator is the M-step backward Euler integrator for prediction and each correction loop. Therefore, we have (48)
u
T
n
+
1
-
U
n
+
1
k
+
1
=
e
A
Δ
T
u
T
n
-
G
T
n
+
1
,
T
n
,
U
n
k
+
1
-
F
D
C
T
n
+
1
,
T
n
,
U
n
k
+
G
T
n
+
1
,
T
n
,
U
n
k
=
e
A
Δ
T
u
T
n
-
e
A
Δ
T
U
n
k
+
e
A
Δ
T
U
n
k
-
F
D
C
T
n
+
1
,
T
n
,
U
n
k
+
G
T
n
+
1
,
T
n
,
u
T
n
-
G
T
n
+
1
,
T
n
,
U
n
k
+
1
-
G
T
n
+
1
,
T
n
,
u
T
n
-
G
T
n
+
1
,
T
n
,
U
n
k
=
e
A
Δ
T
-
I
-
A
Δ
T
-
1
u
T
n
-
U
n
k
+
I
-
A
Δ
T
-
1
u
T
n
-
U
n
k
+
1
+
O
h
k
+
1
,where h is the fine time step. Let α
¯
=
e
A
Δ
T
-
(
I
-
A
Δ
T
)
-
1 and β
¯
=
(
I
-
A
Δ
T
)
-
1; then an upper bound of the error can be obtained from the inductive equation(49)
e
¯
n
+
1
k
+
1
=
α
¯
e
¯
n
k
+
β
¯
e
¯
n
k
+
1
+
O
h
k
+
1
.We rewrite the recurrence relation in form of vectors, which is (50)
e
¯
1
k
+
1
e
¯
2
k
+
1
⋮
e
¯
N
k
+
1
=
1
-
β
¯
1
⋮
⋱
⋱
0
⋯
-
β
¯
1
-
1
0
α
¯
0
⋮
⋱
⋱
0
⋯
α
¯
0
e
¯
1
k
e
¯
2
k
⋮
e
¯
N
k
+
h
O
h
k
O
h
k
⋮
O
h
k
.The fine time step h can be regarded as the convergence factor of deferred correction iterations. If rth order numerical integrator is employed in deferred correction, this factor will be obviously h
r. The constant α
¯, which describes the difference between the true-solution propagator and the coarse propagator, can be regarded as a part of the convergence factor of parareal in contrast to deferred correction. We can choose h close enough to α
¯ to balance the two iterations in the PDC method.Note that the fine time steph is crucial for the final approximation performance of the PDC method. We denote the number of the fine time steps on each time subdomain I
n by s, where s
=
Δ
T
/
h. Polynomials with degree of s are employed in the process of deferred correction, and s is the upper bound of the final approximate order of the PDC method. So there should be a moderate integer r
1 such that s
≥
r
1, which means the final result could reach at least the accuracy O
(
h
r
1
). However, the s can not be very large either because the polynomials of low degree are preferred for interpolations. So there should be a little bigger integer r
2 such that s
≤
r
2 to avoid Runge’s phenomenon.Now, we get a value interval for the fine time step, withΔ
T
/
r
2
≤
h
≤
Δ
T
/
r
1. The choice of the optimal time step is(51)
h
=
Δ
T
r
2
,
α
¯
≤
Δ
T
r
2
,
α
¯
,
Δ
T
r
1
≤
α
¯
≤
Δ
T
r
2
,
Δ
T
r
1
,
α
¯
≥
Δ
T
r
1
,which means that h should be chosen as the nearest point to α
¯ from the value interval [
Δ
T
/
r
2
,
Δ
T
/
r
1
].


## 7. System Efficiency Analysis

Parareal can be implemented in two parallel mechanisms which are called serial-parallel parareal and pipelined parareal in [34]. For the pipelined parareal, each processor alternatively compute the fine and coarse propagators without waiting except the initial wait. It had been shown that the pipelined fashion could be better in terms of the computation speed. In this section, we study the parallel speedup and efficiency of the PDC method based on pipelined parareal method.We first investigate the PDC method for ODEs system (12) from optimal one-step QSC method for PDE. Here we suppose backward Euler methods are used in both the coarse and fine propagators. We also assume that the computing costs of the algebraic equations involved in G and F propagators are the same if they have the same scale and sparse structure. When analyzing the speedup of the classical parareal method, we equalize the computing costs of one time step at each iteration. In this paper we do not give such assumption any more, since interpolations are necessary in correction steps of the PDC method, which is not included in the prediction step. The costs of information transmission and data storage are ignored in this paper.Similar to the speedup and system efficiency analysis shown in [34], we adopt the following notation:(i)
T is the length of time domain.
(ii)
Δ
T is the length of each time subdomain.
(iii)
Δ
t is the time increment for the coarse propagator.
(iv)
δ
t is the time increment for the fine propagator.
(v)
τ
G is the cost of the numerical method at each time step in G.
(vi)
N
G is the number of time steps of the coarse propagator on each time subdomain, and we often choose N
G
=
1.
(vii)
τ
F is the cost of the numerical method at each time step in F.
(viii)
N
F is the number of time steps of the fine propagator on each time subdomain.
(ix)
τ
I is the cost of interpolation averaged at each time step in correction steps.
(x)
N is the number of time subdomains on [
0
,
T
]; that is, T
=
N
Δ
T.
(xi)
K is the number of PDC iterations.
(xii)
S is the speedup.
(xiii)
E is the system efficiency and equals the ratio E
=
S
/
N.According to the pipelined parallel mechanism in [34], the initial step of the PDC method is done with the G propagator and the cost for N steps is N
N
G
τ
G. The cost of the first iteration of the PDC method is N
G
τ
G
+
N
F
τ
F, while the cost of any other iteration of the PDC method is N
G
τ
G
+
N
F
(
τ
F
+
τ
I
).The total cost for the initial andK iterations of the PDC method is (52)
K
+
N
N
G
τ
G
+
K
N
F
τ
F
+
K
-
1
N
F
τ
I
.Let M be the number of deferred correction iterations to compute the same problem sequentially on the fine DC nodes with one processor to achieve the same grade accuracy as K PDC iterations. Then the cost of the serial DC method is N
(
M
N
F
τ
F
+
(
M
-
1
)
N
F
τ
I
). The parallel speedup S of the PDC method is (53)
S
=
N
M
N
F
τ
F
+
M
-
1
N
F
τ
I
K
+
N
N
G
τ
G
+
K
N
F
τ
F
+
K
-
1
N
F
τ
I
,and the system efficiency(54)
E
=
S
N
=
M
N
F
τ
F
+
M
-
1
N
F
τ
I
K
+
N
N
G
τ
G
+
K
N
F
τ
F
+
K
-
1
N
F
τ
I
.If we further assume τ
G
=
τ
F and τ
I
=
β
τ
F for simplicity and let α
=
N
G
/
N
F, the efficiency will be (55)
E
=
M
+
M
-
1
β
K
+
N
α
+
K
+
K
-
1
β
=
M
-
β
/
1
+
β
K
+
N
α
/
1
+
β
+
K
-
β
/
1
+
β
,where β denotes the ratio between the cost for solving algebraic equation and the cost for interpolations; it is often a small constant according to the numerical experiment later.It is reported in [34] that the parallel efficiency of the classical pipelined parareal under our simplifying assumptions is (56)
1
K
+
N
α
+
K
,which is automatically bounded by 1
/
K. In contrast, the efficiency of PDC method is bounded by M
/
K, which is apparently better than that of the classical parareal method. However, we have to pay additional attention to the constant α. The constant α denotes the ratio between the time steps of the coarse and the fine propagators, since we assume τ
G
=
τ
F. For the classical parareal method, it is required that N
F is much greater than N
G, α for classical parareal method is much less than 1, and the resulting efficiency is very close to its upper bound. For the PDC method, the accuracy is improved by iterative corrections, and N
F is not needed to be much greater than N
G; then α for PDC method is not so small. Therefore, the resulting efficiency of PDC method is not so close to its upper bound as the efficiency of the classical parareal method.In order to show the contribution of QSC method to the whole algorithm, we investigate the PDC method for the following ODEs system generated after the classical finite difference for PDE(57)
d
u
t
d
t
=
1
Δ
x
2
A
u
t
+
f
→
t
,
t
∈
0
,
T
u
0
=
u
0
,where Δ
x is the size of spatial grid, the coefficient matrix (58)
A
=
-
2
1
0
⋯
0
1
-
2
1
⋯
0
⋱
⋱
⋱
0
⋯
1
-
2
1
0
⋯
0
1
-
2
,and f
→
(
t
) is the interpolation vector of f
(
x
,
t
) at equispaced grid points.Letτ
¯
G and τ
¯
F be the cost of the numerical method at each time step in G and F, respectively. Let τ
¯
I be the cost of the interpolation averaged at each time step and define β
¯
=
τ
I
¯
/
τ
F
¯. Obviously, τ
¯
G and τ
¯
F are directly related to the scale and structure of the algebraic equations involved. As we know, the finite difference method has the approximate order O
(
Δ
x
2
), while the optimal QSC method has the order O
(
Δ
x
4
). To achieve the same level of error, the dimension of system (57) should be much greater than the dimension of system (12). Let n
1 and n
2 be the dimensions of system (12) and system (57), respectively. From (
b
-
a
/
n
1
)
4
~
(
b
-
a
/
n
2
)
2, we obtain n
2
~
n
1
2.For the PDC method for system (12), the coefficient matrices of the linear equations have sparse structure with bandwidth 4, and the corresponding time complexity is O
(
4
n
1
). For the PDC method for system (57), the coefficient matrices of the linear equations have sparse structure with bandwidth 1, and the corresponding time complexity is O
(
n
2
). Then we have approximate relationships τ
¯
G
≈
n
1
/
4
τ
G and τ
¯
F
≈
n
1
/
4
τ
F, and the approximate system efficiency is (59)
E
¯
=
M
-
β
¯
/
1
+
β
¯
K
+
N
α
/
1
+
β
¯
+
K
-
β
¯
/
1
+
β
¯
,which is almost the same as E, except for a different constant β
¯. Therefore, the application of the QSC method nearly has no effect to increase the system efficiency.The total cost afterK iterations of the PDC method for system (57) is (60)
K
+
N
N
G
τ
¯
G
+
K
N
F
τ
¯
F
+
K
-
1
N
F
τ
¯
I
≈
n
1
4
K
+
N
N
G
τ
G
+
K
N
F
τ
F
+
K
-
1
N
F
β
¯
β
τ
I
,where the ratio β
¯
/
β is a small constant. The cost of PDC method for system (57) is about n
1
/
4 times of the cost of the PDC method for system (12). Since n
1 is a very big integer, the QSC method reduces the total cost greatly. Furthermore, if the optimal two-step QSC method is employed, the reduction factor will be n
1
/
2 rather than n
1
/
4.


## 8. Numerical Experiments

In this section, we take the following parabolic PDE as an example to present the advantages of the QSC-PDC algorithm:(61)
∂
∂
t
u
x
,
t
=
∂
2
∂
x
2
u
x
,
t
+
2
π
2
-
1
e
-
t
/
2
sin
⁡
π
x
,
u
0
,
t
=
u
1
,
t
=
0
,
0
≤
t
≤
5
,
u
x
,
0
=
2
sin
⁡
π
x
,
0
≤
x
≤
1
.The true solution of system (61) is (62)
u
x
,
t
=
2
e
-
t
/
2
sin
⁡
π
x
,
0
≤
x
≤
1
,
0
≤
t
≤
5
.To perform the QSC-PDC method, we first employ a uniform partition of the spatial domain[
0,1
] with 160 equispaced grid points, and the step size is Δ
x
=
1
/
160. Then we take the midpoint of each segment as the collocation points; that is, (63)
τ
i
=
i
-
1
2
Δ
x
,
i
=
1,2
,
…
,
160
.Define the vector g
→
=
(
sin
⁡
π
τ
1
,
sin
⁡
π
τ
2
,
…
,
sin
⁡
(
π
τ
160
)
)
T. The resulting collocation equation is(64)
Q
0
d
c
t
d
t
=
1
Δ
x
2
Q
2
+
1
24
Q
x
x
c
t
+
2
π
2
-
1
e
-
t
/
2
g
→
,
t
∈
0,5
,
c
0
=
2
Q
0
g
→
,where the matrices Q
0, Q
2, and Q
x
x
∈
R
160
×
160 are the same as the matrices in system (12). Here, we only consider the optimal one-step QSC method.We divide the time domain[
0,5
] into 40 subdomains. On each subdomain, we use 1-step backward Euler method for the G propagators and use 5-step backward Euler method for the F
D
C propagators. Polynomials of degree 5 are employed for interpolation in each correction loop. The PDC method is programmed in the pipelined manner, and information is communicated between processors through MPI, and the code is written using C language with GNU Scientific Library. The program is carried out on the Shared Hierarchical Academic Research Computing Network of Canada (https://www.sharcnet.ca/). 41 CPUs are employed in the experiments, including one root CPU and the other 40 CPUs. The root CPU is in charge of sending and receiving data and displaying outputs. The iterative computing of the propagators on the 40 time subdomains are distributed onto the other 40 CPUs.For simplicity, we only examine the error at a special point(
x
,
t
)
=
(
(
80
-
1
/
2
)
Δ
x
,
5
). The resulting PDC error and corresponding running time are shown in Table 2. In order to exhibit the advantages of the QSC-PDC method, we perform the classical parareal method for system (64), and the spatial variable is also discretized by optimal QSC. Here, we call the reference method as QSC-parareal method. The time domain [
0,5
] is also divided into 40 subdomains, and the coarse propagator is also chosen as 1-step backward Euler method on each subdomain. The fine propagator is chosen as 100-step 2-stage Gauss method on each subdomain, which is of high order and easy to carry out.Table 2
The iteration error and running time of QSC-PDC method, withN
=
40, s
=
5.
Iterations
k
=
1
k
=
2
k
=
3
k
=
4
k
=
5
k
=
6
k
=
7
k
=
8
k
=
9
k
=
10
Error
2.05
e
-
6
1.65
e
-
5
1.57
e
-
6
3.73
e
-
7
1.38
e
-
7
3.25
e
-
9
4.89
e
-
9
7.71
e
-
10
1.50
e
-
10
5.71
e
-
11
Time (s)
0.47
0.68
0.87
1.10
1.30
1.50
1.73
1.92
2.11
2.32The classical parareal method for (64) is also programmed in pipelined manner, and the corresponding error at the same point (
x
,
t
)
=
(
(
80
-
1
/
2
)
Δ
x
,
5
) and running time are shown in Table 3.Table 3
The iteration error and running time of classical method for the QSC equations, with Gauss methods ofN
=
40 and 200 steps on each subdomain.
Iterations
k
=
1
k
=
2
k
=
3
k
=
4
k
=
5
k
=
6
k
=
7
k
=
8
k
=
9
k
=
10
Error
8.91
e
-
5
2.83
e
-
5
9.04
e
-
6
2.86
e
-
6
9.29
e
-
7
2.75
e
-
7
1.09
e
-
7
1.24
e
-
8
2.67
e
-
8
1.45
e
-
8
Time (s)
1.24
2.13
3.11
4.10
5.09
6.02
6.97
7.91
8.95
9.83To show the advantages of the deferred correction clearly, we plot the iteration error and running time of the two methods together in Figure6. We can see that the deferred correction technique makes the parareal method achieve better approximations in fixed iterations and costs much less running time.Figure 6
The comparison of error and cost between QSC-PDC method and QSC-parareal method for system (61).


## 9. Conclusions

In this paper, we firstly applied the optimal one-step QSC method to solve linear parabolic PDEs, which lead to linear ODEs with respect to the degrees of freedom. We then combined the parareal algorithm with deferred correction method along the temporal direction to deal with the collocation equations. With the fourth order error on the spatial step size, high accuracy could be obtained without too many collocation points. We also supplied an error estimation of the PDC methods, which illustrates good convergence with small time subdomains. We also analyzed the stabilities of the PDC method with different cases of coarse and fine propagators. We presented the tuning strategy for the time step size of the deferred correction method to balance the parareal and the deferred correction iterations of distinguish individual computational speeds.Compared with the classical parareal method, PDC method has much higher parallel efficiency, and the optimal QSC can reduce the computational cost greatly when achieving a desired accuracy. By the numerical experiments on multicore computers, the advantages of the QSC-PDC algorithm are obvious.



---

*Source: 6943079-2016-05-18.xml*