A note on stochastic exponential integrators for the finite element discretization of semilinear SPDEs driven by additive noise

# A note on stochastic exponential integrators for the finite element discretization of semilinear SPDEs driven by additive noise

Antoine Tambue The African Institute for Mathematical Sciences(AIMS) of South Africa and Stellenbosch University, 6-8 Melrose Road, Muizenberg 7945, South Africa Center for Research in Computational and Applied Mechanics (CERECAM), and Department of Mathematics and Applied Mathematics, University of Cape Town, 7701 Rondebosch, South Africa. Jean Daniel Mukam Fakultät für Mathematik, Technische Universität Chemnitz, 09126 Chemnitz, Germany
###### Abstract

This paper aims to investigate the strong convergence analysis of the exponential Euler method for parabolic stochastic partial differential equation (SPDE) driven by additive noise under more relaxed conditions. The SPDE is discretized in space by the finite element method. In contrast to very restrictive assumptions made in the literature on the drift term to achieve optimal convergence orders with self-adjoint operator, we weaken those assumptions and assume only the standard Lipschitz condition on the drift and deal with not necessarily self-adjoint linear operator. Under this relaxed assumption, we have achieved optimal convergence orders, which depend on the regularity of the noise and the initial data. In particular, for trace class noise we achieve convergence orders of . These optimal convergence orders are due to a new optimal regularity result we have further derived here, which was not yet proved in the literature, even in the case of self-adjoint operator. Numerical experiments to illustrate our theoretical result are provided.

###### keywords:
Stochastic parabolic partial differential equations, Stochastic exponential integrators, Additive noise, Finite element method, Errors estimate.

## 1 Introduction

We consider the numerical approximation of SPDE defined in , , with initial value and boundary conditions of the following type

 dX(t)=[AX(t)+F(X(t))]dt+dW(t),X(0)=X0,t∈(0,T] (1)

on the Hilbert space , where is the final time, is a linear operator which is unbounded and not necessarily self-adjoint. The noise is a Wiener process defined in a filtered probability space and is positive and self-adjoint. The filtration is assumed to fulfill the usual conditions. Precise assumptions on , and will be given in the next section to ensure the existence of the unique mild solution of (1), which has the following representation (see Prato (); Prevot ())

 X(t)=S(t)X0+∫t0S(t−s)F(X(s))ds+∫t0S(t−s)dW(s),t∈[0,T]. (2)

Equations of type (1) are used to model different real world phenomena in different fields such as biology, chemistry, physics etc Shardlow (); SebaGatam (). In more cases, explicit solutions of SPDEs are unknown, therefore numerical methods are the only good tools for their approximations. Numerical approximation of SPDE of type (1) is therefore an active research area and have attracted a lot of attentions since two decades (see for example Antonio1 (); Xiaojie1 (); Xiaojie2 (); Raphael (); Jentzen1 (); Jentzen2 () and references therein). For many numerical schemes, the optimal convergence order in time under only standard Lipschitz conditions is (see for example Raphael (); Antonio1 ()). By incorporating suitable linear functionals of the noise, the accelerated exponential Euler (AEE) methods were proved in Jentzen1 (); Jentzen2 (); Xiaojie1 () to converge strongly to the mild solution of (1). Due to more information incorporated on the noise, the AEE schemes usually achieve convergence orders higher than . The convergence analysis is usually done under restrictive assumptions (see for example (Jentzen1, , Assumption 2.4)). The AEE method in Jentzen1 () was recently analyzed in Xiaojie1 () under more relaxed conditions. These assumptions were also used in Xiaojie2 (); Antjd1 () for implicit scheme and exponential integrators with standard Brownian increments to obtain optimal convergence orders greater than . However assumptions made on the drift function in Xiaojie1 (); Xiaojie2 (); Antjd1 () are still restrictive as they involve first and sometime second order derivatives of the drift function. In many problems the nonlinear function may not be differentiable. An illustrative example is the function , , which is not differentiable at . In this paper, we investigate the strong convergence analysis of standard stochastic exponential integrators Antonio1 () using only standard Lipschitz condition on the drift function. The result indicates how the convergence orders depend on the regularity of the initial data and the noise. In fact, we achieve optimal convergence order , where is the regularity’s parameter of the noise and initial data (see Assumption 2.1). Let us mention that even under restrictive assumptions (namely (Jentzen1, , Assumption 2.4) and (Jentzen2, , Assumption 2)), the convergence orders in time of the AEE schemes analyzed in Jentzen1 (); Jentzen2 () are sub-optimal and have a logarithm reduction. It is worth to mention that the convergence order in time of the AEE scheme studied in Xiaojie1 () is optimal and does not have any logarithmic reduction. This is due to the sharp integral estimate (Stig1, , Lemma 3.2 (iii)) and the strong regularity assumptions on the drift term, namely (Xiaojie1, , Assumption 2.1). Note that (Stig1, , Lemma 3.2 (iii)) is only valid for self-adjoint operator. In this paper, we also extend (Stig1, , Lemma 3.2 (iii)) to the case of not necessarily self-adjoint operator in Lemma 2.1 and we further prove a new optimal regularity result, namely (20), which is not yet proved in the literature to the best of our knowledge. Lemma 2.1 and the regularity result (20) of Theorem 2.1 are therefore key ingredients to achieve optimal convergence orders in both space and time in this work.

The rest of the paper is organized as follows. Section 2 deals with the well posedness, the regularity result, the numerical scheme and the main result. In Section 3, we provide the proof of the main result. In Section 4, we have expanded the discussion on the nonlinear term and the associated challenges. We end the paper in Section 5 by giving some numerical experiments to sustain our theoretical result.

## 2 Mathematical setting and numerical method

### 2.1 Main assumptions and well posedness problem

Let us define functional spaces, norms and notations that will be used in the rest of the paper. Let be a separable Hilbert space, we denote by the Banach space of all equivalence classes of squared integrable -valued random variables, by the space of bounded linear mappings from to endowed with the usual operator norm , by the space of Hilbert-Schmidt operators from to . We equip with the norm

 ∥l∥2L2(U,H):=∞∑i=1∥lψi∥2,l∈L2(U,H), (3)

where is an orthonormal basis of . Note that (3) is independent of the orthonormal basis of . For simplicity, we use the notations and . In the sequel, we take and assume to be second order and is given by

 Au=d∑i,j=1∂∂xi(qij(x)∂u∂xj)−d∑i=1qi(x)∂u∂xi, (4)

where , . We assume that there exists such that

 d∑i,j=1qij(x)ξiξj≥c1|ξ|2,∀ξ∈Rd,x∈¯¯¯¯Λ.

As in Antonio1 (); Antjd1 (), we introduce two spaces and , such that , that depend on the boundary conditions for the domain of the operator and the corresponding bilinear form. For example, for Dirichlet boundary conditions, we take

 V=H=H10(Λ)={v∈H1(Λ):v=0on∂Λ}.

For Robin boundary conditions and Neumann boundary conditions, we take

 H={v∈H2(Λ):∂v/∂vA+α0v=0,on∂Λ},α0∈R.

Let be the bilinear operator associated to , which depends on the boundary condition, see e.g. Antonio1 (); Antjd1 () for details. Note that using the above assumptions on and we can prove that is sectorial and generates an analytic semigroup , the fractional powers of are also well defined and characterized Henry (); Antonio1 (); Antjd1 () for any by

where is the Gamma function, see Henry (). In addition, the following holds Henry (); Antonio1 (); Antjd1 ()

 (−A)αS(t)=S(t)(−A)αonD((−A)α),α≥0. (6)

In order to ensure the existence and the uniqueness of solution of (1) and for the purpose of the convergence analysis, we make the following assumption.

###### Assumption 2.1

We assume that the initial data , , the nonlinear function and the covariance operator satisfy the following properties

 ∥F(0)∥≤C,∥F(Y)−F(Z)∥ ≤ C∥Y−Z∥,Y,Z∈H, (7) ∥∥∥(−A)β−12Q12∥∥∥L2(H) < C, (8)

where is a positive constant.

The following lemma provides a version of (Stig1, , Lemma 3.2 (iii)) for not necessarily self-adjoint operator. Lemma 2.1 is useful to achieve optimal regularity result of the mild solution (2).

###### Lemma 2.1

For any and , there exists a constant such that

 ∫t2t1∥∥(−A)ρ/2S(t2−r)∥∥2L(H)dr≤C(t2−t1)1−ρ,0≤t1≤t2≤T, (9) ∫t2t1∥∥(−A)γ/2S(t2−r)∥∥L(H)dr≤C(t2−t1)1−γ/2,0≤t1≤t2≤T. (10)

Proof. Let us write , where and are respectively the self-adjoint and the non self-adjoint parts of . As in (Antonio2, , (147)), we use the Zassenhaus product formula Magnus (); Zassen2 () to decompose the semigroup as follows.

 S(t)=eAt=e(As+An)t=eAsteAnt∞∏k=2eCk, (11)

where are called Zassenhaus exponents. In (11), let us set

 SN(t):=eAnt∞∏k=2eCk,S(t)=Ss(t)SN(t), (12)

where is the semigroup generated by . Using the Baker-Campbell-Hausdorff representation formula Magnus (); Zassen1 (); Zassen2 (), one can prove as in Antonio2 () that is a linear bounded operator. Therefore using (12) and the boundedness of yields

 ∫t2t1∥∥(−A)ρ/2S(t2−r)∥∥2L(H)dr = ∫t2t1∥∥(−A)ρ/2Ss(t2−r)SN(t2−r)∥∥2L(H)dr (13) ≤ C∫t2t1∥∥(−A)ρ/2Ss(t2−r)∥∥2L(H)dr.

As in Antonio2 (), note that with equivalent norms (see Suzuki ()). Therefore by (Lions, , (3.3)) and by interpolation technique, it holds that , , with equivalent norms. It follows therefore that

 ∥∥(−A)ρ/2Ss(t2−r)∥∥L(H)≤C∥∥(−As)ρ/2Ss(t2−r)∥∥L(H). (14)

Substituting (14) in (13) yields

 ∫t2t1∥∥(−A)ρ/2S(t2−r)∥∥2L(H)dr≤C∫t2t1∥∥(−As)ρ/2Ss(t2−r)∥∥2L(H)dr. (15)

Since is self-adjoint, it follows from (Stig1, , Lemma 3.2 (iii)) that

 (16)

Substituting (16) in (15) complete the proof of (9). To prove (10), we use Cauchy-Schwarz inequality, (6) and (9). This yields

 ∫t2t1∥(−A)γ/2S(t2−r)∥L(H)dr = ∫t2t1∥∥∥(−A)γ/4S(12(t2−r))(−A)γ/4S(12(t2−r))∥∥∥L(H)dr (17) ≤ ∫t2t1∥∥∥(−A)γ/4S(12(t2−r))∥∥∥2L(H)dr ≤ C(t2−t1)1−γ/2.

This completes the proof of Lemma 2.1.

###### Theorem 2.1

Let Assumption 2.1 be satisfied. If is an -measurable -valued random variable, then there exists a unique mild solution of problem (1) represented by (2) and satisfying

 P[∫T0∥X(s)∥2ds<∞]=1, (18)

and for any , there exists a constant such that

 supt∈[0,T]E∥X(t)∥p≤C(1+E∥X0∥p). (19)

Moreover, the following optimal regularity estimate holds

 ∥∥(−A)β/2X(t)∥∥L2(Ω,H) ≤ C(1+∥(−A)β/2X0∥L2(Ω,H)),t∈[0,T]. (20)

Proof. The proof of the existence, the uniqueness and (19) can be found in (Prato, , Theorem 7.2). The proof of (20) for is similar to that of (Arnulf1, , Theorem 1) for multiplicative noise and can also be found in (Antjd1, , Lemma 3). Note that the case is of great importance in numerical analysis, it allows to avoid reduction of convergence order. It corresponds to (Arnulf1, , Theorem 1) with and (Stig1, , Theorem 3.1) with . To the best of our knowledge the case is not treated so far in the scientific literature, even for self-adjoint operators. We fill that gap with the help of Lemma 2.1. Indeed, from the mild solution (2), it follows by using triangle inequality that

 ∥∥(−A)β/2X(t)∥∥L2(Ω,H) ≤ ∥∥(−A)β/2X0∥∥L2(Ω,H)+∫t0∥∥(−A)β/2S(t−s)∥L(H)∥F(X(s))∥∥ds (21) + ∥∥∥∫t0(−A)β/2S(t−s)dW(s)∥∥∥L(H):=I0+I1+I2.

Using Lemma 2.1 we can easily prove that . Using Itô-isometry property, Lemma 2.1 and Assumption 2.1 yields

 I22 = ∫t0∥∥(−A)β/2S(t−s)Q1/2∥∥2L2(H)ds (22) ≤ ∫t0∥∥(−A)1/2S(t−s)∥∥2L(H)∥∥∥(−A)β−12Q1/2∥∥∥2L2(H)ds ≤ C∫t0∥∥(−A)1/2S(t−s)∥∥2L(H)ds ≤ C.

Substituting (22) and the estimate of in (21) complete the proof of the regularity result in (20).

### 2.2 Numerical scheme and main result

Let us first perform the space approximation of problem (1). We start by discretizing our domain by a finite triangulation. Let be a triangulation with maximal length . Let denote the space of continuous functions that are piecewise linear over the triangulation . We consider the projection and the discrete operator defined respectively from to and from to by

 (Phu,χ)=(u,χ),∀χ∈Vh,u∈H,and(Ahϕ,χ)=−a(ϕ,χ),∀ϕ,χ∈Vh. (23)

The discrete operator is also a generator of an analytic semigroup . The semi-discrete in space version of problem (1) consists of finding such that

 dXh(t)=[AhXh(t)+PhF(Xh(t))]dt+PhdW(t),Xh(0)=PhX0,t∈(0,T]. (24)

The mild solution of (24) has the following integral form

 Xh(t)=Sh(t)Xh(0)+∫t0Sh(t−s)PhF(Xh(s))ds+∫t0Sh(t−s)PhdW(s),t∈[0,T]. (25)

Let , where and . The mild solution of (1) at , can be also written as follows.

 X(tm)=S(Δt)X(tm−1)+∫tmtm−1S(tm−s)F(X(s))ds+∫tmtm−1S(tm−s)dW(s). (26)

For the time discretization, we consider the exponential Euler method proposed in Antonio1 (), which gives the numerical approximation of through the following recurrence

 Xhm=eAhΔtXhm−1+A−1h(eAhΔt−I)PhF(Xhm−1)+eAhΔtPh(Wtm+1−Wtm). (27)

with . The scheme (27) can be writen in the following integral form, useful for the error estimate

 Xhm=Sh(Δt)Xhm−1+∫tmtm−1Sh(tm−s)PhF(Xh(s))ds+∫tmtm−1Sh(Δt)PhdW(s). (28)

An equivalent formulation of (27) easy for simulation is given by

 Xhm=Xhm−1+PhΔWm−1+Δtφ1(ΔtAh)[Ah(Xhm−1+PhΔWm−1)+PhF(Xhm−1)], (29)

where

 φ1(ΔtAh)=(ΔtAh)−1(eΔtAh−I)=1Δt∫Δt0e(Δt−s)Ahds. (30)

With the numerical method in hand, we can state our strong convergence result.

Throughout this paper, is a positive constant independent of , , and and that may change from one place to another.

###### Theorem 2.2 (Main Result)

Let be the mild solution of problem (1) and the approximated solution at time with the scheme (27). If Assumption 2.1 is fulfilled, then the following error estimate holds

 ∥X(tm)−Xhm∥L2(Ω,H)≤C(hβ+Δtβ/2),m=0,1,⋯,M. (31)
###### Remark 2.1

Theorem 2.2 also holds if we replace the numerical approximation by the stochastic exponential Rosenbrock-Euler scheme proposed in Antjd1 (). Note that the arguments in this paper can be applied to the numerical schemes in Xiaojie1 (); Xiaojie2 () with the drift term satisfying only the global Lipschitz condition with the linear self-adjoint operator . In this case, the time error estimate in Xiaojie1 () will have the convergence rate 111Note that in our error estimates the deterministic part is of order . To update the proof in Xiaojie1 (), we just need to follow the deterministic part of the current work as the noise term is already high order in Xiaojie1 (). and that in Xiaojie2 () will have the same convergence rate as in (31). Note that in the case of not necessarily self-adjoint operator , the arguments in this paper are not enough to applied to the numerical scheme in Xiaojie2 (). In fact one need in addition some error estimates of the associated deterministic problem, which do not rely on the spectral decomposition of . This can be found in our accompanied paper Antjd2 ().

## 3 Preliminary results and proof of the main result

### 3.1 Preliminary results

The proof of the main result requires three important lemmas. Let us start by considering the following deterministic problem: find such that , , . The corresponding semi-discrete problem in space consists of finding such that , , . The following lemma can be found in (Antonio1, , Lemma 3.1) or (Antjd1, , Lemma 7).

###### Lemma 3.1

Let and . Then the following estimate holds

 ∥u(t)−uh(t)∥=∥(eAt−eAhtPh)v∥≤Chrt−(r−γ)/2∥v∥γ,v∈D((−A)γ/2). (32)
###### Lemma 3.2

Let Assumption 2.1 be fulfilled. Then the following error estimate holds

 ∥∥∥(−Ah)β−12PhQ12∥∥∥L2(H) ≤ C, (33) ∫t0∥(S(s)−Sh(s)Ph)v∥2ds ≤ Chγ∥v∥γ−1,v∈D((−A)(γ−1)/2),0≤γ≤2. (34)

Proof. The proof of (33) can be found in (Antjd1, , Lemma 11) and the proof of (34) can be found in (Antonio2, , Lemma 6.1 (ii)).

###### Lemma 3.3

(Discrete Gronwall lemma) Let and be sequences of non negatives numbers. Let be a non negative constant. If

 xk≤c+∑0≤i

then the following estimate holds

 xk≤c∏0≤j

Proof. As in Stephen (), we use the convention that an empty product is . Applying (Stephen, , Lemma 2) with and using the telescopic identity yields

 xk ≤ c+∑0≤i

where at the last step we used the inequality .

### 3.2 Proof of the main result

Let us now turn to the the proof of Theorem 2.2. Subtracting (28) from (26) and taking the norm in both sides yields

 ∥X(tm)−Xhm∥L2(Ω,H) ≤ ∥S(Δt)X(tm−1)−Sh(Δt)Xhm−1∥L2(Ω,H) (38) + ∫tmtm−1∥S(tm−s)F(X(s))−Sh(tm−s)PhF(Xhm−1)∥L2(Ω,H)ds + ∥∥∥∫tmtm−1(S(tm−s)−Sh(Δt)Ph)dW(s)∥∥∥L2(Ω,H) := II1+II2+II3.

Using triangle inequality, (20) and Lemma 3.1 with , it holds that

 II1 ≤ ∥(S(Δt)−Sh(Δt)Ph)X(tm−1)∥L2(Ω,H)+∥Sh(Δt)(X(tm−1)−Xhm−1)∥L2(Ω,H) (39) ≤ Chβ+C∥X(tm−1)−Xhm−1∥L2(Ω,H).

Using Assumption 2.1, the boundness of and ((Antonio2, , Lemma 4.2)), it holds that

 II2≤C∫tmtm−1∥F(X(s))∥L2(Ω,H)ds+C∫tmtm−1∥PhF(Xhm−1))∥L2(Ω,H)ds≤CΔt. (40)

Using the Itô-isometry property and triangle inequality, we split in two terms.

 II3 ≤ ∫tmtm−1∥∥∥(S(tm−s)−Sh(tm−s)Ph)Q12∥∥∥2L2(H)ds (41) + ∫tmtm−1∥∥∥(Sh(tm−s)−Sh(Δt))PhQ12∥∥∥2L2(H)ds := II31+II32.

Using Lemma 3.2 and Assumption 2.1, it holds that

 II31 ≤ ∫tmtm−1∞∑i=1∥(S(tm−s)−Sh(tm−s)Ph)Q12ψi∥2ds (42) = ∞∑i=1∫tmtm−1∥(S(tm−s)−Sh(tm−s)Ph)Q12ψi∥2ds ≤ ∞∑i=1Ch2β∥Q12ψi∥2β−1=Ch2β∥∥∥(−A)β−12Q12∥∥∥L2(H)≤Ch2β.

Using (33) it holds for small enough that

 II32 ≤ ∫tmtm−1∥∥∥Sh(tm−s)(I−Sh(s−tm−1))(−Ah)1−β2∥∥∥2L(H)∥∥∥(−Ah)β−12PhQ12∥∥∥2L2(H)ds (43) ≤ ≤ C∫tmtm−1(tm−s)−1+ϵ(s−tm−1)β−ϵds ≤ CΔtβ−ϵ∫tmtm−1(tm−s)−1+ϵds≤CΔtβ.

Substituting (43), (42), (40) and (39) in (38) yields

 ∥X(tm)−Xhm∥L2(Ω,H) ≤ Chβ+CΔtβ/2+C∥X(tm−1)−Xhm−1∥L2(Ω,H) (44) ≤ C(hβ+Δtβ/2)+C∥X(tm−1)−Xhm−1∥L2(Ω,H) + CΔtm−2∑k=0∥X(tk)−Xkk∥L2(Ω,H).

Applying the discrete Gronwall’s lemma (Lemma 3.3) to (44) yields

 ∥X(tm)−Xhm∥L2(Ω,H)≤C(hβ+Δtβ/2)exp(C+m−2∑k=0Δt)≤C(hβ+Δtβ/2). (45)

This completes the proof of Theorem 2.2.

## 4 Remark on current error estimates in the literature

In the literature, to obtain the error estimates many authors Xiaojie1 (); Xiaojie2 (); Antonio1 (); Antjd1 (); Jentzen1 (); Jentzen2 (); Raphael () usually iterate the numerical scheme, the mild solution of the problem (1) and the mild solution of the semi-discrete problem in space (24). Indeed we usually have

 Xhm=Sh(tm)PhX0+m−1∑k=0∫tk+1tkSh(tm−s)PhF(Xhk)ds+m−1∑k=0∫tk+1tkSh(Δt)PhdW(s) (46)

and

 Xh(tm)=Sh(tm)PhX0+m−1∑k=0∫tk+1tkS(tm−s)PhF(Xh(s))ds+m−1∑k=0∫tk+1tkSh(tm−s)PhdW(s). (47)

Subtracting (46) from (47), taking the norm and using triangle inequality yields

 ∥Xh(tm)−Xhm∥L2(Ω,H) ≤ m−1∑k=0∫tk+1tk∥∥Sh(tm−s)(PhF(Xh(s))−PhF(Xhk))∥∥L2(Ω,H)ds (48) + m−1∑k=0∥∥∥∫tk+1tk(Sh(tm−s)−Sh(Δt))dW(s)∥∥∥L2(Ω,H) =: III1+III2.

Note that compare to (38), the summation appears in the right hand side of (48). This usually leads to an optimal convergence order in time if the drift term is only Lipschitz continuous (see e.g. Antonio1 ()). In fact, the term involving the drift function (deterministic term) is usually estimated as follows.

 III1 ≤ m−1∑k=0∫tk+1tk∥S(tm−s)(PhF(Xh(s))−PhF(Xh(tk))∥L2(Ω,H) + m−1∑k=0∫tk+1tk∥S(tm−s)(PhF(Xh(tk))−PhF(Xhk))∥L2(Ω,H)ds ≤ Cm−1∑k=0∫tk+1tk∥Xh(s)−Xh(tk)∥L2(Ω,H)ds+Cm−1∑k=0∫tk+1tk∥Xh(tk)−Xhk∥L2(Ω,H)ds ≤ CΔtmin(β,1)/2+CΔtm−1∑k=0∥Xh(tk)−Xhk∥L2(Ω,H), (50)

where we have used the triangle inequality, (7) and (Antjd1, , Lemma 4).

To the best of our knowledge, to obtain convergence rate in time greater more than , almost all authors Xiaojie1 (); Xiaojie2 (); Antjd1 (); Jentzen1 (); Jentzen2 () require to be twice differentiable with derivatives satisfying certain assumptions (e.g. (Jentzen1, , Assumption 2.4) or (Xiaojie1, , Assumption 2.1)), so that they can apply the Taylor expansion to (4). But note that in the case where is not differentiable, this is not feasible. We fill that gap in this paper by a mean of a sharp integral estimate, a new regularity result (namely Lemma 2.1 and Theorem 2.1 and (20)). Note also that, in our convergence proof, we do not iterate the numerical scheme and the mild solution from time to .