On parareal algorithms for semilinear parabolic Stochastic PDEs

TL;DR

This paper investigates the performance of parareal algorithms applied to semilinear parabolic stochastic PDEs, focusing on the choice of coarse integrator and noise regularity, supported by theoretical analysis and numerical experiments.

Contribution

It introduces the use of exponential Euler schemes as coarse integrators in parareal algorithms for stochastic PDEs and analyzes their impact on convergence and efficiency.

Findings

Exponential Euler as coarse integrator improves convergence.

Noise regularity significantly affects algorithm performance.

Number of iterations influences accuracy and computational cost.

Abstract

Parareal algorithms are studied for semilinear parabolic stochastic partial differential equations. These algorithms proceed as two-level integrators, with fine and coarse schemes, and have been designed to achieve a `parallel in real time' implementation. In this work, the fine integrator is given by the exponential Euler scheme. Two choices for the coarse integrator are considered: the linear implicit Euler scheme, and the exponential Euler scheme. The influence on the performance of the parareal algorithm, of the choice of the coarse integrator, of the regularity of the noise, and of the number of parareal iterations, is investigated, with theoretical analysis results and with extensive numerical experiments.

Equations

$$\frac{\partial u(t,x)}{\partial t}=$$

$$u(0,x)=$$

$$u_{n+1}^{(k+1)}$$

$$u_0^{(k)}$$

$$\Vert (-A{)}^{\alpha -\frac{1}{2}}{Q}^{\frac{1}{2}}{\Vert }_{{\mathcal{L}}_2(H)}<\infty ,$$

$$du(t)=Au(t)dt+F(u(t))dt+d{W}^Q(t),$$

$$\{\begin{array}{l}D(A)={\mathbb{H}}^2(0,1)\cap {\mathbb{H}}_0^1(0,1),\\ Au=u^{\prime \prime },u\in D(A).\end{array}$$

$$e^{tA}u=\sum _{p\in \mathbb{N}}{e}^{-{\lambda }_{p}t}⟨u,e_p⟩e_p.$$

$$|u{|}_{\alpha }^2=\sum _{p\in \mathbb{N}}{\lambda }_p^{2\alpha }⟨u,e_p{⟩}^2\in [0,\infty ].$$

$$\Vert (-A{)}^{\alpha }{e}^{tA}{\Vert }_{\mathcal{L}(H)}\le C_{\alpha }\min (t,1{)}^{-\alpha },\Vert (-A{)}^{-\alpha }(e^{tA}-I){\Vert }_{\mathcal{L}(H)}\le C_{\alpha }\min (t,1{)}^{\alpha }.$$

$$|DF(u).h|_{-\alpha}\leq C_{F,\alpha,\kappa}\bigl{(}1+|u|_{\alpha+\kappa}\bigr{)}|h|_{-\alpha},$$

$$|DF(u).h|_{\alpha-\kappa}\leq C_{F,\alpha,\kappa}\bigl{(}1+|u|_{\alpha}\bigr{)}|h|_{\alpha}.$$

$$\big{|}\bigl{(}DF(u_{2})-DF(u_{1})\bigr{)}.h\big{|}_{-\alpha}\leq C_{F,\alpha,\kappa}\bigl{(}1+|u_{1}|_{\alpha+\kappa}+|u_{2}|_{\alpha+\kappa}\bigr{)}|u_{2}-u_{1}|_{-\alpha}|h|_{\alpha+\kappa}.$$

$$W(t)=\sum _{p\in \mathbb{N}}{\beta }_p(t){ϵ}_p.$$

$$W^Q(t)=\sum _{p\in \mathbb{N}}\sqrt{{\gamma }_p}{\beta }_p(t){ϵ}_p,$$

$$Q^{\frac{1}{2}}u=\sum _{p\in \mathbb{N}}\sqrt{{\gamma }_p}⟨x,{ϵ}_p⟩{ϵ}_p,Qu=\sum _{p\in \mathbb{N}}{\gamma }_p⟨x,{ϵ}_p⟩{ϵ}_p,\forall u\in H.$$

$$\overline{\alpha }=\sup \left\{\alpha \in (0,\infty ),\Vert (-A{)}^{\alpha -\frac{1}{2}}{Q}^{\frac{1}{2}}{\Vert }_{{\mathcal{L}}_2(H)}<\infty \right\},$$

$$\vvvert U\vvvert_{q}=\bigl{(}\mathbb{E}[|U|^{q}]\bigr{)}^{\frac{1}{q}},~{}\vvvert U\vvvert_{q,\alpha}=\bigl{(}\mathbb{E}[|U|_{\alpha}^{q}\bigr{)}^{\frac{1}{q}}.$$

$$u(t)=e^{tA}{u}_0+{\int }_0^t{e}^{(t-s)A}F(u(s))ds+{\int }_0^t{e}^{(t-s)A}d{W}^Q(s).$$

$$\underset{0\leq t\leq T}{\sup}~{}\vvvert u(t)\vvvert_{q}\leq C_{T,q}\bigl{(}1+|u_{0}|\bigr{)}.$$

$$\text{vvvert}u(t){\text{vvvert}}_{q,\alpha }$$

$$\text{vvvert}u(t)-u(s){\text{vvvert}}_q$$

$$t_n=n\Delta T,t_{n,j}=t_n+j\delta t=(nJ+j)\delta t.$$

$$\mathcal{G}(u,t_{n},t_{n+1})=\hat{S}_{\Delta T}u+\Delta T\hat{S}_{\Delta T}F(u)+\hat{S}_{\Delta T}\bigl{(}W^{Q}(t_{n+1})-W^{Q}(t_{n})\bigr{)},$$

$$\mathcal{F}_{\rm aux}(u,t_{n,j},t_{n,j+1})=e^{\delta tA}u+\delta te^{\delta tA}F(u)+e^{\delta tA}\bigl{(}W^{Q}(t_{n,j+1})-W^{Q}(t_{n,j})\bigr{)}.$$

$$\mathcal{F}(\cdot ,t_n,t_{n+1})={\mathcal{F}}_{\mathrm{aux}}(\cdot ,t_{n,J-1},t_{n,J})\circ \cdots \circ {\mathcal{F}}_{\mathrm{aux}}(\cdot ,t_{n,0},t_{n,1}).$$

$$v_{n,j+1}$$

$$v_{n,0}$$

$$u_{n+1}^{(0)}$$

$$u_0^{(0)}$$

$$u_{n+1}^{(k+1)}$$

$$u_0^{(k)}$$

$$u_{n+1}^{\mathrm{ref}}$$

Full text

The influence on the performance of the parareal algorithm, of the choice of the coarse integrator, of the regularity of the noise, and of the number of parareal iterations, is investigated, with theoretical analysis results and with extensive numerical experiments.

1. Introduction

In the last two decades, numerical methods for Stochastic Partial Differential Equations (SPDEs) have been extensively studied, see for instance the monographs [13], [15], [18] and references therein. The rate of convergence of the schemes used for temporal and spatial discretization is related to the regularity of the noise, which may be arbitrarily low, and in such situation effective numerical methods are difficult to construct. For instance, for a one-dimensional semilinear parabolic SPDE, driven by Gaussian space-time white noise, trajectories are only $\alpha$-Hölder continuous in time and $2\alpha$-Hölder continuous in space, for $\alpha<\frac{1}{4}$, and standard Euler and finite difference schemes thus have a (strong) rate of convergence equal to $\alpha$ and $2\alpha$ respectively.

In this article, we will only focus on the temporal discretization and consider a semilinear parabolic SPDE of the type

[TABLE]

that is, a one-dimensional semilinear heat equation with homogeneous Dirichlet boundary conditions, and with additive Gaussian noise, which is white in time and colored in space. A rigorous interpretation as a stochastic evolution equation (in the sense of [5]), driven by a Wiener process, with values in an infinite dimensional Hilbert space is given by (1), see also Section 2, where precise assumptions for the nonlinear operator $f$ and the covariance operator $Q$ are stated. It is well-known that the rate of convergence (studied only in the strong sense in this article) of the error depends on the properties of the covariance operator $Q$ of the Wiener process, and can be arbitrarily small. In particular, in the case of space-time white noise ($Q$ is the identity), the order of convergence of Euler schemes is essentially $\frac{1}{4}$. If the noise is of trace-class ($Q$ has finite trace), the order of convergence is essentially $\frac{1}{2}$.

The objective of this article is to study the applicability of so-called parareal algorithms to improve the performance of temporal discretization schemes. Implementing parallel algorithms to solve time-dependent evolution equations is not natural, and parallel-in-time integration methods have been extensively studied, based on multigird or multiple shooting techniques. We refer to the pionneering contributions [3, 4, 12, 14, 19, 22, 23], and the review [6]. The parareal algorithm has been introduced in [16], see also [21]. It is a parallel in real time technique, based on the use of two integrators with two different time-step sizes: a coarse integrator, denoted by $\mathcal{G}$, with coarse time-step size $\Delta T$, and a fine integrator, denoted by $\mathcal{F}$, using $J$ steps with time-step size $\delta t$, such that $\Delta T=J\delta t$. In practice, $\mathcal{G}$ may be less accurate but cheaper than $\mathcal{F}$. The parareal algorithm is an iterative method, using a predictor-corrector strategy, in which computations of the fine integrator at each iteration are performed in parallel: the recursion is given by

[TABLE]

where $k$ represents the index for parareal iterations. For $k=0$, the scheme is initialized using the coarse integrator. We refer to Equation (5) and to Section 3 below for a detailed presentation of the algorithm.

Since the pioneering work [16], where the parareal algorithm has been introduced for a class of ordinary differential equations, several extensions have been considered, see for instance [1, 7, 8, 9, 10, 20]. In particular, [1] deals with parabolic PDEs, and studies the stability and convergence properties, which may require regularity properties, depending on the choice of integrators. The application of the parareal algorithm for stochastic systems has been considered first in [2], and more recently in [11] for stochastic Schrödinger PDEs and in [24] for a class of stochastic differential equations. More precisely, in [11], parareal algorithms for stochastic Schrödinger equation with damping are studied with $\mathcal{F}$ being the exact solver and $\mathcal{G}$ being the exponential-$\theta$ scheme. The longterm convergence is obtained for the case $\theta>\frac{1}{2}$ or sufficient large $\alpha$, which ensures sufficient exponential decay of the coarse integrator.

Let us now describe the contributions of this article. The parareal algorithm is applied to the SPDE above. The exponential Euler scheme is chosen as the fine integrator. The main contribution of this article is to reveal that the parareal algorithm behaves differently, depending on the choice of the coarse integrator, when applied to semilinear parabolic SPDEs. Such results, based on both theoretical analysis and numerical experiments, have not been reported before, up to our knowledge. The error of the parareal algorithm, considered in this article, is the distance between the solution computed by the proposed algorithm and a reference solution generated by the fine integrator (which is not computationally expensive and is not computed in practice). The parareal algorithm is useful to reduce computational cost only if the order of convergence of this error (with respect to the coarse time-step size $\Delta T$) is strictly larger for some $k\geq 1$ than for $k=0$.

To state the main results of this article (see Theorems 4.1, 5.2 and 5.3 below for precise statements), let us assume that, for some $\alpha>0$, the covariance operator $Q$ satisfies a condition of the type

[TABLE]

where $\|\cdot\|_{\mathcal{L}_{2}(H)}$ denotes the Hilbert-Schmidt norm for operators from $H$ to $H$.

First, assume that the linear implicit Euler scheme is chosen as the coarse integrator, i.e. $e^{\Delta TA}$ is approximated by $(I-\Delta TA)^{-1}$. In addition, assume that $F=0$. Then the order of convergence of the parareal algorithm, with $k$ iterations, is essentially $\min(\alpha,k+1)$, and saturates at $\alpha$ when $k$ increases, see Theorem 4.1 and numerical experiments in Section 4.2. In particular, if $\alpha\leq 1$, which includes the space-time white noise case with essentially $\alpha=\frac{1}{4}$, the application of the parareal algorithm is useless, since this order of convergence does not depend on $k$. The way the error behaves in terms of $k$ thus depends a lot on the regularity of the noise.

Second, assume that the exponential Euler scheme is chosen as the coarse integrator. Then it is proved that the order of convergence of the parareal algorithm, with $k$ iterations, is at least essentially of size $(k+1)\alpha$, if $\alpha\in(0,\frac{1}{2}]$, see Theorem 5.2. Contrary to the first case, this order is linear in $k$, thus applying the parareal algorithm always reduces the computational cost, whatever the regularity of the noise. Numerical experiments, see Section 5.2, reveal that the order of convergence is sharp for $k=0$ (the parareal algorithm is not applied) and $k=1$ (one iteration is applied), but is larger when $k\geq 2$. In fact, the choice $k=1$ is optimal when considering the final computational cost. Theorem 5.3 proves that for $\alpha\in(0,\frac{1}{4}]$ and $k\geq 2$, the order of convergence is at least of size $2k\alpha$, which is indeed larger than $(k+1)\alpha$.

In conclusion, the parareal algorithm may offer an effective strategy to reduce computational cost for the simulation of trajectories of SPDEs. Several questions remain open: for instance, generalizations in higher dimension, algorithms for equations with multiplicative noise, or using other integrators, are left for future works.

This article is organized as follows. Precise assumptions on the operators $A$ and $F$, and on the covariance operator $Q$, are provided in Section 2. Section 3 is devoted to introducing the parareal algorithm, to presenting the possible choices of coarse integrators, and to defining the error. The study of the behavior of the parareal algorithm when the linear implicit Euler scheme is chosen as the coarse integrator is provided in Section 4: more precisely, see Theorem 4.1 for the theoretical error estimates and Section 4.2 for the numerical experiments. The study of the behavior of the parareal algorithm when the exponential Euler scheme is chosen as the coarse integrator is provided in Section 5: more precisely, see Theorems 5.2 and 5.3 for the theoretical error estimates and Section 5.2 for the numerical experiments.

2. Setting

Let $H=L^{2}(0,1)$, with norm and inner product denoted by $|\cdot|$ and $\langle\cdot,\cdot\rangle$ respectively. In this work, stochastic evolution equations, with additive noise, of the following type are considered (see [5], [17]: given an initial condition $u(0)=u_{0}\in H$,

[TABLE]

where the solution \bigl{(}u(t)\bigr{)}_{t\geq 0} is a stochastic process with values in $H$. The operator $A$ is a linear parabolic differential operator with homogeneous Dirichlet boundary conditions (Section 2.1), the operator $F$ is a globally Lipschitz, non-linear operator (Section 2.2), and \bigl{(}W^{Q}(t)\bigr{)}_{t\geq 0} is a $Q$-Wiener process. Appropriate assumptions to justify the global well-posedness of (1) are introduced below.

In the sequel, the initial condition $u_{0}$ is assumed to be deterministic, however the extension to a random initial condition $u_{0}$ (independent of the Wiener process, and satisfying appropriate moment conditions) is straightforward by a conditioning argument.

2.1. Linear operator $A$

The linear operator $A$ is defined as the unbounded linear operator on $H=L^{2}(0,1)$, such that

[TABLE]

In other words, $A$ is the Laplace operator with homogeneous Dirichlet boundary conditions. Recall that $A$ is an unbounded, self-adjoint, linear operator, and that $Ae_{p}=-\lambda_{p}$ for all $p\in\mathbb{N}$, where the eigenvalues are given by $\lambda_{p}=(\pi p)^{2}$, and the eigenfunctions e_{p}=\sqrt{2}\sin\bigl{(}p\pi\cdot\bigr{)} form a complete orthonormal system of $H$.

The linear operator $A$ generates an analytic and strongly continuous semigroup on $H$, denoted by \bigl{(}e^{tA}\bigr{)}_{t\geq 0}. Note that for all $u\in H$, one has

[TABLE]

For any $\alpha\in[0,1]$ and $u\in H$, let

[TABLE]

For $u\in D((-A)^{\alpha})=\{u\in H:|u|_{\alpha}<\infty\}$, set $(-A)^{\alpha}u=\sum_{p\in\mathbb{N}}\lambda_{p}^{\alpha}\langle u,e_{p}\rangle e_{p}\in H$, and note that $|u|_{\alpha}=|(-A)^{\alpha}u|$. In addition, for any $\alpha\in[0,1]$ and $u\in H$, let $(-A)^{-\alpha}u=\sum_{p\in\mathbb{N}}\lambda_{p}^{-\alpha}\langle u,e_{p}\rangle e_{p}\in H$, and $|u|_{-\alpha}^{2}=\sum_{p\in\mathbb{N}}\lambda_{p}^{-2\alpha}\langle u,e_{p}\rangle^{2}$.

Regularization properties of the semigroup \bigl{(}e^{tA}\bigr{)}_{t\geq 0} are stated in Proposition 2.1 below. The following notation is used. First, $\mathcal{L}(H)$ is the space of linear bounded operators from $H$ to $H$, with the operator norm denoted by $\|\cdot\|_{\mathcal{L}(H)}$. Second, $\mathcal{L}_{2}(H)$ is the space of Hilbert-Schmidt operators from $H$ to $H$, with the Hilbert-Schmidt norm denoted by $\|\cdot\|_{\mathcal{L}_{2}(H)}$.

Proposition 2.1**.**

For all $t\geq 0$, $\|e^{tA}\|_{\mathcal{L}(H)}\leq e^{-\lambda_{1}t}$. Moreover, for all $\alpha\in[0,1]$, there exists $C_{\alpha}\in(0,\infty)$ such that for all $t\in(0,\infty)$,

[TABLE]

2.2. Nonlinear operator $F$

The analysis of the rate of convergence for parareal algorithm below proceeds in a simplified, abstract, framework, whereas numerical experiments are performed in the more natural framework of Nemytskii operators. The abstract framework does not encompass this case. Indeed, the treatment of Nemytskii would require the introduction of further concepts (such as $\gamma$-Radonifying operators, in order to work in Banach spaces $L^{p}(0,1)$). Instead of increasing the technical level of the presentation, the choice made for this article is to study the main features of the parareal algorithms applied to SPDEs in a more pedagogical way owing to the simplified framework.

Let us first state the assumptions on the nonlinear operator $F$ which are employed for the theoretical analysis.

Assumption 1**.**

The nonlinear operator $F:H\to H$ is globally Lipschitz continuous, and is twice Fréchet differentiable, with bounded first and second order derivatives. Moreover, for any $\alpha\in(0,\frac{1}{2})$ and any arbitrarily small $\kappa\in(0,\frac{1}{2}-\alpha)$, there exists $C_{F,\alpha,\kappa}\in(0,\infty)$ such that, for all $u\in D((-A)^{\alpha+\kappa})$ and $h\in H$,

[TABLE]

and for all $u,h\in D((-A)^{\alpha})$,

[TABLE]

Finally, for all $u_{1},u_{2},h\in D((-A)^{\alpha+\kappa})$,

[TABLE]

Let us now recall that a Nemytskii operator $F:H\to H$ is defined such that $F(u)=f\circ u$ for all $u\in H$, where $f:\mathbb{R}\to\mathbb{R}$ is a real-valued mapping, assumed to be at least globally Lipschitz continuous. As explained above, even if $f$ is assumed of class $\mathcal{C}^{2}$ with bounded first and second order derivatives, the associated nonlinear Nemytskii operator $F$ does not satisfy the conditions of Assumption 1. Indeed, the appropriate generalization requires estimates in $L^{p}(0,1)$ spaces, for $p\in(2,\infty)$ (using Hölder inequality).

2.3. Wiener process

Let \bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)} be a probability space, equipped with a filtration \bigl{(}\mathcal{F}_{t}\bigr{)}_{t\geq 0} satisfying the usual conditions. The expectation operator is denoted by $\mathbb{E}[\cdot]$.

Let \left(\bigl{(}\beta_{p}(t)\bigr{)}_{t\geq 0}\right)_{p\in\mathbb{N}} denote a sequence of independent standard real-valued Wiener processes, and let \bigl{(}\epsilon_{p}\bigr{)}_{p\in\mathbb{N}} be a complete orthonormal system of $H$, and \bigl{(}\gamma_{p}\bigr{)}_{p\in\mathbb{N}} be a sequence of nonnegative real numbers. The cylindrical Wiener process is defined as

[TABLE]

The $Q$-Wiener process is then defined as

[TABLE]

and can be written as $W^{Q}(t)=Q^{\frac{1}{2}}W(t)$, where the linear self-adjoint operators $Q^{\frac{1}{2}}$ and $Q$ satisfy

[TABLE]

Note that the $Q$-Wiener process $W^{Q}(t)$ takes values in $H$ if and only if $Q$ is a trace-class linear operator, i.e., ${\rm Tr}(Q)=\|Q^{\frac{1}{2}}\|_{\mathcal{L}_{2}(H)}^{2}=\sum_{p\in\mathbb{N}}\gamma_{p}<\infty$.

Assumption 2 states the conditions on $Q$ required to ensure the well-posedness of (1).

Assumption 2**.**

Assume that there exists $\alpha>0$ such that $\|(-A)^{\alpha-\frac{1}{2}}Q^{\frac{1}{2}}\|_{\mathcal{L}_{2}(H)}<\infty$.

Define the parameter $\overline{\alpha}$ as follows:

[TABLE]

then $\overline{\alpha}>0$ if and only if Assumption 2 is satisfied. For instance, if $Q=I$ (space-time white noise), then $\overline{\alpha}=\frac{1}{4}$. If $Q$ is a trace-class operator, then $\overline{\alpha}=\frac{1}{2}$. In the article, we are mostly interested in the regime $\overline{\alpha}\in(0,\frac{1}{2}]$.

The numerical experiments are performed with the following example: for all $p\in\mathbb{N}$, $\epsilon_{p}=e_{p}$ (thus the operators $A$ and $Q$ commute), and $\gamma_{p}=\lambda_{p}^{\frac{1}{2}-2\overline{\alpha}}$ (observe that in that case the notation is consistent with the definition of $\overline{\alpha}$ in the general case).

To conclude this section, let us introduce the following notation: if $U$ is a $H$-valued random variable, for all $\alpha\in[0,1]$, and $q\in\mathbb{N}$,

[TABLE]

2.4. Well-posedness and regularity properties

Solutions of (1) are understood in the mild sense: for all $t\geq 0$,

[TABLE]

Under the assumptions stated above, this problem is globally well-posed. We quote without proof the following standard result.

Proposition 2.2**.**

Let Assumption 2 be satisfied. For any initial condition $u_{0}\in H$, there exists a unique mild solution (2) of the SPDE (1). Moreover, for any $T>0$ and $q\in\mathbb{N}$, there exists $C_{T,q}\in(0,\infty)$ such that

[TABLE]

Moreover, for any \alpha\in\bigl{(}0,\min(\overline{\alpha},\frac{1}{2})\bigr{)}, there exists $C_{T,q,\alpha}\in(0,\infty)$ such that

[TABLE]

3. Parareal algorithms

3.1. Fine and coarse integrators

Let $T\in(0,\infty)$ be given. Introduce the so-called coarse and fine time-step sizes $\Delta T$ and $\delta t$. It is assumed that $T=N\Delta T$ and $\Delta T=J\delta t$, where $N$ and $J$ are integers. For all $n\in\{0,\ldots,N\}$ and $j\in\{0,\ldots,J\}$, let

[TABLE]

Note that the coarse and the fine integrators introduced below are random mappings. Precisely, for all $n\in\mathbb{N}$, the mappings $\mathcal{G}_{n}=\mathcal{G}(\cdot,t_{n},t_{n+1})$ and $\mathcal{F}_{n}=\mathcal{F}(\cdot,t_{n},t_{n+1})$ depend on the Wiener increments \bigl{(}W^{Q}(t)-W^{Q}(t_{n})\bigr{)}_{t_{n}\leq t\leq t_{n+1}}.

3.1.1. Coarse integrator

The coarse integrator is a numerical scheme with time-step size $\Delta T$. In this work, it is defined as follows: for all $n\in\{0,\ldots,N-1\}$ and all $u\in H$, let

[TABLE]

with

• •

either $\hat{S}_{\Delta T}=e^{\Delta TA}$ (exponential Euler scheme),

• •

or $\hat{S}_{\Delta T}=S_{\Delta T}=(I-\Delta TA)^{-1}$ (linear implicit Euler scheme).

The notation $\mathcal{G}_{\rm expo}$, resp. $\mathcal{G}_{\rm imp}$, is often used below, to refer to the coarse integrator with the exponential Euler scheme, resp. with the linear implicit Euler scheme. As will be seen below, these two coarse integrators have very different behaviors when applied to the SPDEs considered in this article.

3.1.2. Fine integrator

The fine integrator consists of $J$ steps of a numerical scheme with time-step size $\delta t$. In this work, this numerical scheme is obtained by the exponential Euler scheme. More precisely, introduce the auxiliary integrator $\mathcal{F}_{\rm aux}$: for all $n\in\{0,\ldots,N-1\}$, all $j\in\{0,\ldots,J-1\}$ and all $u\in H$,

[TABLE]

The fine integrator $\mathcal{F}$, at the coarse time scale, is then defined as follows:

[TABLE]

In other words, the solution $v_{n,J}=\mathcal{F}(u,t_{n},t_{n+1})$. is computed using the following recursion formula

[TABLE]

for $j\in\{0,\cdots,J-1\}$ with $J\delta t=\Delta T$.

3.2. The parareal algorithm

The initialization step of the parareal algorithm consists in applying the coarse integrator: for $n\in\{0,\ldots,N-1\}$

[TABLE]

Let $K\in\mathbb{N}$ denote the number of parareal iterations. Iterations for $k\in\{0,\ldots,K-1\}$ are defined as follows: given the values \bigl{(}u_{m}^{(k)}\bigr{)}_{0\leq m\leq N} at iteration $k$, then compute, for all $n\in\{0,\ldots,N-1\}$,

[TABLE]

The core of the approach lies in the ability, at each iteration in $k$, to perform in parallel the computations in (5) for different values of $n$, hence the terminology of “parareal algorithms” for “parallelization in real time”.

3.3. The reference solution

The reference solution is defined using the fine integrator:

[TABLE]

Observe that following (4), this reference solution is in fact defined in terms of the integrator $\mathcal{F}_{\rm aux}$, applied with the fine time-step size $\delta t$. Precisely, for all $n\in\{0,\ldots,N-1\}$, one has $u_{n+1}^{\rm ref}=v_{n,J}^{\rm ref}$ defined by

[TABLE]

This may be rewritten as follows: $u_{n}^{\rm ref}=v_{nJ}^{\rm ref}$, where for all $\ell\in\{0,\ldots,NJ-1\}$,

[TABLE]

In addition, note that

[TABLE]

To conclude this section, we state without proof two standard results concerning, first, the qualitative properties of the reference solution (moment estimates), second, the rate of convergence of the error $u_{n}^{\rm ref}-u(t_{n})$.

Proposition 3.1**.**

Let $T>0$, \alpha\in\bigl{(}0,\min(\overline{\alpha},\frac{1}{2})\bigr{)} and $q\in\mathbb{N}$. There exists $C_{T,q,\alpha}\in(0,\infty)$ such that for all $u_{0}\in D((-A)^{\alpha})$,

[TABLE]

and such that the following error estimate holds true:

[TABLE]

Remark 3.2**.**

The assumption that $u_{0}\in D((-A)^{\alpha})$ may be weakened using Proposition 2.1.

Remark 3.3**.**

Since noise is additive in the SPDE (1), the order of convergence in Proposition 3.1 may be larger than $\frac{1}{2}$ when $\overline{\alpha}$ is sufficiently large. This type of estimate is not considered in this article.

3.4. Error and residual operators

In the implementation of parareal algorithms, the reference solution $u_{n}^{\rm ref}$ defined above is not computed in practice. Instead, the quantity $u_{n}^{(k)}$ defined in (5) is computed. To estimate the error between $u_{n}^{(k)}$ and $u(t_{n})$, due to Proposition 3.1, it is sufficient to study the error between $u_{n}^{(k)}$ and $u_{n}^{\rm ref}$. Observe that the parareal solution $u_{n}^{(k)}$ may be computed with a lower computational cost than the reference solution $u_{n}^{\rm ref}$, using parallel computations in (5) (expect for the initialization). Error estimates are required to determine the choice of time-step sizes $\Delta T$ and $\delta t$, and of the number of parareal iterations $K$, to achieve a given error criterion, with minimal computational cost.

For all $n\in\{0,\ldots,N\}$ and $k\in\{0,\ldots,K\}$ (where $K$ is the number of parareal iterations), let the error be defined by

[TABLE]

Note that $\epsilon_{0}^{(k)}=0$ for all $k\in\mathbb{N}_{0}$.

Moreover, by construction, one get $\epsilon_{n}^{(k)}=0$ for all $n\leq k$, which indicates that the numerical solution $\{u_{n}^{(k)}\}_{n=0,\cdots,N}$ will definitely converge to the reference solution $\{u_{n}^{\rm ref}\}_{n=0,\cdots,N}$, if the iterated number $K$ is sufficiently large, i.e. $K\geq N$. However, to get a speedup, in practice, the iterated number $K$ will be chosen significantly smaller than $N$, which is further discussed in Section 3.5.

It is convenient to introduce the residual operators defined by

[TABLE]

for all $n\in\{0,\ldots,N-1\}$, where the notation $\mathcal{F}_{n}(u):=\mathcal{F}(u,t_{n},t_{n+1})$ and $\mathcal{G}_{n}(u):=\mathcal{G}(u,t_{n},t_{n+1})$ is used.

Then the error defined by (7) satisfies the recursion formula, where the residual operators $\mathcal{R}_{n}$ defined by (8) appear:

[TABLE]

where the linear operator $\hat{S}_{\Delta T}$ depends on the choice of the coarse integrator, see (3).

Up to this point, the choice of the coarse integrator plays no role in the presentation. The major finding of this article is that the behavior of the parareal algorithm applied for SPDEs (1) differs when choosing the exponential Euler scheme or the linear implicit Euler scheme as the coarse integrator. Indeed, the theoretical results and the numerical experiments reveal that, as the number of parareal iterations $k$ increases, the evolution of the order of convergence of the error $\epsilon_{n}^{(k)}$ has a different behavior depending on the choice of coarse integrator.

3.5. Analysis of the computational cost

The objective of this section is to compare the costs for computing $u_{n}^{(k)}$ using the parareal algorithm 5, and for computing the reference solution $u_{n}^{\rm ref}$.

The computational advantage of using the parareal algorithm is due to the possibility to compute the quantities $u_{n}^{(k+1)}$ in parallel, for fixed $k\geq 0$, see (5). Let $N_{\rm proc}$ denote the number of available processors.

Let $T\in(0,\infty)$ denote the final time, and consider $n=N$ such that $N\Delta T=T$. Denote by $\tau_{\mathcal{G}}$ the computational time for one evaluation of $\mathcal{G}(u,t_{n},t_{n+1})$. It is assumed that $\tau_{\mathcal{G}}$ does not depend on $\Delta T$, $n\in\mathbb{N}_{0}$ and on $u\in H$. Denote also by $\tau_{\mathcal{F},{\rm aux}}$ denote the computational time for one realization of $\mathcal{F}_{\rm aux}(u,t_{n,j},t_{n,j+1})$. Then the computational time for one realization of $\mathcal{F}(\cdot,t_{n},t_{n+1})$, denoted by $\tau_{\mathcal{F}}$, is

[TABLE]

It is also assumed that $\tau_{\mathcal{F},{\rm aux}}(u,t_{n,j},t_{n,j+1})$ does not depend on $\delta t$, $n,j$ and on $u$.

3.5.1. Parareal algorithm

For the initialization step, the computational cost is equal to $N\tau_{\mathcal{G}}$, since at this stage no parallelization procedure is applied.

For each iteration of the algorithm, observe that in (5), the third term $\mathcal{G}(u_{n}^{(k)},t_{n},t_{n+1})$ has already been computed at the previous iteration, and that the values of the second term $\mathcal{F}(u_{n}^{(k)},t_{n},t_{n+1})$ may be computed in parallel. A sequential computation remains to be done, thus the computational cost of one iteration of the parareal algorithm is

[TABLE]

If $K$ iterations of the parareal algorithm are performed, the associated computational cost is thus equal to

[TABLE]

3.5.2. Reference solution

The reference solution $u_{n}^{\rm ref}$ is computed using the fine integrator $\mathcal{F}_{\rm aux}$ with time-step size $\delta t$, see (6). The associated computational cost is equal to

[TABLE]

3.5.3. Efficiency

The efficiency of the parareal algorithm, compared with a direct simulation using the fine integrator only, is thus studied in terms of the ratio

[TABLE]

Note that since the efficiency never goes to infinity as $\Delta T$ goes to [math], the parareal algorithm does not improve the rate of convergence with respect to the time-step size. Instead, the parareal algorithm may improve the computational efficiency. To go further in the analysis of the efficiency of the parareal algorithm, it is essential to study the rate of convergence of the error $\epsilon_{n}^{(k)}$ in terms of $\Delta T$.

On the one hand, assume that the rate of convergence does not depend on $k$. Then, to balance the errors $u_{n}^{(k)}-u_{n}^{\rm ref}$ and $u_{n}^{\rm ref}-u(n\Delta T)$, it is necessary to choose $\Delta T$ and $\delta t$ of the same size, hence

[TABLE]

To maximize the efficiency $\mathcal{E}$ above, the optimal choice is apparently $K=0$: parareal iterations increase the computational cost (linearly in $K$), in spite of the use of parallelization.

On the other hand, assume that the rate of convergence of $\epsilon_{n}^{(k)}$ is strictly larger than the rate for $\epsilon_{n}^{(0)}$. Then the time-step sizes are chosen such that $\delta t={\rm o}(\Delta T)$ as $\Delta T\to 0$, and the efficiency is then of size $\frac{N_{\rm proc}}{K}$. This means that parallelization ($N_{\rm proc}\geq 2$) reduces the cost, and that the optimal choice is $K=1$.

In Sections 4 and 5, the rates of convergence of the error $\epsilon_{n}^{(k)}$ with respect to $\Delta T$ are studied, depending on the choice of the coarse integrator, and on the regularity of the noise.

4. Linear-implicit Euler scheme as the coarse integrator

The objective of this section is to prove that, when the coarse integrator is chosen as the linear implicit Euler scheme, i.e. $\mathcal{G}=\mathcal{G}_{\rm imp}$ with $\hat{S}_{\Delta T}=S_{\Delta T}$, then the behavior of the parareal algorithm depends a lot on the regularity of the noise. More precisely, it is proved that the order of convergence of the error $\epsilon_{n}^{(k)}$ to [math] cannot exceed $\overline{\alpha}$, and in particular saturates when $k$ increases. Essentially, the order of convergence (in the framework studied below) is equal to $\min(\overline{\alpha},k+1)$. The theoretical and numerical results are consistent, and show that the rates obtained by the theoretical analysis are sharp.

For the theoretical analysis developped in this section, the framework is as follows. First, the initial condition $u_{0}=0$ and the nonlinear operator $F=0$ are set equal to [math]. Second, the covariance $Q$ of the noise commutes with $A$, i.e. $Qe_{p}=\gamma_{p}e_{p}$, for all $p\in\mathbb{N}$, and eigenvalues satisfy $\gamma_{p}=\lambda_{p}^{\frac{1}{2}-2\overline{\alpha}}$.

In this case, the solution is a Gaussian process, and $u(t)=\int_{0}^{t}e^{(t-s)A}dW^{Q}(s)$. Moreover, the recursion formula (9) for the error yields the equalities

[TABLE]

4.1. Theoretical error estimates

The main theoretical result of this section is Theorem 4.1.

Theorem 4.1**.**

Assume that $\gamma_{p}=\lambda_{p}^{\frac{1}{2}-2\overline{\alpha}}$, with $\overline{\alpha}>0$. Let $T\in(0,\infty)$ and $k\in\mathbb{N}_{0}$.

• •

If $k+1<\overline{\alpha}$, then there exists $C_{T,k,\overline{\alpha}}\in(0,\infty)$ such that

[TABLE]

• •

If $k+1\geq\overline{\alpha}$, then for all $\alpha\in(0,\overline{\alpha})$, there exists $C_{T,k,\alpha}\in(0,\infty)$ such that

[TABLE]

In particular, for $\overline{\alpha}=\frac{1}{4}$ ($Q=I$, space-time white noise), or $\overline{\alpha}=\frac{1}{2}$ ($Q$ is trace-class), the rate of convergence does not depend on $k$, and performing parareal iterations does not increase the order of convergence, see Section 3.5.

Proof.

For all $p\in\mathbb{N}$, let $\epsilon_{n}^{(k)}(p)=\langle\epsilon_{n}^{(k)},e_{p}\rangle$ denote the $p$-th component of the error $\epsilon_{n}^{(k)}$. Then the expression above is rewritten as

[TABLE]

where, for $z\in(-\infty,0]$, one has $V(z):=\frac{1}{1-z}$ and $R(z):=e^{z}-\frac{1}{1-z}$. The inequality $|R(z)|\leq 1\wedge|z|^{2}$ yields

[TABLE]

Observe that

[TABLE]

Thus for all $n$ and all $k$, one obtains (using a recursion argument)

[TABLE]

It remains to study the error $\epsilon_{n}^{(0)}$ at the initialization step. One has the identity

[TABLE]

Let us prove the following claim: there exists $C\in(0,\infty)$ such that

[TABLE]

Using the Itô isometry formula, and the fact $b^{m}-a^{m}\leq mb^{m-1}(b-a)$ for all $0\leq a\leq b$ and $m\in\mathbb{N}$, (with $b=V(-\lambda_{p}(\Delta T)$ and $a=e^{-\lambda_{p}\Delta T}$), one obtains

[TABLE]

This concludes the proof of the claim (10). Then, using the expression $\gamma_{p}=\lambda_{p}^{\frac{1}{2}-2\overline{\alpha}}$ of the eigenvalues of the covariance operator $Q$, one obtains

[TABLE]

It remains to identify the orders of convergence. On the one hand, assume that $k+1<\overline{\alpha}$. Using \bigl{(}1\wedge(\lambda_{p}\Delta T)\bigr{)}^{2(k+1)}\leq(\lambda_{p}\Delta T)^{2(k+1)} yields

[TABLE]

since $\overline{\alpha}-k-1>0$.

On the other hand, assume that $k+1\geq\overline{\alpha}$, and let $\alpha\in(0,\overline{\alpha})$. Using the inequality \bigl{(}1\wedge(\lambda_{p}\Delta T)\bigr{)}^{2(k+1)}\leq(\lambda_{p}\Delta T)^{2\alpha} yields

[TABLE]

This concludes the proof of Theorem 4.1. ∎

4.2. Numerical experiments

The objective of this section is to demonstrate that the orders of convergence obtained in Theorem 4.1 are sharp. In addition, experiments in the semilinear case ($F\neq 0$) are also provided.

First, the SPDE $du(t)=Au(t)+dW^{Q}(t)$, with $u(0)=0$, is considered, where the covariance $Q$ is given as above ($Qe_{p}=\gamma_{p}e_{p}$, with $\gamma_{p}=\lambda_{p}^{\frac{1}{2}-2\overline{\alpha}}$). Spatial discretization is performed using finite differences, with mesh size $h=0.01$. In addition, the noise is truncated, i.e. the $Q$-Wiener process $W^{Q}(t)$ is replaced by $\sum_{p=1}^{P}\gamma_{p}^{\frac{1}{2}}\beta_{p}(t)e_{p}$, with $P=100$. Numerical parameters are chosen as follows: the final time is $T=1$, the fine time-step size is $\delta t=2^{-13}$, and the coarse time-step size is $\Delta T=J\delta t$ with $J=2^{j}$, $j=4,\cdots,9$. An average over $M=100$ independent Monte-Carlo samples is used to approximate the expectations.

First, Figure 1 reports numerical simulations for $\overline{\alpha}=0.25$ (space-time white noise) and $\overline{\alpha}=0.5$ (trace-class noise). For all values of $k$, the order of convergence is equal to $\overline{\alpha}$, as predicted by Theorem 4.1.

Second, Figures 2 (two fixed values of $k$ and $\overline{\alpha}$ varies) and 3 (two fixed values of $\overline{\alpha}$ and $k$ varies) allow us to check that the orders of convergence in Theorem 4.1 are sharp: the order is indeed equal to $\min(\overline{\alpha},k+1)$.

To conclude this section, we report numerical simulations in the semilinear case.

Figures 4 and 5 show the order for semilinear equation (1) with $F(u)=\cos(u)$ and $F(u)=5\cos(u)$, respectively. The order for the additive noise case with $\overline{\alpha}=4$ (on the left) is limited to $\frac{3}{2}$ when $k\geq 2$, which is the same as the deterministic case (on the right).

5. Exponential Euler scheme as the coarse integrator

The objective of this section is to prove that, when the exponential Euler scheme is chosen as the coarse integrator, i.e. $\mathcal{G}=\mathcal{G}_{\rm expo}$ with $\hat{S}_{\Delta T}=e^{\Delta TA}$, then parareal iterations improve the rate of convergence of the error $\epsilon_{n}^{(k)}$ to [math], in terms of $\Delta T$. Contrary to the situation of Section 4, this effect holds true without restrictions on the regularity parameter $\overline{\alpha}$, in particular for $\overline{\alpha}=\frac{1}{4}$ (space-time white noise).

The analysis in this section is performed for the SPDE (1), with the nonlinear coefficient $F$ satisfying Assumption 1. Let also Assumption 2 be satisfied.

The content of this section is organized as follows. The main results of this sections are the error estimates stated in Theorems 5.2 – which gives an order of convergence $(k+1)\min(\overline{\alpha},\frac{1}{2})$ for all $k\in\mathbb{N}_{0}$ – and 5.3 – which gives an improved order of convergence when $k\geq 2$, see Section 5.1. Numerical experiments in Section 5.2 illustrate that the result in Theorem 5.2 is sharp when $k=0$ and $k=1$, and that indeed better convergence rates are obtained for $k\geq 2$. Proofs of the results are provided in Section 5.3, based on auxiliary results which are proved in Section 5.4.

Observe that, when the coarse integrator is the exponential Euler scheme, then the recursion formula (9) for the error yields the equality

[TABLE]

where we recall that the residual operator $\mathcal{R}_{n}$ are defined by (8). As will be clear below, if $F=0$ then $\epsilon_{n}^{(k)}=0$ for all $n$, as soon as $k\geq 1$. This property reveals why the choice of the exponential Euler scheme as the coarse and the fine integrator provides better results.

Note also that the fact that the noise is additive in (1) is fundamental in the analysis. In addition, since noise is additive, it is expected that the order of convergence for $k=0$ may be equal to $1$ (instead of $\frac{1}{2}$) if $\overline{\alpha}$ is sufficiently large (at least larger than $\frac{1}{2}$). This effect is not considered below, since the objective is mainly to study the increase in the order of convergence produced by parareal iterations, and in particular in situations where the noise is not very regular, i.e. for space-time white noise, with $\overline{\alpha}=\frac{1}{4}$.

5.1. Statement of error estimates

For the analysis, it is important to first state moment bounds for the solution $u_{n}^{(k)}$, which are similar to Propositions 2.2 and 3.1 for the exact solution $u(t_{n})$ and the reference solution $u_{n}^{\rm ref}$ respectively.

Proposition 5.1**.**

For $T\in(0,\infty)$, $k\in\mathbb{N}_{0}$, $q\in\mathbb{N}$. There exists $C_{T,k,q}\in(0,\infty)$ such that, for all $u_{0}\in H$ and $\Delta T\in(0,1)$,

[TABLE]

Moreover, let $\alpha\in[0,\min(\overline{\alpha},\frac{1}{2}))$. There exists $C_{T,k,q,\alpha}\in(0,\infty)$ such that, for all $u_{0}\in D((-A)^{\alpha})$ and $\Delta T\in(0,1)$,

[TABLE]

The first error estimates are stated in Theorem 5.2, which may be interpreted as follows: each parareal iteration improves the rate of convergence, proportionally to $\min(\overline{\alpha},\frac{1}{2})$. Let us stress that Theorem 5.2 is optimal for $k=0$ and $k=1$, as illustrated by the numerical experiments reported in Section 5.2.

Theorem 5.2**.**

Let $T\in(0,\infty)$, $k\in\mathbb{N}_{0}$, $q\in\mathbb{N}$, $\alpha\in[0,\min(\overline{\alpha},\frac{1}{2}))$, and arbitrarily small $\kappa>0$. For all $u_{0}\in D((-A)^{\alpha})$, there exists $C_{T,k,q,\alpha,\kappa}(u_{0})\in(0,\infty)$ such that, for all $\Delta T\in(0,1)$,

[TABLE]

Following the discussion in Section 3.5 concerning the computational cost, Theorem 5.2 shows that applying the parareal algorithm may be used to reduce the computational cost, if the coarse integrator is the exponential Euler scheme, whatever the regularity of the noise.

The second result, Theorem 5.3, states that the estimate from Theorem 5.3 can be improved when $k\geq 2$. The practical relevance of this result is questionable: it requires $k\geq 2$, whereas it is expected (see Section 3.5) that choosing $k=1$ is optimal. Nevertheless, the study of the phenomenon stated in Theorem 5.3 is motivated by the numerical experiments reported in Section 5.2, which exhibit that indeed the order of convergence is larger than $(k+1)\min(\overline{\alpha},\frac{1}{2})$, namely it is equal at least equal to $(k-1)\min(2\overline{\alpha},\frac{1}{2})+2\alpha\leq k\min(2\overline{\alpha},\frac{1}{2})$.

Theorem 5.3**.**

Let $T\in(0,\infty)$, $k\in\mathbb{N}\setminus\{1\}$, $q\in\mathbb{N}$, $\alpha\in(0,\min(\overline{\alpha},\frac{1}{2}))$, and arbitrarily small $\kappa>0$. For all $u_{0}\in D((-A)^{\alpha})$, there exists $C_{T,k,q,\alpha,\kappa}(u_{0})\in(0,\infty)$ such that, for all $\Delta T\in(0,1)$,

[TABLE]

The proofs of Proposition 5.1 and Theorems 5.2 and 5.3 are postponed to Section 5.3.

To simplify the exposition, the way the constants $C_{T,k,q,\alpha,\kappa}(u_{0})$ above depend on $|u_{0}|_{\alpha}$ is not made precise.

Remark 5.4**.**

Theorem 5.2 (and Lemma 5.5–$({\rm i})$) can be proved under less restrictive conditions on the linear operator $F$, instead of Assumption 1: there exists $\eta\in(0,\frac{1}{2})$ such that one has estimates of the type

[TABLE]

and

[TABLE]

This setting encompasses the case of Nemytskii operators, with $\eta\in(\frac{1}{4},\frac{1}{2})$.

5.2. Numerical experiments

The objective of this section is to illustrate Theorem 5.2 and 5.3.

The SPDE (1), with nonlinear operator $F(u)=5\cos(u)$ and initial condition $u(0)=0$, is considered, with covariance operator given by $Qe_{p}=\gamma_{p}e_{p}$, with $\gamma_{p}=\lambda_{p}^{\frac{1}{2}-2\overline{\alpha}}$. Spatial discretization is performed using finite differences, with mesh size $h=0.01$. In addition, the noise is truncated, i.e. the $Q$-Wiener process $W^{Q}(t)$ is replaced by $\sum_{p=1}^{P}\gamma_{p}^{\frac{1}{2}}\beta_{p}(t)e_{p}$, with $P=100$. Numerical parameters are chosen as follows: the final time is $T=1$, the fine time-step size is $\delta t=2^{-15}$, and the coarse time-step size is $\Delta T=J\delta t$ with $J=2^{j}$, $j=5,\cdots,10$. An average over $M=100$ independent Monte-Carlo samples is used to approximate the expectations.

Figure 6 demonstrates that the orders of convergence in Theorem 5.2 are sharp when $k=0$ and $k=1$, for different values of $\overline{\alpha}$.

Figure 7 then illustrates that for $k\geq 2$, the orders of convergence in Theorem 5.2 are not sharp for $\overline{\alpha}=0.25$ and $\overline{\alpha}=0.4$. In fact, even the improved error estimates from Theorem 5.3 seem to be sub-optimal in the numerical experiment, especially when $\overline{\alpha}=0.4$ (figure on the right) – in this case $\min(2\alpha,\frac{1}{2})=\frac{1}{2}$. When $\overline{\alpha}=0.25$ (figure on the left), the improved theoretical rate in Theorem 5.3 is closer to the observed numerical rate.

5.3. Proof of error estimates

The objective of this section is to provide proofs for the results stated in Section 5.1. In fact, they are all based on the following Lemma concerning Lipschitz continuity properties, in appropriate norms, of the residual operators.

Lemma 5.5**.**

Let Assumptions 1 and 2 hold, and $\alpha\in[0,1)$.

$({\rm i})$* There exists $C_{F,\alpha}\in(0,\infty)$ such that for all $\Delta T\in(0,1]$ and all $u_{1},u_{2}\in H$, one has*

[TABLE]

$({\rm ii})$* For all $q\in\mathbb{N}$, $\alpha\in(0,\min(\overline{\alpha},\frac{1}{2}))$, $\beta\in[0,\alpha)$ and arbitrarily small \kappa\in\bigl{(}0,\min(\alpha-\beta,1-\alpha-\beta)\bigr{)}, there exists $C_{T,q,\alpha,\beta,\kappa}\in(0,\infty)$ such that for all $\Delta T\in(0,1)$ and all $u_{1},u_{2}\in D((-A)^{\alpha})$, one has*

[TABLE]

The proof of Lemma 5.5 is postponed to Section 5.4.

Observe that Statement $({\rm i})$ (in the case $\alpha=0$) is not sufficient to exhibit a positive rate of convergence $\gamma$ for the error $\epsilon_{n}^{(k)}$, in terms of $\Delta T$. Indeed, that would required the Lipschitz constant of the residual operator to be of size $\Delta T^{1+\gamma}$. Statement $({\rm i})$, with $\alpha=0$, is used only to establish stability properties of the algorithm, i.e. moment bounds for $u_{n}^{(k)}$ or to apply a Gronwall lemma. With $\alpha>0$, this statement is employed to establish moment bounds in $|\cdot|_{\alpha}$ norms, and to deal with some terms for which regularization properties of the semigroup cannot be used.

Statement $({\rm ii})$ is the key result in order to get the rates of convergence given in Theorem 5.2 and 5.3. It requires to choose appropriate norms to analyze the Lipschitz constant of the residual operators: estimation in a weaker norm for $\alpha>0$, or with higher regularity for $\beta>0$. More precisely, Theorem 5.2 follows from Statement $({\rm ii})$ with $\beta=0$, and then Theorem 5.3 follows from Statement $({\rm ii})$ with $\beta>0$.

Finally, observe that the need to choose appropriate norms is specific to the infinite dimensional situation. Moreover, Lemma 5.5 is not satisfied if the linear implicit Euler scheme is chosen as the coarse integrator, even if $F=0$: indeed, in that case, \mathcal{R}_{n}(u_{2})-\mathcal{R}_{n}(u_{1})=\bigl{(}S_{\Delta T}-e^{\Delta T}\bigr{)}(u_{2}-u_{1}), and $\|S_{\Delta T}-e^{\Delta T}\|_{\mathcal{L}(H)}$ does not converge to [math] as $\Delta T\to 0$. This observation explains why the behaviors and the analysis of the parareal algorithm depends a lot on the choice of the coarse integrator for parabolic semilinear SPDEs.

Remark 5.6**.**

Note that the estimates stated in Lemma 5.5 require $u_{1}$ and $u_{2}$ to be deterministic, however below they are applied to random elements, which are measurable with respect to \sigma\bigl{(}W(t);t\leq t_{n}\bigr{)}, whereas the noise component in the residual operator $\mathcal{R}_{n}$ is measurable with respect to \sigma\bigl{(}W(t);t_{n}\leq t\leq t_{n+1}\bigr{)}. Applying a straightforward conditioning argument, and interpreting the expectation in Lemma 5.5 as a conditional expectation, yield the required estimates below.

Proof of Proposition 5.1.

First, moment bounds for the reference solution $u_{n}^{\rm ref}$ are provided by Proposition 3.1. In addition, in this section, the initialization step of the parareal algorithm consists in applying the fine integrator (which is the exponential Euler scheme) with time-step size $\Delta T$, thus the result of Proposition 3.1 also applies to obtain moment bounds for $u_{n}^{(0)}$.

Since $u_{n}^{(k)}=\epsilon_{n}^{(k)}+u_{n}^{\rm ref}$, it thus only remains to prove moment bounds for $\epsilon_{n}^{(k)}$, when $k\geq 1$.

Using the expression of the error (11), and the global Lipschitz continuity of $F$, moments are treated as follows: for all $q\in\mathbb{N}$ and $\alpha\in[0,1)$,

[TABLE]

Assume first that $\alpha=0$. Then the claim follows from the application of Lemma 5.5–$({\rm i})$, the use of the discrete Gronwall Lemma, and from the use of a recursion argument with respect to $k$. When $\alpha>0$, it remains to apply Lemma 5.5–$({\rm i})$, and to apply the regularization estimate from Proposition 2.1 to conclude the proof. ∎

Proof of Theorem 5.2.

First, assume that $k=0$. The claim then follows from Proposition 3.1, since both the coarse and the fine integrators are given by the exponential Euler scheme, thus \epsilon_{n}^{(k)}=\bigl{(}u_{n}^{(k)}-u(t_{n})\bigr{)}-\bigl{(}u_{n}^{\rm ref}-u(t_{n})\bigr{)}, and the initialization of the parareal algorithm consists in applying the coarse integrator.

Let now $k\geq 1$. Owing to the expression of the error (11), applying Lemma 5.5–$({\rm ii})$ with $\beta=0$ and Hölder inequality (see Remark 5.6 for the conditioning argument), one obtains

[TABLE]

Due to the moment estimates of Proposition 5.1, applying the discrete Gronwall lemma yields

[TABLE]

Then a straightforward recursion argument yields

[TABLE]

owiing to the result when $k=0$. This concludes the proof of Theorem 5.2. ∎

Proof of Theorem 5.3.

The proof consists of two steps.

Step $1$. If $k\geq 1$, the error evaluated in the $|\cdot|_{\alpha}$ norm is also of the order $(k+1)(\alpha-\kappa)$:

[TABLE]

Indeed, using the expression of the error (11), Proposition 2.1, and the two statements of Lemma 5.5, one obtains

[TABLE]

owing to Theorem 5.2. With the assumptions $\alpha\leq\frac{1}{2}$, it is straightforward to check that $1-\alpha+k(\alpha-\kappa)\geq(k+1)(\alpha-\kappa)$. This concludes Step $1$.

Step $2$. It remains to establish the error estimate of Theorem 5.3.

Using the expression of the error (11), Proposition 2.1, the two statements in Lemma 5.5, with $\beta=\alpha-\kappa$, and the moment bounds from Proposition 5.1, one obtains

[TABLE]

using the fact that $1-\alpha\geq\frac{1}{2}\geq\min(2\alpha,\frac{1}{2})-2\kappa$.

Applying the discrete Gronwall lemma and a recursion argument,

[TABLE]

owing to Step $1$ with $k=1$.

This concludes the proof of Theorem 5.3. ∎

Remark 5.7**.**

The claim in Step $1$ in the proof of Theorem 5.3 requires $k\geq 1$, and is not correct when $k=0$. In fact, when $k\geq 1$, in the expression of the error, contributions of stochastic terms, which would have low regularity properties, do not appear explicitly since noise is additive.

The application of the claim of Step $2$ explains why Theorem 5.3 holds true only if $k\geq 2$.

5.4. Proof of Lemma 5.5

The proof of Lemma5.5 consists in proving bounds on the derivative $D\mathcal{R}_{n}(u).h=D\mathcal{F}_{n}(u).h-DG_{n}(u).h$ (recall the residual operator $\mathcal{R}_{n}$ is defined by (8)). Let us introduce notation which is used below. On the one hand, due to (3), it is straightforward to compute

[TABLE]

On the other hand, the derivative of $\mathcal{F}_{n}$ is computed as follows: one has $D\mathcal{F}_{n}(u).h=\eta_{n,J}^{h}$, where $J\delta t=\Delta T$, and one has the following recursion formulae

[TABLE]

Then a straightforward computation yields

[TABLE]

Observe that the same term $e^{\Delta TA}h$ appears in expressions of the derivatives $D\mathcal{F}_{n}(u).h$ and $D\mathcal{G}_{n}(u).h$. This term thus does not appear in $D\mathcal{R}_{n}(u)$. This is a crucial property (which is not satisfied when the coarse integrator is the linear implicit Euler scheme). Observe also that $\mathcal{R}_{n}=0$ if $F=0$, since the linear part is solved exactly.

Proof of Lemma 5.5–$({\rm i})$.

Since $F$ is globally Lipschitz continuous, applying the discrete Gronwall Lemma yields the following (almost sure) inequality:

[TABLE]

Owing to Proposition 2.1 and using the global Lipschitz continuity property of $F$, it is straightforward to check that

[TABLE]

using the inequality $\delta t\sum_{j=0}^{J-1}\frac{1}{(J-j)\delta t)^{\alpha}}\leq C_{\alpha}\Delta T^{1-\alpha}$. This concludes the proof. ∎

Observe that in the proof above, only the Lipschitz continuity property for $F$ was employed. In the proof below, the estimates stated in Assumption 1 are crucial.

Proof of Lemma 5.5–$({\rm ii})$.

Using the notation above, and expressions of $D\mathcal{F}_{n}(u).h$ and of $D\mathcal{G}_{n}(u).h$, the derivative $D\mathcal{R}_{n}(u).h=D\mathcal{F}_{n}(u).h-D\mathcal{G}_{n}(u).h$ is decomposed as follows:

[TABLE]

The first term is treated as follows: using Assumption 1,

[TABLE]

since $j\delta t\leq J\delta t=\Delta T$, and $\alpha+\beta\leq 2\alpha\leq 1$.

The second term is treated as follows: observe that $\eta_{n,0}^{h}=h$, thus, using Assumption 1,

[TABLE]

Note that the increment $\eta_{n,j}^{h}-\eta_{n,0}^{h}$ satisfies

[TABLE]

with the inequality

[TABLE]

Using the inequality $|\eta_{n,j}^{h}|\leq C|h|$, and using that $\kappa$ is chosen such that $\alpha+\beta+\kappa<1$, one obtains

[TABLE]

It remains to treat the third term: using $\beta+\kappa\leq\alpha$, the last inequality from Assumption 1 yields

[TABLE]

with $|\eta_{n,j}^{h}|_{\beta}\leq C|h|_{\beta}$. Note that the following moment bound holds true: there exists $C_{q,\beta}\in(0,\infty)$ such that $\vvvert v_{n,j}\vvvert_{q,\beta}\leq C_{q,\beta}$. Moreover, the proof of the following regularity result is straightforward:

[TABLE]

see Proposition 2.2 for a similar estimate when $\beta=0$, for the process \bigl{(}u(t)\bigr{)}_{t\geq 0}. Then applying Hölder inequality yields

[TABLE]

Gathering the estimates then concludes the proof. ∎