# Robust Convergence of Parareal Algorithms with Arbitrarily High-order Fine Propagators

The aim of this paper is to analyze the robust convergence of a class of parareal algorithms for solving parabolic problems. The coarse propagator is fixed to the backward Euler method and the fine propagator is a high-order single step integrator. Under some conditions on the fine propagator, we show that there exists some critical J_* such that the parareal solver converges linearly with a convergence rate near 0.3, provided that the ratio between the coarse time step and fine time step named J satisfies J ≥ J_*. The convergence is robust even if the problem data is nonsmooth and incompatible with boundary conditions. The qualified methods include all absolutely stable single step methods, whose stability function satisfies |r(-∞)|<1, and hence the fine propagator could be arbitrarily high-order. Moreover, we examine some popular high-order single step methods, e.g., two-, three- and four-stage Lobatto IIIC methods, and verify that the corresponding parareal algorithms converge linearly with a factor 0.31 and the threshold for these cases is J_* = 2. Intensive numerical examples are presented to support and complete our theoretical predictions.

## Authors

• 13 publications
• 2 publications
• 35 publications
• ### High order linearly implicit methods for evolution equations: How to solve an ODE by inverting only linear systems

This paper introduces a new class of numerical methods for the time inte...
11/14/2019 ∙ by Guillaume Dujardin, et al. ∙ 0

• ### Arbitrarily High-order Linear Schemes for Gradient Flow Models

We present a paradigm for developing arbitrarily high order, linear, unc...
10/16/2019 ∙ by Yuezheng Gong, et al. ∙ 0

• ### Irksome: Automating Runge–Kutta time-stepping for finite element methods

While implicit Runge–Kutta methods possess high order accuracy and impor...
06/29/2020 ∙ by Patrick E. Farrell, et al. ∙ 0

• ### On "Optimal" h-Independent Convergence of Parareal and MGRIT using Runge-Kutta Time Integration

Parareal and multigrid-reduction-in-time (MGRIT) are two popular paralle...
06/16/2019 ∙ by Stephanie Friedhoff, et al. ∙ 0

• ### A Two-Level Fourth-Order Approach For Time-Fractional Convection-Diffusion-Reaction Equation With Variable Coefficients

This paper develops a two-level fourth-order scheme for solving time-fra...
04/06/2021 ∙ by Eric Ngondiep, et al. ∙ 0

• ### Convergence analysis of multi-level spectral deferred corrections

The spectral deferred correction (SDC) method is class of iterative solv...
02/18/2020 ∙ by Gitte Kremling, et al. ∙ 0

• ### Does Higher Order LSTM Have Better Accuracy in Chunking and Named Entity Recognition?

Current researches usually employ single order setting by default when d...
11/22/2017 ∙ by Yi Zhang, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The main focus of this paper is to study the convergence of a class of parareal solver for the parabolic problems. Specifically, we let and consider the initial value problem of seeking satisfying

 {u′(t)+Au(t)=f(t),0

where is a positive definite, selfadjoint, linear operator with a compact inverse, defined in Hilbert space with domain dense in . Here is a given initial condition and is a given forcing term. Throughout the paper, denotes the norm of the space .

Parallel-in-time (PinT) methods, dating back to the work of Nievergelt in 1964 [23], have attracted a lot of interest in the last several decades. The parareal method, introduced in 2001 [20], is perhaps one of the most popular PinT algorithms. This method is relatively simple to implement, and can be employed for any single step integrators. In recent years, the parareal algorithm and some relevant algorithms, have been applied in many fields, such as turbulent plasma [25, 26], structural (fluid) dynamics [9, 12], molecular dynamics [4], optimal control [21, 22] Volterra integral equations and fractional models [19, 34], etc. We refer the interested reader to survey papers [13, 24] and references therein.

The parareal algorithm is defined by using two time propagators, and , associated with the large step size and the small step size respectively, where we assume that the ratio is an integer greater than . The fine time propagator is operated with small step size in each coarse sub-interval parallelly, after which the coarse time propagator is operated with large step size sequentially for corrections. In general, the coarse propagator is assumed to be much cheaper than the fine propagator . Therefore, throughout this paper, we fix to the backward-Euler method and study the choices of . Then a natural question arises related to convergence of the parareal algorithm. For parabolic type problems, in the pioneer work [5], Bal proved a fast convergence of the parareal method with a strongly stable coarse propagator and the exact fine propagator, provided some regularity assumptions on the problem data. The analysis works for both linear and nonlinear problems. This convergence behavior is clearly observed in numerical experiments, see e.g. Figure 5.1. However, without those regularity assumptions, the convergence observed from the empirical experiments will be much slower than expected, cf. Figure 5.2.

This interesting phenomenon motivates the current work, where we aim to study the convergence of parareal algorithm which is expected to be robust in the case of nonsmooth / incompatible problem data, that is related to various applications, e.g., optimal control, inverse problems, and stochastic models. There have existed some case studies. In [22], Mathew, Sarkis and Schaerer considered the backward Euler method as the fine propagator and proved the robust convergence of the parareal algorithm with a convergence factor (for all ); see also [14, 29] for some related discussion. In [32], Wu showed that the convergence factors for the second-order diagonal implicit Runge–Kutta method and a single step TR/BDF2 method (i.e., the ode23tb solver for ODEs in MATLAB) are (with ) and (with ), respectively. These error bounds might be slightly improved by increasing . See also [33] for the analysis for a third-order diagonal implicit Runge–Kutta method with a convergence factor (). For fourth-order Gauss–Runge–Kutta integrator, in [33]

Wu and Zhou showed that the threshold depends on both the largest eigenvalue of operator

and the step size . Note that the eigenvalues of may approach infinity, e.g. with homogeneous boundary conditions. Therefore, this kind of integrators might not be suitable for the parareal algorithm.

Then a natural question arises: in what case there exists a threshold (independent of step sizes , , terminal time , problem data and , as well as the distribution of spectrum of the elliptic operator ), such that for any , the parareal algorithm for solving the parabolic equation (1) converges robustly? Our study provides a positive answer to this question: if the fine propagator is strongly stable, in sense that the stability function satisfies , then there must exist such a positive threshold so that for all the parareal algorithm converges linearly with convergence factor close to . The convergence is robust even if the initial data is nonsmooth or incompatible with boundary conditions. Noting that all L-stable Runge–Kutta schemes satisfy that condition, so the fine propagator can be arbitrarily high-order. As examples, we analyzed three popular L-stable schemes, i.e., two-, three-, four-stage Lobatto IIIC schemes. We show that for all these cases the parareal algorithm converges linearly with factor less than and . Our theoretical results are fully supported by numerical experiments.

The rest of the paper is organized as follows. In Section 2, we introduce singe step integrators and parareal algorithms for solving the parabolic problem. Then we show the convergence of the algorithm in Section 3 by using the spectrum decomposition. Moreover, in Section 4

, we present case studies on three popular L-stable Runge–Kutta schemes, and show a sharper estimate for the threshold

. Finally, in Section 5, we present some numerical results to illustrate and complement the theoretical analysis.

## 2 Single step methods and parareal algorithm

In this section, we present the basic setting of the single step time stepping methods for solving the parabolic equation (1) and the parareal algorithm. See more detailed discussion in the monograph [30, Chapter 7-9] and the comprehensive survey paper [13].

### 2.1 Single step integrators for solving parabolic equations

To begin with, we consider the time discretization for the parabolic equation (1). We split the interval into subintervals with the uniform mesh size , and set , . Then a framework of a single step scheme approximating reads:

 un+1=r(−ΔtA)un+Δtm∑i=1pi(−ΔtA)f(tn+ciΔt),for all 0≤n≤N−1, (2)

Here, and are rational functions and are distinct real numbers in . Throughout the paper, we assume that the scheme (2) satisfies the following assumptions.

• and , for all , uniformly in and . Besides, the numerator of is of lower degree than its denominator.

• The time stepping scheme (2) is accurate of order in sense that

 r(−λ)=e−λ+O(λq+1),as λ→0.

and for

 m∑i=1cjipi(−λ)−j!(−λ)j+1(e−λ−j∑ℓ=0(−λ)ℓℓ!)=O(λq−j),as λ→0.
• The rational function is strongly stable in sense that .

###### Remark 1

Condition (P3) is essential for the convergence of parareal iteration. If , e.g., Crank-Nicolson method and implicit Runge-Kutta methods of Gauss type, the parareal method converges only if the eigenvalues of is bounded from above (which is not true for parabolic equations) and the ratio between the coarse step size and the fine step size is sufficiently large (depending on the upper bound of eigenvalues of ). Besides, this condition is also important in case that problem data is nonsmooth, e.g., . Time stepping schemes violating this condition may lose the optimal convergence rate in the nonsmooth data case [30, Chapter 8].

Practically, it is convenient to choose that share the same denominator of :

 r(λ)=a0(λ)g(λ),andpi(λ)=ai(λ)g(λ),for i=1,2,…,m,

where and are polynomials. Then the integrator (2) could be written as

 g(−ΔtA)un+1=a0(−ΔtA)un+Δtm∑i=1ai(−ΔtA)f(tn+ciΔt),for all 1≤n≤N.

See e.g. [30, pp. 131] for the construction of such rational functions satisfying (P1)-(P3).

Under those conditions, there holds the following error estimate for the time stepping scheme (2). The proof is given in [30, Theorems 7.2 and 8.3].

###### Lemma 1

Suppose that the Conditions (P1)-(P3) are fulfilled. Let be the solution to parabolic equation (1), and be the solution to the time stepping scheme (2). Then there holds

 ∥un−u(tn)∥≤c(Δt)q(t−qn∥u0∥+tnq−1∑ℓ=0sups≤tn∥Aq−ℓf(ℓ)(s)∥+∫tn0∥f(q)(s)∥ds),

provided that , with and for all .

###### Remark 2

Lemma 1 indicates that, under Conditions (P1)-(P3), the solution of the time stepping scheme (2) converges to the exact solution with order provided that the source term and initial condition satisfy certain compatibility conditions. For example, if we consider the parabolic equation where with homogeneous Dirichlet boundary condition, it requires on the boundary for . In order to avoid the restrictive compatibility conditions, we shall assume that the time discretization scheme (2) is strictly accurate of order in sense that

 m∑i=1cjipi(−λ)−j!(−λ)j+1(r(−λ)−j∑ℓ=0(−λ)ℓℓ!)=0,% for all 0≤j≤q−1.

It is well-known that a single step method with a given could be accurate of order (Gauss–Legendre method) [11, Section 2.2], but at most strictly accurate of order [6, Lemma 5].

###### Remark 3

The error estimate in Lemma 1 could be slightly improved if the time integrator is L-stable, i.e. ; see e.g., [30, Theorem 7.2].

### 2.2 Parareal algorithm

Next, we state the parareal solver for the single step scheme (2). Let , with a positive integer , be the coarse step size. Without loss of generality, we assume that is an integer, and let Then, two numerical propagators and are assigned to the coarse and fine time grids, where is usually a low-order and inexpensive numerical method (such as backward Euler scheme), and is given by the single step integrator (2). Specifically, for and , letting denote the identity operator, we define the coarse and finer propagator as

 G(Tn,ΔT,v,f)=(I+ΔTA)−1(v+ΔTf(Tn)).

and

 F(tn,Δt,v,f)=r(−ΔtA)v+Δtm∑i=1pi(−ΔtA)f(tn+ciΔt).

respectively. Then, the parareal solver is described in Algorithm 1.

The aim of this paper is to show that the iterative solution , generated by the parareal algorithm, linearly converges to the exact time stepping solution of the single step integrator (2) with fine time step , i.e.,

 max1≤n≤N∥Unk−unJ∥≤cγk, (3)

with some convergence factor strictly smaller than . We shall prove that there exists a positive threshold , independent of , and the upper bound of spectrum of , such that if , then (3) is true with close to , under conditions (P1)-(P3).

## 3 Convergence analysis

Next, we briefly test the convergence factor of parareal iteration. Taking comparision with the exact time stepping solution in (2), we arrive at

 Un+1k+1−u(n+1)J =(I+ΔTA)−1[(Unk+1−unJ)−(Unk−unJ)] +F(Tn+(J−1)Δt,Δt,˜Un,J−1,f)−F(Tn+(J−1)Δt,Δt,unJ−1,f)) =(I+ΔTA)−1[(Unk+1−unJ)−(Unk−unJ)]+r(−ΔtA)(˜Un,J−1−u(n+1)J−1) =⋯ =(I+ΔTA)−1[(Unk+1−unJ)−(Unk−unJ)]+r(−ΔtA)J(˜Un,0−unJ) =(I+ΔTA)−1[(Unk+1−unJ)−(Unk−unJ)]+r(−ΔtA)J(Unk−unJ),

For the sake of simplicity, we define and rewrite the above equation as

 En+1k+1=(I+ΔTA)−1(Enk+1−Enk)+r(−ΔtA)JEnk

Recall that the operator is a positive definite, selfadjoint, linear operator with a compact inverse, defined in Hilbert space . Then by the spectral theory, has positive eigenvalues , where and

, and the corresponding eigenfunctions

form an orthonormal basis of the Hilbert space . Then, letting , by means of spectrum decomposition, we derive

 en+1k+1,j=enk+1,j−enk,j1+ΔTλj+r(−Δtλj)Jenk,j.

By letting , we have

 en+1k+1,j≤(1+dj)−1enk+1,j+(r(−dj/J)J−(1+dj)−1)enk,j.

We apply the recursion and use the fact that , and hence obtain

 en+1k+1,j= (r(−dj/J)J−(1+dj)−1)enk,j+(1+dj)−1enk+1,j = (r(−dj/J)J−(1+dj)−1)(enk,j+(1+dj)−1en−1k,j)+(1+dj)−2en−1k+1,j = … = (r(−dj/J)J−(1+dj)−1)(enk,j+(1+dj)−1en−1k,j+⋯+(1+dj)−(n−1)e1k,j).

Now taking the absolute value on the both sides yields

 |en+1k+1,j|≤ |r(−dj/J)J−(1+dj)−1|⋅(1+(1+dj)−1+⋯+(1+dj)−(n−1))max1≤n≤N|enk,j| (4) ≤ |r(−dj/J)J−(1+dj)−1|1−(1+dj)−1max1≤n≤N|enk,j| = ∣∣∣(1+dj)r(−dj/J)J−1dj∣∣∣max1≤n≤N|enk,j| ≤ sups∈(0,∞)∣∣∣(1+s)r(−s/J)J−1s∣∣∣max1≤n≤N|enk,j|.

If the the leading factor is strictly smaller than one, i.e.,

 sups∈(0,∞)∣∣∣(1+s)r(−s/J)J−1s∣∣∣≤γ<1, (5)

then converges to zero linearly with a factor (smaller than) , and hence the parareal iteration converges linearly to the time stepping solution (2) in sense of (3).

Our convergence analysis in this and next sections heavily depend on the constant defined by

 κα:=sups∈(0,∞)∣∣∣(1+s)e−αs−1s∣∣∣,for  α∈[0,2]. (6)

To begin with, we establish an simple upper bound for the constant .

###### Lemma 2

Let , and be the constant defined in (6). Then there holds

 κα≤{eα−2,α∈[1,2];max(eα−2,1−α),α∈[0,1).

###### Proof

First of all, we show the claim that

 (1+s)e−αs−1s≥−eα−2. (7)

To this end, we define the auxiliary function

 g(s)=(1+s)e−αs+eα−2s.

Then a simple computation yields

 g′(s)=(1−α−αs)e−αs+eα−2andg′′(s)=(α2+α2s−2α)e−αs.

It is easy to observe that admits a single root at , and

 g′(x)≥g′((2−α)/α)=0.

Therefore is increasing in . As a result, , and hence

 (1+s)e−αs−1≥−eα−2s∀  s≥0,

which implies (7). Moreover, for , we observe that

 1+s≤es≤eαs∀s≥0,

which immediately leads to for all . This completes the proof for the desired results in case that .

Now we turn to the case that . Let and define

 g(x)=(1+s)e−αs−κ∗αs.

Then the simple computation yields

 g′(s)=(1−α−αs)e−αs−κ∗α% andg′′(s)=(α2+α2s−2α)e−αs.

Noting that , and . These imply and hence is decreasing function in . Therefore , which further implies

 (1+s)e−αs−κ∗αs≤1.

This leads to the desired assertion for the case that .

Lemma 2 only provides a rough upper bound for . In fact, for a fixed we can further improve the upper bound via a more careful computation. In Figure 3.1, we numerically compute the constant for and plot those values.

The next lemma provides a sharper estimate for .

###### Lemma 3

Let be the constant defined in (6). Then .

###### Proof

To show the sharp estimate, we note that for , there holds and . Therefore for all . This further implies for .

Meanwhile, we note that

 φ(s):=s2esψ′(s)=s2+s+1−es.

It is easy to verify that has a unique root, denoted by , in , so does . Therefore, . Using the Newton’s algorithm, we find and hence .

Now we state our main theorem which verifies the desired result (5) with .

###### Theorem 1

Let conditions (P1)-(P3) hold valid. Then there exists a threshold such that for all

 sups∈(0,∞)∣∣∣(1+s)r(−s/J)J−1s∣∣∣≤0.3,

###### Proof

First of all, we aim to show that

 limJ→∞sups≥0(1+s)s∣∣r(−s/J)J−e−s∣∣=0. (8)

For a given , for any , it is obvious that is bounded by a constant . Meanwhile, note that conditions (P1)-(P3) are fullfilled. Then by means of the nonsmooth data error estimate [30, Theorem 7.2], there holds

 ∣∣r(−s/J)J−e−s∣∣≤cJ−q,

where is independent of . Then we derive

 limJ→∞sups>δ(1+s)s∣∣r(−s/J)J−e−s∣∣=0.

For , conditions (P1) and (P2) imply

 (1+s)s∣∣r(−s/J)J−e−s∣∣ =(1+s)s|r(−s/J)−1−e−s/J|∣∣ ∣∣J−1∑i=0r(−s/J)−ie−(J−1−i)s/J∣∣ ∣∣ ≤C(α,δ)J−q.

Therefore we arrive at

 limJ→∞sup0

which completes the proof of the (8). This together with Lemma 3 implies that

which completes the proof of the lemma.

Next, using Theorem 1, we are able to show the linear convergence of the parareal iteration 1.

###### Theorem 2

Let conditions (P1)-(P3) be fullfilled and the data regularity in Lemma 1 hold valid. Let be the solution to the time stepping scheme (2), and be the solution obtained from the parareal algorithm 1. Then there exists a threshold such that for all , we have

 max1≤n≤Nc∥Unk−unJ∥≤cγkwithγ=0.3.

###### Proof

In Algorithm 1, the initial guess is obtained by the coarse propagator, i.e., the backward Euler scheme. Then Lemma 1 (with ) implies the estimate

 ∥Un0−unJ∥ ≤∥Un0−u(Tn)∥+∥u(Tn)−unJ∥ (9) ≤c((ΔT)T−1n+(Δt)t−1nJ)≤cn−1.

Let and . The the relation (4) and Theorem 1 imply

 sup1≤n≤Nc∥Enk∥2 ≤∞∑j=1sup1≤n≤Nc|enk,j|2≤γ2∞∑j=1sup1≤n≤Nc|enk−1,j|2 ≤⋯≤γ2k∞∑j=1sup1≤n≤Nc|en0,j|2

with . This together with the estimate leads to

 sup1≤n≤Nc∥Enk∥2 ≤cγ2k∞∑j=1Nc∑n=1|en0,j|2≤cγ2kNc∑n=1∥En0∥2≤cγ2kNc∑n=1n−2≤cγ2k,

where in the second last inequality we apply the estimate (9). This completes the proof of the theorem.

###### Remark 4

Theorem 2 provides an useful upper bound of the convergence factor for all single step integrators (satisfying (P1)-(P3)), which might not be sharp for specific one. For example, in [22, Lemma 4.3], Mathew, Sarkis and Schaerer considered the backward Euler method and proved that the convergence factor of the Parareal algorithm is around (with ). In [32], Wu showed that convergence factors are (with ) and (with ) for the second-order diagonal implicit Runge–Kutta method and a single step TR/BDF2 method (i.e., the ode23tb solver for ODEs in MATLAB), respectively. These error bounds might be slightly improved by increasing . See also [33] for the analysis for a third-order diagonal implicit Runge–Kutta method with a convergence factor and .

###### Remark 5

Theorem 2 only provides the existence of the threshold without any upper bound estimate. It is obvious that a huge may destroy the parallelism of the algorithm. Then a question arise naturally: is it possible to find for a given scheme satisfying conditions (P1)-(P3)? This is the focus of Section 4.

## 4 Case studies for several high-order single step integrators

In this section, we shall study some popular single step methods. As we mentioned in Remark 5, Theorem 2 did not provide an sharp estimate for the threshold . In fact, there is no universal estimate for all single step methods. Fortunately, for any given single step integrator satisfying conditions (P1)-(P3) and fixed convergence rate , we have a regular routine to find a sharper estimate for .

We consider three time-stepping methods, namely the the two-, three-, four-stage Lobatto IIIC methods, which are respectively second-, fourth- and sixth-order accurate, to the initial and boundary value problem (1). For the reader’s convenience, we present the Butcher tableaus of the two-, three-, four-stage Lobatto IIIC methods, respectively,

 (10)

,

 (11)

and

 (12)

Let us also briefly recall some well-known facts about Lobatto IIIC; for details we refer to [15]. These methods can be viewed as discontinuous collocation methods. The order of the -stage Lobatto IIIC methods is . In particular, the methods are algebraically stable and L-stable, that makes them suitable for stiff problems. The stability functions ,

 r(z):=1+zb⊤(I−zOι)−11with1:=(1,…,1)⊤∈Rq,

is given by the -Padé approximation to and vanishes at infinity, i.e., Note that the computational cost of implicit Runge-Kutta methods increases fast with the stage number, and we refer to [7, 17, 16] and the reference therein for some efficient implementations.

The following argument highly depends on the upper bound for the constant defined in (6). From Figure 3.1, we observe that Lemma 2 gives an sharp estimate for for , while the estimate for could be further improved. The next lemma provides an estimate for , which is useful in the analysis of convergence rate.

###### Lemma 4

Let be the constant defined in (6). Then .

###### Proof

With and , we observe that and . Meanwhile, since for , we derive that . Now we intend to show that admits a unique root in , denoted as . Then . Noting that

 ψ′(s)=1−(1+βs+βs2)e−βss2,

it suffices to show that has a unique root in . It is straightforward to see that the function

 g′(s)=(−2+β+βs)βse−βs

admits a unique root in . Then by the fact that and Rolle’s theorem, we conclude that has at most one root in . Meanwhile, we observe and . Therefore, there exists a unique root of in , named as , which lies in . Using the Newton’s algorithm, we find and hence .

###### Proposition 1

Let be the solution to the time stepping scheme (2) using the two-stage Lobatto IIIC method (10), and be the solution obtained from the parareal algorithm 1. Then for all , there holds

 max1≤n≤Nc∥Unk−un∥≤cγkwithγ=0.31.

###### Proof

It suffices to show that for any , there holds

 sups∈(0,∞)∣∣∣(1+s)r(−s/J)J−1s∣∣∣≤0.31,where  r(−s)=2s2+2s+2. (13)

To this end, we define and . Then Lemma 2 implies that , and meanwhile Lemma 4 indicates .

Next, we aim to show that for all with . First of all, using the fact that for all , we derive the first inequality . Then we turn to the second inequality , equivalent to in . Noting that , which admits a unique root at . Meanwhile, we observe that in , and in . Besides, since and , we conclude that in , and has a unique root in . Moreover, the facts and implies that there exists a constant , s.t. , and in , in . Then we note that and conclude that in , which implies in . As a result, we arrive at for all which implies

 sups∈(0,Js∗)∣∣∣