Parareal Algorithms Implemented with IMEX Runge-Kutta Methods

Parareal algorithm is a very powerful parallel computation method and has received much interest frommany 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 RKmethods 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.


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-called resonance or beating 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][7][8].Nowadays, this algorithm, as well as its variants [3,[9][10][11][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,   () =  (,  ()) ,  ∈ [0, ] , to formulate the parareal algorithm, the time domain [0, ] 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 Δ , 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 Δ , 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 Δ and F Δ 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 define G Δ and F Δ .For example, Guibert and Tromeur-Dervout [22] considered the adaptive time-slices (i.e., the length of timeslices is not equal) to define G Δ and F Δ 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 Δ and F Δ 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 Δ 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 function  in (1) can be written into two parts where   is the nonstiff or mildly stiff part of  and   is the stiff part.A common example for such case arises from semidiscretization of reaction diffusion equations   =   + (, , ).In this case the function  in (1) can be written as (, ()) = () +   (, ()) with  being a tridiagonal matrix and stiff.In the case that  can be partitioned into stiff and nonstiff parts, it is more advisable to use the implicitexplicit (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][24][25][26][27][28][29][30][31][32][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 Δ and/or F Δ 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 Section 2, 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.

The Parareal Algorithm and the IMEX RK Methods
In this section, we revisit the parareal algorithm and the IMEX RK methods.

The Parareal
satisfy some termination criteria, stop the iteration; else go to Step 1.
Here, the finer propagation operator F Δ and the coarse propagation operator G Δ 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 where  is the iteration index.Note that for  → +∞ method (3) will generate a series of values   which satisfy  +1 = F Δ (  ,   , Δ).This implies that the converged solution   will achieve the accuracy of the propagator F Δ .

The IMEX RK Methods. Let us consider a pair of two
Runge-Kutta methods defined by the arrays with the same abscissae 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 ℎ > 0 be a step size and define the step point   = ℎ for integer .Consider the following stiff/nonstiff partitioned ODEs: and by applying the top part of (4) to stiff component   and the bottom part to the nonstiff component, we obtain the following scheme: where ∑  =1 (⋅) = 0 with  < .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  = 1, 2, 3. We will use the triplet IMEX( im ,  ex , ) to identify scheme (4), where  im is the number of stages of the implicit part,  ex is the number of stages of the explicit part, and  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: Second Order IMEX RK Methods.The second order IMEX RK methods considered here are the IMEX trapezoidal scheme [34]: The 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) where  = 0.4358665215 and  = 0.552929147.Both the implicit and explicit schemes of IMEX(4, 4, 3) are third order RK methods, and the implicit scheme is -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.

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 matrix M := M() of size .Its elements are defined by The next necessity for our analysis is the estimation of ‖M  ‖ ∞ , which is proved by Gander and Vandewalle [7].
Lemma 1 (see [7]).Let  > 0 be an integer.Then the infinite norm of M  can be estimated by Assume that the stability functions of the underlying numerical methods used to define the two propagators G Δ and F Δ are  F (, ) and  G (, ), respectively.Then it follows by applying the parareal algorithm to the model equation (11) that where we recall  = Δ/Δ.
Proof.We denote by    the error at iteration  of the parareal algorithm at coarse time point   ; that is,    = (  ) −    .Then with an induction argument on , it follows by applying the iterative process (3) to the model equation (11) that this error satisfies where we have used the fact that   0 = 0 for any  ≥ 0. Set   = (  1 ,   2 , . . .,    )  .Then relation ( 19) can be written in matrix form as where the matrix M is defined by (15) with  =  G (, ).This implies Hence, the proof of this theorem is completed by letting  → +∞ in the second inequality of (22).
The factor (, ) := |  F (/, /) −  G (, )|/(1 − | G (, )|) is called by Gander and Vandewalle [7] the linear 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 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.From Remark 3, we see that the underlying numerical method used to define propagator G Δ 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 Δ and F Δ 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 Δ should be implicit and the one for the propagator F Δ should be explicit.The reason is obviousthe 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, IMEX G and IMEX F , we consider here  = 50 and the following four special combinations which leads to four parareal algorithms: For the convenience of plotting the stability region of the parareal algorithm, we consider the case (, ) = (, ), ,  ∈ R. The stability functions of the above four parareal algorithms are denoted by R IM-EX (, ), R IM-IMEX (, ), R IMEX-EX (, ), and R IMEX-IMEX (, ), 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 IM-EX 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 when IMEX G is chosen as the IMEX Euler method (see (8)).For IMEX F = IMEX Euler and IMEX F = IMEX Trapezoidal and IMEX F = IMEX(4, 4, 3), the stability region of the algorithms P IM-EX , P IM-IMEX , P IMEX-EX , and P IMEX-IMEX is plotted in Figures 1(a), 1(b), and 1(c), respectively.For IMEX G = IMEX Euler, the stability functions of the parareal algorithms shown in Figures 1(a    where  = 50.For the case IMEX F = IMEX(4, 4, 3), the stability functions of the four parareal algorithms are very complex and the presentations are omitted.
From Figure 1, we see that the stability region of the P IMEX-IMEX algorithm seems only a little smaller than that of the P IM-IMEX algorithm, and both are significantly larger than that of the P IM-EX algorithm.We note that the computation cost of P IM-IMEX is obviously more expensive than that of P IMEX-IMEX , since in some cases solving nonlinear equations is not involved in the P IMEX-IMEX algorithm.Therefore, we get the conclusion that, for IMEX G = IMEX Euler, it is better to use the P IMEX-IMEX algorithm instead of P IM-IMEX (in the sense of computational cost) and P IM-EX (in the sense of stability).Moreover, it is clear that the computational cost of the P IMEX-EX 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 IM-EX algorithm.However, from the point of computational cost of view, this algorithm still stands comparison with P IM-EX .
We next consider the case IMEX G = IMEX Trapezoidal scheme (see the left scheme of ( 9 where  = 50.Again, for the case IMEX F = IMEX(4, 4, 3) we omit the presentations of the four stability functions because of high complexity.Based on the results shown in Figures 1 and 4, we get the following two conclusions: (1) it is more advisable to use IMEX Euler method as the underlying numerical method for the G Δ propagator instead of using the IMEX Trapezoidal method;
defined by the implicit part of IMEX G and F Δ is defined by the explicit part IMEX F ; P IM-IMEX : G Δ is defined by the implicit component of IMEX G and F Δ is defined by IMEX F ; P IMEX-EX : G Δ is defined by IMEX G and F Δ is defined by the explicit component of IMEX F ; P IMEX-IMEX : G Δ is defined by IMEX G and F Δ is defined by IMEX F .

Figure 4 :
Figure 4: Measured error diminishing of the solutions  and V when the P IMEX-EX algorithm is used.

Theorem 2 .
Let Δ be given and   = Δ for  = 0, 1, . ... Let the underlying numerical method used to define the propagator G Δ running on the coarse time-grids be in its region of absolute stability; that is, | G | < 1.Then one has the following estimates: sup >0       (  ) −  0       ,  ≥ 1.