# Parareal Algorithms Implemented with IMEX Runge-Kutta Methods

Authors: Zhiyong Wang; Shulin Wu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Parareal algorithm is a very powerful parallel computation method and has received much interest from many researchers over the past few years. The aim of this paper is to investigate the performance of parareal algorithm implemented with IMEX Runge-Kutta (RK) methods. A stability criterion of the parareal algorithm coupled with IMEX RK methods is established and the advantage (in the sense of stability) of implementing with this kind of RK methods is numerically investigated. Finally, numerical examples are given to illustrate the efficiency and performance of different parareal methods.


---

## Body

## 1. Introduction

The parareal algorithm was presented by Lions et al. in [1] as a numerical method to solve time dependent problems in parallel. The peculiarity of this algorithm is that it approximates successfully the solution later in time before having fully accurate approximations from earlier time points. This algorithm is becoming more and more popular in scientific and engineering computation and many excellent results have been obtained.As mentioned above, the parareal algorithm was first introduced in [1] and an improved version which aims to solve nondifferentiable PDEs was discussed by Bal and Maday in [2]. Some further modifications and improvements can be found in [3], in which the authors use the “filtering” technique to overcome the so-calledresonance orbeating phenomenon (see also [4, 5]) that triggers numerical instability when the algorithm is applied to linear structure dynamics. The stability and convergence are the main subjects of theoretical analysis of the algorithm which have been investigated widely by many researchers; see, for example, [6–8]. Nowadays, this algorithm, as well as its variants [3, 9–12], has been used in many fields by many researchers, such as morphological transformation simulations [13], structural (fluid) dynamics simulations [4, 5, 14], optimal control [15, 16], Hamiltonian simulations [17, 18], turbulent plasma simulations [19, 20], and solution of Volterra integral equations [21]. (There is an increasing interest in parareal, and it is very possible that some important references are not mentioned here.)For ODEs,(1)
y
′
(
t
)
=
f
t
,
y
t
,
t
∈
[
0
,
T
]
,
y
(
0
)
=
y
0
,to formulate the parareal algorithm, the time domain [
0
,
T
] of interest is first partitioned into several time-slices whose boundary points are treated as coarse time-grids. And then, the algorithm consists of three basic steps. First, using some numerical propagator, say G
Δ
T, the solution is approximated on each coarse time-grid to provide a seed—that is, an initial condition—to the time-slice. Second, another propagator, denoted by F
Δ
t, is applied independently and therefore concurrently in each time-slice to advance the solution from the starting point of this time-slice to its end point. Finally, an iterative process is invoked to improve the accuracy of the seeds and eliminate the jumps of the solution on the coarse time-grids. In most cases, the numerical propagators G
Δ
T and F
Δ
t are defined by traditional RK method with equal length of the time-slices, such as Radau methods, Lobatto methods, and Gauss methods, and there are lots of investigations about theoretical properties and performance of practical application of the algorithm in such case. See, for example, [4, 5, 7, 9, 14] and references therein.There is also some research about special choice of the underlying numerical methods that are used to defineG
Δ
T and F
Δ
t. For example, Guibert and Tromeur-Dervout [22] considered the adaptive time-slices (i.e., the length of time-slices is not equal) to define G
Δ
T and F
Δ
t and the derived parareal algorithm can be applied to strong stiff ODEs and differential algebraic equations (DAEs). Another example is the one investigated by Minion [11], where the authors defined the numerical propagators G
Δ
T and F
Δ
t by spectral deferred correction based on the Gaussian quadrature spectral integration. It was shown that the iterative process based on spectral deferred correction rather than traditional methods within a parareal framework results in a significant decrease in the overall computational cost of the algorithm. Wu et al. [9] suggested that the propagator F
Δ
t was defined by local Richardson extrapolation, and it turns out that the Parareal-Richardson method has advantages of higher accuracy and better stability.In some cases, the functionf in (1) can be written into two parts (2)
f
=
f
s
+
f
n
s
,where f
n
s is the nonstiff or mildly stiff part of f and f
s is the stiff part. A common example for such case arises from semidiscretization of reaction diffusion equations u
t
=
d
u
x
x
+
r
(
t
,
x
,
u
). In this case the function f in (1) can be written as f
(
t
,
y
(
t
)
)
=
A
y
(
t
)
+
f
n
s
(
t
,
y
(
t
)
) with A being a tridiagonal matrix and stiff. In the case that f can be partitioned into stiff and nonstiff parts, it is more advisable to use the implicit-explicit (IMEX) RK methods (or additive RK methods called sometimes) to solve ODEs (1). For this type of RK methods, the stiff part is assigned with the implicit component of the IMEX RK method to satisfy stability requirement and the nonstiff or mildly stiff part is assigned with the explicit component to reduce computational cost. For more details about the RK methods and the IMEX RK methods, the interested reader may refer to [23–33] and references therein.Therefore, under condition (2), we think that it is valuable to adopt the IMEX RK methods instead of the traditional RK methods to define G
Δ
T and/or F
Δ
t used in parareal framework. This is the main motivation of our work. In certain aspects, this paper should be viewed as an experimental one. While the stability of the derived algorithm we present is easily proven by the results given by Gander and Vandewalle [7], stability region of the parareal algorithm defined by the IMEX RK methods is numerically shown.The structure of this paper is as follows. In Section2, the parareal algorithm and the IMEX RK methods are revisited. Section 3 is concerned with the stability of the derived parareal algorithm. Several combinations of different IMEX RK methods and the corresponding advantages/disadvantages are numerically shown. Finally, in Section 4, we test the Gray-Scott model arising in chemical reaction to illustrate our results.


## 2. The Parareal Algorithm and the IMEX RK Methods

In this section, we revisit the parareal algorithm and the IMEX RK methods.


### 2.1. The Parareal Algorithm

To introduce the parareal algorithm, let us first partition time interval[
0
,
T
] into N subintervals S
n
=
[
T
n
,
T
n
+
1
], n
=
0,1
,
…
,
N
-
1, with equal length Δ
T and T
n
=
n
Δ
T. We call T
n’s the coarse time-grids. We then use some finer step size Δ
t
=
Δ
T
/
M (M
>
1 is an integer) to partition each relative large interval S
n into M finer intervals s
n
,
m
=
T
n
+
m
/
M
,
T
n
+
(
m
+
1
)
/
M with T
n
+
m
/
M
=
T
n
+
m
Δ
t and m
=
0,1
,
…
,
M
-
1. Now, the parareal algorithm can be described as follows. We designate by symbol ⊖ the time-sequential steps performed on the coarse time-grids and by symbol ⊕ the parallel steps performed on the decomposed finer time-grids.The Parareal Algorithm. Consider the following.⊖
Step
0 (initialization). Perform sequential computation Y
n
+
1
0
=
G
Δ
T
(
T
n
,
Y
n
0
,
Δ
T
) with Y
0
0
=
y
0, n
=
0,1
,
…
,
N
-
1, for k
=
0,1
,
….⊕
Step
1. Perform in subinterval S
n
=
[
T
n
,
T
n
+
1
] the computation Y
^
n
+
(
m
+
1
)
/
M
=
F
Δ
t
T
n
+
m
/
M
,
Y
^
n
+
m
/
M
,
Δ
t with initial value Y
n
k, m
=
0,1
,
…
,
M
-
1.⊖
Step
2. Perform sequential corrections Y
n
+
1
k
+
1
=
G
Δ
T
T
n
,
Y
n
k
+
1
,
Δ
T
+
Y
^
n
+
1
-
G
Δ
T
T
n
,
Y
n
k
,
Δ
T with Y
0
k
+
1
=
y
0, n
=
0,1
,
…
,
N
-
1.⊖
Step
3. If for n
=
0,1
,
…
,
N
-
1, Y
n
k
+
1 satisfy some termination criteria, stop the iteration; else go to Step 1.Here, the finer propagation operatorF
Δ
t and the coarse propagation operator G
Δ
T relay on one-step numerical methods and if the underlying methods are implicit, in most cases (linear ODEs is an exception), nonlinear algebraic equations need to be solved.The above algorithm can be written compactly as(3)
Y
n
+
1
k
+
1
=
G
Δ
T
T
n
,
Y
n
k
+
1
,
Δ
T
+
F
Δ
t
T
n
,
Y
n
k
,
Δ
t
-
G
Δ
T
T
n
,
Y
n
k
,
Δ
T
,where k is the iteration index. Note that for k
→
+
∞ method (3) will generate a series of values Y
n which satisfy Y
n
+
1
=
F
Δ
t
T
n
,
Y
n
,
Δ
t. This implies that the converged solution Y
n will achieve the accuracy of the propagator F
Δ
t.


### 2.2. The IMEX RK Methods

Let us consider a pair of two Runge-Kutta methods defined by the arrays(4)
0
0
0
0
⋯
0
c
2
a
21
a
22
0
⋯
0
c
3
a
31
a
32
a
33
0
⋮
⋮
⋮
⋮
⋮
⋱
0
c
s
a
s
1
a
s
2
⋯
a
s
,
s
-
1
a
s
s
b
1
b
2
⋯
b
s
-
1
b
s
0
0
0
0
⋯
0
c
2
a
^
21
0
0
⋯
0
c
3
a
^
31
a
^
32
0
0
⋮
⋮
⋮
⋮
⋮
⋱
0
c
s
a
^
s
1
a
^
s
2
⋯
a
^
s
,
s
-
1
0
3
r
h
t
b
^
1
b
^
2
⋯
b
^
s
-
1
b
^
swith the same abscissae (5)
c
i
=
∑
j
=
1
i
‍
a
i
j
=
∑
j
=
1
i
-
1
‍
a
^
i
j
,
i
=
2
,
…
,
s
.The top formula of (4) determines a diagonally implicit (semi-implicit) RK method and the bottom formula is an explicit RK method. In addition, let h
>
0 be a step size and define the step point t
n
=
n
h for integer n. Consider the following stiff/nonstiff partitioned ODEs: (6)
y
′
(
t
)
=
f
s
(
t
,
y
(
t
)
)
+
f
n
s
(
t
,
y
(
t
)
)
,
t
∈
[
0
,
T
]
,
y
(
0
)
=
y
0
,and by applying the top part of (4) to stiff component f
s and the bottom part to the nonstiff component, we obtain the following scheme: (7)
K
i
=
y
n
+
h
∑
j
=
1
i
‍
a
i
j
f
s
t
n
+
c
j
h
,
K
j
+
h
∑
j
=
1
i
-
1
‍
a
^
i
j
f
n
s
t
n
+
c
j
h
,
K
j
,
i
=
1,2
,
…
,
s
,
y
n
+
1
=
y
n
+
h
∑
j
=
1
s
‍
b
j
f
s
t
n
+
c
j
h
,
K
j
+
h
∑
j
=
1
s
‍
b
^
j
f
n
s
t
n
+
c
j
h
,
K
j
,where ∑
j
=
1
i
‍
·
=
0 with i
<
j.The order conditions for IMEX RK methods are very complicated and beyond the topic of this paper. For this aspect, we refer the interested reader to, for example, [26, 30]. At the moment, we introduce some frequently used IMEX methods with order p
=
1,2
,
3. We will use the triplet IMEX
(
s
im
,
s
ex
,
p
) to identify scheme (4), where s
im is the number of stages of the implicit part, s
ex is the number of stages of the explicit part, and p is the order of the IMEX scheme.First Order IMEX RK Methods. The most popular first order IMEX RK method is the IMEX Euler method: (8)
0
0
0
1
0
1
0
1
0
0
0
1
1
0
1
0
︸
IMEX
Euler
Second Order IMEX RK Methods. The second order IMEX RK methods considered here are the IMEX trapezoidal scheme [34]: (9)
0
0
0
1
1
2
1
2
1
0
1
2
1
1
0
2
0
0
0
1
1
0
1
2
1
1
0
2
︸
IMEX
trapezoidalThe IMEX trapezoidal scheme is a combination of the trapezoidal rule and Heun’s second order method (the explicit trapezoidal rule).Third Order IMEX RK Methods. For third order IMEX RK method, we consider the one constructed by Ascher et al. [23] (this scheme is denoted by IMEX
(
4,4
,
3
)):(10)
0
0
0
0
0
a
0
a
0
0
0.7179332607
0
0.2820667392
a
0
1
0
1.2084966490
-
0.644363171
a
0
1.2084966490
-
0.644363171
a
0
0
0
0
a
0
0
0
0.3212788860
0.3966543747
0
0
-
0.1058858296
b
b
0
0
1.2084966490
-
0.644363171
a
︸
IMEX
(
4,4
,
3
)where a
=
0.4358665215 and b
=
0.552929147. Both the implicit and explicit schemes of IMEX
(
4,4
,
3
) are third order RK methods, and the implicit scheme is L-stable.To finish this section, we introduce the concept of stability of the IMEX RK methods which will be used in the stability analysis of the parareal algorithm.The asymptotic behavior of the IMEX RK methods applied to (6) is analyzed on the basis of the scalar model equation: (11)
y
′
(
t
)
=
λ
y
+
μ
y
,
λ
,
μ
∈
C
,which was proposed by Frank et al. [35] and Verwer and Sommeijer [36] and adopted by Koto [28, 37]. For model equation (11), from (7) we arrive at (12)
K
=
E
y
n
+
h
λ
A
K
+
h
μ
A
^
K
,
y
n
+
1
=
y
n
+
h
λ
b
K
+
h
μ
b
^
K
,where we have used vector and matrix notations to denote the values shown in the Butcher tableau (4); that is, A
=
(
a
i
j
), A
^
=
(
A
^
i
j
), b
=
(
b
1
,
b
2
,
…
,
b
s
), and b
^
=
(
b
^
1
,
b
^
2
,
…
,
b
^
s
). Moreover, E
=
(
1,1
,
…
,
1
)
T.Letα
=
h
λ, β
=
h
μ. We then have (13)
y
n
+
1
=
R
(
α
,
β
)
y
n
,where R
(
α
,
β
) is calledstability function of the IMEX RK method which is defined by (14)
R
(
α
,
β
)
=
1
+
α
b
+
β
b
^
I
-
α
A
-
β
A
^
-
1
E
,where I is an identity matrix of size s.


## 2.1. The Parareal Algorithm

To introduce the parareal algorithm, let us first partition time interval[
0
,
T
] into N subintervals S
n
=
[
T
n
,
T
n
+
1
], n
=
0,1
,
…
,
N
-
1, with equal length Δ
T and T
n
=
n
Δ
T. We call T
n’s the coarse time-grids. We then use some finer step size Δ
t
=
Δ
T
/
M (M
>
1 is an integer) to partition each relative large interval S
n into M finer intervals s
n
,
m
=
T
n
+
m
/
M
,
T
n
+
(
m
+
1
)
/
M with T
n
+
m
/
M
=
T
n
+
m
Δ
t and m
=
0,1
,
…
,
M
-
1. Now, the parareal algorithm can be described as follows. We designate by symbol ⊖ the time-sequential steps performed on the coarse time-grids and by symbol ⊕ the parallel steps performed on the decomposed finer time-grids.The Parareal Algorithm. Consider the following.⊖
Step
0 (initialization). Perform sequential computation Y
n
+
1
0
=
G
Δ
T
(
T
n
,
Y
n
0
,
Δ
T
) with Y
0
0
=
y
0, n
=
0,1
,
…
,
N
-
1, for k
=
0,1
,
….⊕
Step
1. Perform in subinterval S
n
=
[
T
n
,
T
n
+
1
] the computation Y
^
n
+
(
m
+
1
)
/
M
=
F
Δ
t
T
n
+
m
/
M
,
Y
^
n
+
m
/
M
,
Δ
t with initial value Y
n
k, m
=
0,1
,
…
,
M
-
1.⊖
Step
2. Perform sequential corrections Y
n
+
1
k
+
1
=
G
Δ
T
T
n
,
Y
n
k
+
1
,
Δ
T
+
Y
^
n
+
1
-
G
Δ
T
T
n
,
Y
n
k
,
Δ
T with Y
0
k
+
1
=
y
0, n
=
0,1
,
…
,
N
-
1.⊖
Step
3. If for n
=
0,1
,
…
,
N
-
1, Y
n
k
+
1 satisfy some termination criteria, stop the iteration; else go to Step 1.Here, the finer propagation operatorF
Δ
t and the coarse propagation operator G
Δ
T relay on one-step numerical methods and if the underlying methods are implicit, in most cases (linear ODEs is an exception), nonlinear algebraic equations need to be solved.The above algorithm can be written compactly as(3)
Y
n
+
1
k
+
1
=
G
Δ
T
T
n
,
Y
n
k
+
1
,
Δ
T
+
F
Δ
t
T
n
,
Y
n
k
,
Δ
t
-
G
Δ
T
T
n
,
Y
n
k
,
Δ
T
,where k is the iteration index. Note that for k
→
+
∞ method (3) will generate a series of values Y
n which satisfy Y
n
+
1
=
F
Δ
t
T
n
,
Y
n
,
Δ
t. This implies that the converged solution Y
n will achieve the accuracy of the propagator F
Δ
t.


## 2.2. The IMEX RK Methods

Let us consider a pair of two Runge-Kutta methods defined by the arrays(4)
0
0
0
0
⋯
0
c
2
a
21
a
22
0
⋯
0
c
3
a
31
a
32
a
33
0
⋮
⋮
⋮
⋮
⋮
⋱
0
c
s
a
s
1
a
s
2
⋯
a
s
,
s
-
1
a
s
s
b
1
b
2
⋯
b
s
-
1
b
s
0
0
0
0
⋯
0
c
2
a
^
21
0
0
⋯
0
c
3
a
^
31
a
^
32
0
0
⋮
⋮
⋮
⋮
⋮
⋱
0
c
s
a
^
s
1
a
^
s
2
⋯
a
^
s
,
s
-
1
0
3
r
h
t
b
^
1
b
^
2
⋯
b
^
s
-
1
b
^
swith the same abscissae (5)
c
i
=
∑
j
=
1
i
‍
a
i
j
=
∑
j
=
1
i
-
1
‍
a
^
i
j
,
i
=
2
,
…
,
s
.The top formula of (4) determines a diagonally implicit (semi-implicit) RK method and the bottom formula is an explicit RK method. In addition, let h
>
0 be a step size and define the step point t
n
=
n
h for integer n. Consider the following stiff/nonstiff partitioned ODEs: (6)
y
′
(
t
)
=
f
s
(
t
,
y
(
t
)
)
+
f
n
s
(
t
,
y
(
t
)
)
,
t
∈
[
0
,
T
]
,
y
(
0
)
=
y
0
,and by applying the top part of (4) to stiff component f
s and the bottom part to the nonstiff component, we obtain the following scheme: (7)
K
i
=
y
n
+
h
∑
j
=
1
i
‍
a
i
j
f
s
t
n
+
c
j
h
,
K
j
+
h
∑
j
=
1
i
-
1
‍
a
^
i
j
f
n
s
t
n
+
c
j
h
,
K
j
,
i
=
1,2
,
…
,
s
,
y
n
+
1
=
y
n
+
h
∑
j
=
1
s
‍
b
j
f
s
t
n
+
c
j
h
,
K
j
+
h
∑
j
=
1
s
‍
b
^
j
f
n
s
t
n
+
c
j
h
,
K
j
,where ∑
j
=
1
i
‍
·
=
0 with i
<
j.The order conditions for IMEX RK methods are very complicated and beyond the topic of this paper. For this aspect, we refer the interested reader to, for example, [26, 30]. At the moment, we introduce some frequently used IMEX methods with order p
=
1,2
,
3. We will use the triplet IMEX
(
s
im
,
s
ex
,
p
) to identify scheme (4), where s
im is the number of stages of the implicit part, s
ex is the number of stages of the explicit part, and p is the order of the IMEX scheme.First Order IMEX RK Methods. The most popular first order IMEX RK method is the IMEX Euler method: (8)
0
0
0
1
0
1
0
1
0
0
0
1
1
0
1
0
︸
IMEX
Euler
Second Order IMEX RK Methods. The second order IMEX RK methods considered here are the IMEX trapezoidal scheme [34]: (9)
0
0
0
1
1
2
1
2
1
0
1
2
1
1
0
2
0
0
0
1
1
0
1
2
1
1
0
2
︸
IMEX
trapezoidalThe IMEX trapezoidal scheme is a combination of the trapezoidal rule and Heun’s second order method (the explicit trapezoidal rule).Third Order IMEX RK Methods. For third order IMEX RK method, we consider the one constructed by Ascher et al. [23] (this scheme is denoted by IMEX
(
4,4
,
3
)):(10)
0
0
0
0
0
a
0
a
0
0
0.7179332607
0
0.2820667392
a
0
1
0
1.2084966490
-
0.644363171
a
0
1.2084966490
-
0.644363171
a
0
0
0
0
a
0
0
0
0.3212788860
0.3966543747
0
0
-
0.1058858296
b
b
0
0
1.2084966490
-
0.644363171
a
︸
IMEX
(
4,4
,
3
)where a
=
0.4358665215 and b
=
0.552929147. Both the implicit and explicit schemes of IMEX
(
4,4
,
3
) are third order RK methods, and the implicit scheme is L-stable.To finish this section, we introduce the concept of stability of the IMEX RK methods which will be used in the stability analysis of the parareal algorithm.The asymptotic behavior of the IMEX RK methods applied to (6) is analyzed on the basis of the scalar model equation: (11)
y
′
(
t
)
=
λ
y
+
μ
y
,
λ
,
μ
∈
C
,which was proposed by Frank et al. [35] and Verwer and Sommeijer [36] and adopted by Koto [28, 37]. For model equation (11), from (7) we arrive at (12)
K
=
E
y
n
+
h
λ
A
K
+
h
μ
A
^
K
,
y
n
+
1
=
y
n
+
h
λ
b
K
+
h
μ
b
^
K
,where we have used vector and matrix notations to denote the values shown in the Butcher tableau (4); that is, A
=
(
a
i
j
), A
^
=
(
A
^
i
j
), b
=
(
b
1
,
b
2
,
…
,
b
s
), and b
^
=
(
b
^
1
,
b
^
2
,
…
,
b
^
s
). Moreover, E
=
(
1,1
,
…
,
1
)
T.Letα
=
h
λ, β
=
h
μ. We then have (13)
y
n
+
1
=
R
(
α
,
β
)
y
n
,where R
(
α
,
β
) is calledstability function of the IMEX RK method which is defined by (14)
R
(
α
,
β
)
=
1
+
α
b
+
β
b
^
I
-
α
A
-
β
A
^
-
1
E
,where I is an identity matrix of size s.


## 3. Stability of the Parareal Algorithm

In our analysis of the stability of the parareal algorithm implemented with IMEX RK methods, an important role is a strictly lower triangular Toeplitz matrixM
:
=
M
(
σ
) of size N. Its elements are defined by (15)
M
i
j
=
0
,
if
i
≤
j
,
σ
i
-
j
-
1
,
if
i
>
j
.The next necessity for our analysis is the estimation of M
k
∞, which is proved by Gander and Vandewalle [7].Lemma 1 (see [7]).
Letk
>
0 be an integer. Then the infinite norm of M
k can be estimated by (16)
M
k
σ
∞
≤
min
⁡
1
-
σ
N
-
1
1
-
σ
k
,
N
-
1
k
,
i
f
σ
<
1
,
σ
N
-
k
-
1
N
-
1
k
,
i
f
σ
≥
1
.Assume that the stability functions of the underlying numerical methods used to define the two propagatorsG
Δ
T and F
Δ
t are R
F
(
α
,
β
) and R
G
(
α
,
β
), respectively. Then it follows by applying the parareal algorithm to the model equation (11) that (17)
G
Δ
T
T
n
,
Y
n
k
,
Δ
T
=
R
G
(
α
,
β
)
Y
n
k
,
F
Δ
t
T
n
,
Y
n
k
,
Δ
T
=
R
F
M
α
M
,
β
M
Y
n
k
,
k
=
1,2
,
…
,where we recall M
=
Δ
T
/
Δ
t.Theorem 2.
LetΔ
T be given and T
n
=
n
Δ
T for n
=
0,1
,
….  Let the underlying numerical method used to define the propagator G
Δ
T running on the coarse time-grids be in its region of absolute stability; that is, |
R
G
|
<
1. Then one has the following estimates: (18)
sup
⁡
n
>
0
y
t
n
-
Y
n
k
≤
R
F
M
α
/
M
,
β
/
M
-
R
G
α
,
β
1
-
R
G
α
,
β
k
sup
⁡
n
>
0
y
t
n
-
Y
n
0
,
m
m
m
m
m
m
m
m
m
m
m
m
m
i
m
m
m
m
m
m
m
m
m
m
m
m
m
k
≥
1
.Proof.
We denote bye
n
k the error at iteration k of the parareal algorithm at coarse time point T
n; that is, e
n
k
=
y
(
t
n
)
-
Y
n
k. Then with an induction argument on n, it follows by applying the iterative process (3) to the model equation (11) that this error satisfies (19)
e
n
+
1
k
+
1
=
R
G
(
α
,
β
)
e
n
k
+
1
+
R
F
α
M
,
β
M
-
R
G
α
,
β
e
n
k
=
R
F
α
M
,
β
M
-
R
G
α
,
β
∑
j
=
1
n
‍
R
G
n
-
j
(
α
,
β
)
e
j
k
,where we have used the fact that e
0
k
=
0 for any k
≥
0. Set E
k
=
e
1
k
,
e
2
k
,
…
,
e
N
k
T. Then relation (19) can be written in matrix form as (20)
E
k
+
1
=
R
F
α
M
,
β
M
-
R
G
(
α
,
β
)
M
R
G
(
α
,
β
)
E
k
,where the matrix M is defined by (15) with σ
=
R
G
(
α
,
β
). This implies (21)
E
k
=
R
F
α
M
,
β
M
-
R
G
(
α
,
β
)
k
M
k
R
G
(
α
,
β
)
E
0
.Therefore, from Lemma 1 and the assumption |
R
G
|
<
1, we have (22)
E
k
∞
≤
R
F
α
M
,
β
M
-
R
G
α
,
β
k
×
M
k
R
G
(
α
,
β
)
E
0
∞
≤
R
F
α
M
,
β
M
-
R
G
α
,
β
k
×
1
-
R
G
α
,
β
N
-
1
1
-
R
G
α
,
β
k
E
0
∞
.Hence, the proof of this theorem is completed by letting N
→
+
∞ in the second inequality of (22).The factorρ
(
α
,
β
)
:
=
|
R
F
M
(
α
/
M
,
β
/
M
)
-
R
G
(
α
,
β
)
|
/
(
1
-
R
G
α
,
β
) is called by Gander and Vandewalle [7] thelinear convergence factor of the parareal algorithm performed on unbounded time intervals. According to the concepts introduced by Bal [6] and Staff and Rønquist [8], it is also called the stability function of the parareal algorithm. In this paper, we use the latter name.Define(23)
D
=
(
α
,
β
)
∣
ρ
(
α
,
β
)
<
1
,
α
,
β
∈
C
.Then if (
α
,
β
)
∈
D, the parareal algorithm is convergent on unbounded time intervals. In what follows of this paper, the set D is called stability region of the parareal algorithm.Remark 3.
Ifb
=
b
^ and A
=
A
^ (i.e., the IMEX RK methods reduce to the traditional RK methods), it follows from (14) that R
(
α
,
β
)
=
R
(
z
)
=
1
+
z
b
(
I
-
z
A
)
-
1
E with z
=
α
+
β. Therefore, if both the propagators G
Δ
T and F
Δ
t are defined by the traditional RK methods, the stability function ρ
(
α
,
β
) can be rewritten as ρ
~
(
z
)
=
|
R
F
M
(
z
/
M
)
-
R
G
(
z
)
|
/
(
1
-
|
R
G
(
z
)
|
), which is the one obtained by Gander and Vandewalle in [7].From Remark3, we see that the underlying numerical method used to define propagator G
Δ
T that is in its absolute stability region is necessary for the stability of the parareal algorithm. In fact, in the previous analysis and application of the parareal algorithm (i.e., both the propagators G
Δ
T and F
Δ
t are defined by the traditional RK methods), it has been proposed by Farhat et al. [4, 5, 14] that the underlying RK method for G
Δ
T should be implicit and the one for the propagator F
Δ
t should be explicit. The reason is obvious—the implicit RK method is used to guarantee stability and the explicit RK method is used to reduce computation cost as much as possible.The greatest advantage of using an IMEX RK method is that the stability of the combined scheme approaches to that of its implicit part while the computational cost can bear comparison with simply using the explicit part. Hence, in the parareal framework, it is naturally desired that the adoption of the IMEX RK methods will provide a far more stable parareal algorithm with computational cost sharply reduced. To see whether our desirability can be carried through, for two IMEX RK methods,IME
X
G and IME
X
F, we consider here M
=
50 and the following four special combinations which leads to four parareal algorithms:P
I
M
-
E
X: G
Δ
T is defined by the implicit part of IME
X
G and F
Δ
t is defined by the explicit part IME
X
F;
P
I
M
-
I
M
E
X:  G
Δ
T is defined by the implicit component of IME
X
G and F
Δ
t is defined by IME
X
F;
P
I
M
E
X
-
E
X:  G
Δ
T is defined by IME
X
G and F
Δ
t is defined by the explicit component of IME
X
F;
P
I
M
E
X
-
I
M
E
X:  G
Δ
T is defined by IME
X
G and F
Δ
t is defined by IME
X
F.For the convenience of plotting the stability region of the parareal algorithm, we consider the case (
α
,
β
)
=
(
x
,
i
y
), x
,
y
∈
R. The stability functions of the above four parareal algorithms are denoted by R
IM-EX
(
x
,
y
), R
IM-IMEX
(
x
,
y
), R
IMEX-EX
(
x
,
y
), and R
IMEX-IMEX
(
x
,
y
), respectively. The stability regions are denoted by D
IM-EX, D
IM-IMEX, D
IMEX-EX, and D
IMEX-IMEX, respectively. We remark that the P
I
M
-
E
X algorithm is a conventional parareal algorithm and we will compare its stability region with the ones of the other three algorithms.We first illustrate the stability regions of the parareal algorithms whenIME
X
G is chosen as the IMEX Euler method (see (8)). For IME
X
F = IMEX Euler and IME
X
F = IMEX Trapezoidal and IME
X
F
=
IMEX
(
4,4
,
3
), the stability region of the algorithms P
I
M
-
E
X, P
I
M
-
I
M
E
X, P
I
M
E
X
-
E
X, and P
I
M
E
X
-
I
M
E
X is plotted in Figures 1(a), 1(b), and 1(c), respectively. For IME
X
G
=
IMEX Euler, the stability functions of the parareal algorithms shown in Figures 1(a) and 1(b) are (24)
Figure
1(a):
R
IM-EX
(
x
,
y
)
=
1
/
(
x
+
i
y
-
1
)
-
1
+
(
x
+
i
y
)
/
M
M
1
-
1
/
(
x
+
i
y
-
1
)
,
R
IM-IMEX
(
x
,
y
)
=
1
/
(
x
+
i
y
-
1
)
-
(
M
+
i
y
)
/
M
-
x
M
1
-
1
/
(
x
+
i
y
-
1
)
,
R
IMEX-EX
(
x
,
y
)
=
(
1
+
i
y
)
/
(
1
-
x
)
-
1
+
(
x
+
i
y
)
/
M
M
1
-
(
1
+
i
y
)
/
(
1
-
x
)
,
R
IMEX-IMEX
(
x
,
y
)
=
(
1
+
i
y
)
/
(
1
-
x
)
-
(
M
+
i
y
)
/
M
-
x
M
1
-
(
1
+
i
y
)
/
(
1
-
x
)
,
Figure
1(b):
R
IM-EX
(
x
,
y
)
=
1
x
+
i
y
-
1
-
1
+
x
+
i
y
M
+
(
x
+
i
y
)
2
2
M
2
M
×
1
-
1
x
+
i
y
-
1
-
1
,
R
IM-IMEX
(
x
,
y
)
=
1
x
+
i
y
-
1
-
1
+
2
+
i
y
/
M
x
+
i
y
2
M
-
x
M
×
1
-
1
x
+
i
y
-
1
-
1
,
R
IMEX-EX
(
x
,
y
)
=
1
+
i
y
1
-
x
-
1
+
x
+
i
y
M
+
x
+
i
y
2
2
M
2
M
×
1
-
1
+
i
y
1
-
x
-
1
,
R
IMEX-IMEX
(
x
,
y
)
=
1
+
i
y
1
-
x
-
1
+
2
+
i
y
/
M
x
+
i
y
2
M
-
x
M
×
1
-
1
+
i
y
1
-
x
-
1
,where M
=
50. For the case IME
X
F
=
IMEX
(
4,4
,
3
), the stability functions of the four parareal algorithms are very complex and the presentations are omitted.Figure 1
(a)(
IME
X
G
,
IME
X
F
) = (IMEX Euler, IMEX Euler). (b) (
IME
X
G
,
IME
X
F
) = (IMEX Euler, IMEX Trapezoidal). (c) (
IME
X
G
,
IME
X
F
) = (IMEX Euler, IMEX
(
4,4
,
3
)).
(a)
(b)
(c)From Figure1, we see that the stability region of the P
I
M
E
X
-
I
M
E
X algorithm seems only a little smaller than that of the P
I
M
-
I
M
E
X algorithm, and both are significantly larger than that of the P
I
M
-
E
X algorithm. We note that the computation cost of P
I
M
-
I
M
E
X is obviously more expensive than that of P
I
M
E
X
-
I
M
E
X, since in some cases solving nonlinear equations is not involved in the P
I
M
E
X
-
I
M
E
X algorithm. Therefore, we get the conclusion that, for IMEX
G = IMEX Euler, it is better to use the P
I
M
E
X
-
I
M
E
X algorithm instead of P
I
M
-
I
M
E
X (in the sense of computational cost) and P
I
M
-
E
X (in the sense of stability). Moreover, it is clear that the computational cost of the P
I
M
E
X
-
E
X algorithm is the least one among the four algorithms, while from Figure 1 we see that the stability region of this algorithm is the smallest one, and particularly it is smaller than that of the P
I
M
-
E
X algorithm. However, from the point of computational cost of view, this algorithm still stands comparison with P
I
M
-
E
X.We next consider the caseIMEX
G = IMEX Trapezoidal scheme (see the left scheme of (9)). For IMEX
F = IMEX Euler and IMEX
F = IMEX Trapezoidal and IMEX
F
=
IMEX
(
4,4
,
3
), the stability regions of the four algorithms P
I
M
-
E
X, P
I
M
-
I
M
E
X, P
I
M
E
X
-
E
X, and P
I
M
E
X
-
I
M
E
X are plotted in Figures 2(a), 2(b), and 2(c), respectively. For IME
X
G = IMEX Trapezoidal, the stability functions of the parareal algorithms shown in Figures 2(a) and 2(b) are (25)
Figure
2(a):
R
IM-EX
(
x
,
y
)
=
2
+
x
+
i
y
2
-
x
-
i
y
-
1
+
x
+
i
y
M
M
×
1
-
2
+
x
+
i
y
2
-
x
-
i
y
-
1
,
R
IM-IMEX
(
x
,
y
)
=
2
+
x
+
i
y
2
-
x
-
i
y
-
M
+
i
y
M
-
x
M
×
1
-
2
+
x
+
i
y
2
-
x
-
i
y
-
1
,
R
IMEX-EX
(
x
,
y
)
=
1
+
2
+
i
y
x
+
i
y
2
-
x
-
1
+
x
+
i
y
M
M
×
1
-
1
+
2
+
i
y
x
+
i
y
2
-
x
-
1
,
R
IMEX-IMEX
(
x
,
y
)
=
1
+
2
+
i
y
x
+
i
y
2
-
x
-
M
+
i
y
M
-
x
M
×
1
-
1
+
2
+
i
y
x
+
i
y
2
-
x
-
1
,
Figure
2(b):
R
IM-EX
x
,
y
=
2
+
x
+
i
y
2
-
x
-
i
y
-
1
+
x
+
i
y
M
+
(
x
+
i
y
)
2
2
M
2
M
×
1
-
2
+
x
+
i
y
2
-
x
-
i
y
-
1
,
R
IM-IMEX
(
x
,
y
)
=
2
+
x
+
i
y
2
-
x
-
i
y
-
1
+
2
+
i
(
y
/
M
)
x
+
i
y
2
M
-
x
M
×
1
-
2
+
x
+
i
y
2
-
x
-
i
y
-
1
,
R
IMEX-EX
(
x
,
y
)
=
1
+
2
+
i
y
x
+
i
y
2
-
x
-
1
+
x
+
i
y
M
+
(
x
+
i
y
)
2
2
M
2
M
×
1
-
1
+
2
+
i
y
x
+
i
y
2
-
x
-
1
,
R
IMEX-IMEX
(
x
,
y
)
=
1
+
2
+
i
y
x
+
i
y
2
-
x
1
+
2
+
i
y
/
M
x
+
i
y
2
M
-
x
M
m
m
m
n
m
-
1
+
2
+
i
y
/
M
x
+
i
y
2
M
-
x
M
×
1
-
1
+
2
+
i
y
x
+
i
y
2
-
x
-
1
,where M
=
50. Again, for the case IME
X
F
=
IMEX
(
4,4
,
3
) we omit the presentations of the four stability functions because of high complexity.Figure 2
(a)(
IME
X
G
,
IME
X
F
) = (IMEX Trapezoidal, IMEX Euler). (b) (
IME
X
G
,
IME
X
F
) = (IMEX Trapezoidal, IMEX Trapezoidal), (
IME
X
G
,
IME
X
F
) = (IMEX Trapezoidal, IMEX
(
4,4
,
3
)).
(a)
(b)
(c)ForIME
X
G = IMEX Trapezoidal, from Figure 2 we see that the results are different from the previous case IME
X
G = IMEX Euler and are more interesting: (a) the stability region of each algorithm shrinks to a smaller region; (b) for each l algorithm, the stability regions for different IME
X
F are very similar; (c) it seems that there is no difference between P
I
M
-
E
X (P
I
M
E
X
-
E
X) and P
I
M
-
I
M
E
X (P
I
M
E
X
-
I
M
E
X).Based on the results shown in Figures1 and 4, we get the following two conclusions:(1)
it is more advisable to use IMEX Euler method as the underlying numerical method for theG
Δ
T propagator instead of using the IMEX Trapezoidal method;
(2)
ifIME
X
G = IMEX Trapezoidal, it seems that the P
I
M
-
E
X algorithm is the best one, since it outperforms P
I
M
E
X
-
E
X and P
I
M
E
X
-
I
M
E
X in the sense of stability and outperforms P
I
M
-
I
M
E
X in the sense of computational cost.


## 4. Numerical Results

In this section, we show some numerical results to illustrate the good performance of the parareal algorithm implemented with the IMEX RK method. We test the well-known Gray-Scott model arising from chemical reaction (see [38, 39]): (26)
u
t
=
ϵ
1
u
x
x
-
u
v
2
+
η
1
-
u
,
v
t
=
ϵ
2
v
x
x
+
u
v
2
-
(
η
+
θ
)
v
,where t
∈
[
0,20
], x
∈
[
0,1
] and ϵ
1
=
1
0
-
4, ϵ
2
=
1
0
-
6, η
=
0.024, and θ
=
0.06. The initial-boundary condition for (26) is chosen as (27)
u
(
0
,
x
)
=
1
-
1
2
sin
3
π
x
100
,
v
0
,
x
=
1
4
sin
3
π
x
100
,
u
(
t
,
1
)
=
u
(
t
,
0
)
=
1
,
v
(
t
,
1
)
=
v
(
t
,
0
)
=
0
.With these conditions, we expect three pulses in the solutions of (26). Figure 3 shows a typical behaviour of the solutions.Figure 3
Behavior of the solutions of the Gray-Scott model (26)-(27).
(a)
(b)Figure 4
Measured error diminishing of the solutionsu and v when the P
I
M
E
X
-
E
X algorithm is used.The PDE system (26) is first discretized spatially as u
x
x
≈
(
u
j
+
1
-
2
u
j
+
u
j
-
1
)
/
Δ
x
2 and v
x
x
≈
(
v
j
+
1
-
2
v
j
+
v
j
-
1
)
/
Δ
x
2 with Δ
x
=
0.01, j
=
1,2
,
…
,
100, and then a nonlinear ordinary differential system that consists of 200 ODEs is solved by the parareal algorithm. We use the IMEX Euler method to define both the G
Δ
T propagator and F
Δ
t propagator. We consider here M
=
50, Δ
T
=
1
/
10, and Δ
t
=
Δ
T
/
M. In Figure 4, we show the first 6 iterations of the parareal P
I
M
E
X
-
E
X algorithm, where we see that the error diminishing for the solutions u and v is very rapid.We also test the convergence speed of the pararealP
I
M
E
X
-
E
X algorithm and the parareal algorithm which is denoted by P
E
X
-
I
M
E
X with propagators G
Δ
T and F
Δ
t being defined by the explicit Euler and the IMEX Euler methods, respectively. The convergence curves corresponding to these three parareal algorithms are shown in Figure 5(a). We see that the P
E
X
-
I
M
E
X algorithm is not convergent and both the other two algorithms converge rapidly. This obversion can be explained by the stability region shown in Figure 5(b), where we see that the stability region of P
E
X
-
I
M
E
X is significantly smaller than that of the other two algorithms. Moreover, as has been shown in Figure 1, we see that the stability region of the P
I
M
E
X
-
E
X algorithm is nearly contained in the one of P
I
M
E
X
-
I
M
E
X. This means that, for the Gray-Scott model (26)-(27), if the P
I
M
E
X
-
E
X algorithm converges, so does the P
I
M
E
X
-
I
M
E
X algorithm. Moreover, it is interesting that the P
I
M
E
X
-
E
X algorithm converges a little sharper than P
I
M
E
X
-
I
M
E
X; see Figure 5(a).Figure 5
Measured error diminishing (a) and stability region (b) of the algorithmsP
E
X
-
I
M
E
X, P
I
M
E
X
-
E
X, and P
I
M
E
X
-
I
M
E
X.
(a)
(b)



---

*Source: 395340-2015-03-05.xml*