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.

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-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, [68]. Nowadays, this algorithm, as well as its variants [3, 912], 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, to formulate the parareal algorithm, the time domain 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 , 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 , 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 and 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 and . For example, Guibert and Tromeur-Dervout [22] considered the adaptive time-slices (i.e., the length of time-slices is not equal) to define and 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 and 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 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 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 [2333] 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 and/or 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.

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 into subintervals , , with equal length and . We call ’s the coarse time-grids. We then use some finer step size ( is an integer) to partition each relative large interval into finer intervals with and . 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 with , , for .

Step 1. Perform in subinterval the computation with initial value , .

Step 2. Perform sequential corrections with , .

Step 3. If for , satisfy some termination criteria, stop the iteration; else go to Step 1.

Here, the finer propagation operator and the coarse propagation operator 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 . This implies that the converged solution will achieve the accuracy of the propagator .

2.2. 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 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 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 . We will use the triplet to identify scheme (4), where is the number of stages of the implicit part, 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 ):where and . Both the implicit and explicit schemes of 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.

The asymptotic behavior of the IMEX RK methods applied to (6) is analyzed on the basis of the scalar model equation: 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 where we have used vector and matrix notations to denote the values shown in the Butcher tableau (4); that is, , , , and . Moreover, .

Let , . We then have where is called stability function of the IMEX RK method which is defined by where is an identity matrix of size .

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 matrix of size . Its elements are defined by The next necessity for our analysis is the estimation of , which is proved by Gander and Vandewalle [7].

Lemma 1 (see [7]). Let be an integer. Then the infinite norm of can be estimated by

Assume that the stability functions of the underlying numerical methods used to define the two propagators and are and , respectively. Then it follows by applying the parareal algorithm to the model equation (11) that where we recall .

Theorem 2. Let be given and for .  Let the underlying numerical method used to define the propagator running on the coarse time-grids be in its region of absolute stability; that is, . Then one has the following estimates:

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 for any . Set . Then relation (19) can be written in matrix form as where the matrix is defined by (15) with . This implies Therefore, from Lemma 1 and the assumption , we have Hence, the proof of this theorem is completed by letting in the second inequality of (22).

The factor 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.

DefineThen if , the parareal algorithm is convergent on unbounded time intervals. In what follows of this paper, the set is called stability region of the parareal algorithm.

Remark 3. If and (i.e., the IMEX RK methods reduce to the traditional RK methods), it follows from (14) that with . Therefore, if both the propagators and are defined by the traditional RK methods, the stability function can be rewritten as , which is the one obtained by Gander and Vandewalle in [7].

From Remark 3, we see that the underlying numerical method used to define propagator 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 and are defined by the traditional RK methods), it has been proposed by Farhat et al. [4, 5, 14] that the underlying RK method for should be implicit and the one for the propagator 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, and , we consider here and the following four special combinations which leads to four parareal algorithms::  is defined by the implicit part of and is defined by the explicit part ;:   is defined by the implicit component of and is defined by ;:   is defined by and is defined by the explicit component of ;:   is defined by and is defined by .For the convenience of plotting the stability region of the parareal algorithm, we consider the case , . The stability functions of the above four parareal algorithms are denoted by , , , and , respectively. The stability regions are denoted by , , , and , respectively. We remark that the 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 is chosen as the IMEX Euler method (see (8)). For = IMEX Euler and = IMEX Trapezoidal and , the stability region of the algorithms , , , and is plotted in Figures 1(a), 1(b), and 1(c), respectively. For Euler, the stability functions of the parareal algorithms shown in Figures 1(a) and 1(b) are where . For the case , 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 algorithm seems only a little smaller than that of the algorithm, and both are significantly larger than that of the algorithm. We note that the computation cost of is obviously more expensive than that of , since in some cases solving nonlinear equations is not involved in the algorithm. Therefore, we get the conclusion that, for = IMEX Euler, it is better to use the algorithm instead of (in the sense of computational cost) and (in the sense of stability). Moreover, it is clear that the computational cost of the 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 algorithm. However, from the point of computational cost of view, this algorithm still stands comparison with .

We next consider the case = IMEX Trapezoidal scheme (see the left scheme of (9)). For = IMEX Euler and = IMEX Trapezoidal and , the stability regions of the four algorithms , , , and are plotted in Figures 2(a), 2(b), and 2(c), respectively. For = IMEX Trapezoidal, the stability functions of the parareal algorithms shown in Figures 2(a) and 2(b) are where . Again, for the case we omit the presentations of the four stability functions because of high complexity.

For = IMEX Trapezoidal, from Figure 2 we see that the results are different from the previous case = 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 are very similar; (c) it seems that there is no difference between () and ().

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 propagator instead of using the IMEX Trapezoidal method;(2)if = IMEX Trapezoidal, it seems that the algorithm is the best one, since it outperforms and in the sense of stability and outperforms 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]): where , and , , , and . The initial-boundary condition for (26) is chosen as With these conditions, we expect three pulses in the solutions of (26). Figure 3 shows a typical behaviour of the solutions.

The PDE system (26) is first discretized spatially as and with , , 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 propagator and propagator. We consider here , , and . In Figure 4, we show the first 6 iterations of the parareal algorithm, where we see that the error diminishing for the solutions and is very rapid.

We also test the convergence speed of the parareal algorithm and the parareal algorithm which is denoted by with propagators and 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 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 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 algorithm is nearly contained in the one of . This means that, for the Gray-Scott model (26)-(27), if the algorithm converges, so does the algorithm. Moreover, it is interesting that the algorithm converges a little sharper than ; see Figure 5(a).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors are very grateful to the anonymous referees for the careful reading of a preliminary version of the paper and their valuable suggestions and comments, which really improve the quality of this paper. This work was supported by Science & Technology Bureau of Sichuan Province (2014JQ0035), Non-Destructive Inspection of Sichuan Institutes of Higher Education (2013QZY01), the project sponsored by OATF-UESTC, and the NSF of China (11301058, 11301362, 11371157, and 91130003).