Higher order stable schemes for stochastic convection–reaction–diffusion equations driven by additive Wiener noise

In this paper, we investigate the numerical approximation of stochastic convection–reaction–diffusion equations using two explicit exponential integrators. The stochastic partial differential equation (SPDE) is driven by additive Wiener process. The approximation in space is done via a combination of the standard finite element method and the Galerkin projection method. Using the linear functional of the noise, we construct two accelerated numerical methods, which achieve higher convergence orders. In particular, we achieve convergence rates approximately 1 for trace class noise and 12 for space‐time white noise. These convergence orders are obtained under less regularity assumptions on the nonlinear drift function than those used in the literature for stochastic reaction–diffusion equations. Numerical experiments to illustrate our theoretical results are provided.


INTRODUCTION
This article is devoted to the space-time approximation of the following semilinear parabolic stochastic PDEs (SPDEs) dX(t) = [AX(t) + (x, X(t))]dt + dW(t), t ∈ (0, T], x ∈ Λ, with initial value X(0) = X 0 and Dirichlet boundary conditions or Robin boundary conditions. In (1), Λ is a bounded domain of R d (d ∈ {1, 2, 3}) with smooth boundary or is a convex polygon. T > 0 is a fixed final time. The linear operator A is given by where q ij , q j ∈ L ∞ (Λ), and q i, j satisfy the following ellipticity condition: where c > 0 is a uniform constant. Precise conditions on the nonlinear f will be given in the next section. Let us introduce the Nemytskij operator F ∶ L 2 (Λ) → L 2 (Λ) associated to f in (1), defined by We denote by (Ω,  , P) a probability space with a filtration ( t ) t∈[0,T] ⊂  that fulfills the usual conditions; see, for example, Prévôt and Röckner. 1, Definition 2.1. 11 The noise term W(t) in (1) is assumed to be a Q-Wiener process defined on the filtered probability space (Ω,  , P, { t } t∈[0,T] ), with covariance operator Q ∶ L 2 (Λ) → L 2 (Λ), which is assumed to be linear, self-adjoint, and positive definite. It is well-known (see, e.g., Prévôt and Röckner 1 ) that the noise W(t) can be represented as follows: where (q i , e i ) i∈N d are the eigenvalues and eigenfunctions of the covariance operator Q and ( i ) i∈N d are independent and identically distributed standard Brownian motions. It is well-known that the linear operator A generates an analytic semi-group S(t) =∶ e At , t ≥ 0; see, for example, previous studies. [1][2][3][4] SPDEs of type (1) are used to model many real world phenomena such as convection-reaction-diffusion processes. Since explicit solutions of many SPDEs are usually unknown, numerical approaches are good alternatives to provide their realistic approximations. Having a numerical approximation in hand, one main question is whether it converges toward the mild solution or not. Another interesting information is to know the rate with which it converges to the true solution. There are mainly two types of convergence: namely, strong convergence and weak convergence. Our interest here is on strong convergence. There are numerous numerical methods designed to approximate (1) with linear self-adjoint operator; see, for example, previous studies. [5][6][7][8] To classify a numerical method, one also takes in consideration the rate of convergence. In Jentzen and Kloeden, 6 an exponential Euler scheme achieving higher convergence order by exploiting the linear functional of the noise was introduced for semilinear SPDEs driven by space-time white noise. However, assumptions made in Jentzen and Kloeden 6 to achieve higher convergence order are too restrictive and exclude many nonlinear operators such as F(v) = 1 − v∕1 + v 2 , v ∈  ∶= L 2 (Λ); see the introduction of Jentzen et al. 7 for more details. In Jentzen et al., 7 a modified version of the abovementioned exponential Euler scheme which achieves higher convergence order under more relaxed conditions on F was introduced. In Wang and Ruisheng, 9 an accelerated exponential integrator was investigated and proved to achieve convergence order 1 for trace class noise under more relaxed assumptions than in Jentzen et al. 7 Note that the works in other studies 6,7,9 are only for stochastic diffusion-reaction equations and heavily use the fact that the linear operator A is self-adjoint. Also, such schemes are numerically implementable only if the linear operator is self-adjoint. Here, we are interested in the case of stochastic convection-reaction-diffusion equations, which are more realistic and play a key role in subsurface processes. For such SPDEs, we are interested on building alternative stable numerical schemes achieving higher convergence order. Recently in Lord and Tambue, 10,11 some exponential integrators and implicit schemes achieving convergence order 1 in time were introduced for stochastic convection-reaction-diffusion equations. The idea behind such schemes consists of keeping the convection term in the nonlinear part F. As a consequence, these schemes may lose their stability if the velocity field is very high. Moreover, the convergence results in Lord and Tambue 10,11 exclude the space-time white noise case. Here, we propose novel numerical schemes for stochastic convection-diffusion-reaction equations with general noise (including space-time white noise), which achieve higher convergence order. The idea behind our novel numerical schemes consists of splitting the semi-group appearing in the noise component in two semi-groups, one semi-group generated by the self adjoint part, and another semi-group generated by the advection part (see Section 2 for details). For the convergence proofs of our numerical schemes toward the mild solution, one key argument consists of using an argument based on Miyadera-Voigt perturbation theorem. 2, Chapter III, Corollary 3. 16 In addition, we use relaxed conditions on the nonlinear function F than the ones used in Lord and Tambue 10,11 and in the literature for stochastic reaction diffusion equations. The rest of this paper is structured as follows. In Section 2, the well posedness and the numerical schemes are introduced. In Section 3, we prove the strong convergence of our fully discrete schemes toward the mild solution. In Section 4, we provide some numerical experiments to illustrate our results.

Notations and main assumptions
We denote by ⟨., .⟩ the inner product in the Hilbert space  = L 2 (Λ, R) =∶ L 2 (Λ). We denote by ||.|| the norm associated to the inner product ⟨., .⟩. For all p ≥ 2, L p (Ω, ) stands for the Banach space of all equivalence classes of p integrable -valued random variables. Let () be the space of bounded linear mappings on  endowed with the usual operator norm ||.|| () . The space of Hilbert-Schmidt operators from  to  is denoted by  2 () ∶= HS() and is equipped with the norm ||l|| 2 the following relation called Itô isometry holds: 1,4 For the seek of the convergence analysis of our numerical schemes, we make the following assumptions.
We assume the covariance operator Q to satisfy the following estimate: The nonlinear function F is assumed to be differentiable, and there exist C ≥ 0 and ∈ As a consequence, F satisfies the following Lipschitz condition: Remark 1. Let A 1 and A 2 be, respectively, the self-adjoint and the nonself-adjoint parts of A. Using the equivalence of norms (see, e.g., Fujita and Suzuki and Larsson 3,12 or and v ∈ ((−A) ), it follows that (7) remains true if A is replaced by A 1 . The following equivalence of norms holds:

Proposition 1. Under the hypothesis that the function f in (1) is differentiable and there exists a constant C ≥ 0 such that
the Nemystkii operator F satisfies the desired properties in Assumption 2.1. Note the the derivation in (8) is respect to that second variable.
Proof. We only prove that since proofs of other estimates are similar. The derivative of F is given by Hence, for w ∈ , from the definition of the norm ||.|| L 1 (Λ,R) , using (8), Hölder's inequality, and the fact that Λ is bounded, it follows that Using Hölder's inequality, the Sobolev embedding Using (9) and Tambue and Ngnotchouye, 15, Lemma 3.1 it follows from (10) that Moreover, the following hold:

Fully discrete schemes and main results
Let  h be a set of disjoint intervals of Λ (for d = 1), a triangulation of Λ (for d = 2), or a set of tetrahedra (for d = 3) satisfying the standard regularity assumptions (see Fujita and Suzuki 3 ). Let V h ⊂ V denotes the space of continuous functions that are piecewise linear over the triangulation  h . To discretize in space, we introduce two projections. Our first projection operator P h is the L 2 (Λ) projection onto V h defined for u ∈ L 2 (Λ) by where a(., .) is the corresponding bilinear form associated to the operator A. We denote by S h the semi-group generated by A h . The second projection P N , N ∈ N is the projection onto a finite number of spectral modes e i * defined for u ∈ L 2 (Λ) by The semi-discrete version of the problem (1) consists of finding The mild solution of (14) is given by Let S 1 (t) =∶ e A 1 t the semi-group generated by A 1 and S 2 (t) =∶ e A 2 t the semi-group generated by A 2 . Let A 1h and A 2h be the semi-discrete versions of A 1 and A 2 , respectively, defined as in (13). We denote by S 1h (t) and S 2h (t) the semi-groups generated A 1h and A 2h , respectively. The semi-groups S 1 (t), S 2 (t), S 1h (t), and S 2h (t) satisfy the smoothing properties of Mukam and Tambue 16, Proposition 2.2 ; see also Fujita and Suzuki and Larsson. 3, 12 We introduce the following stochastic convolutions: To build our numerical schemes, we use the following approximation of the noise: To build our first numerical scheme, we use the approximation F(X h (s)) ≈ F(X h (t m )), for s ∈ [t m , t m + 1 ). This yields the following scheme, called accelerated SETD1 (ASETD1): X h 0 = P h X 0 and recursively by Note that the numerical method (17) can be written in the following form, efficient for simulation: where the linear operator 1 is given by (21). To obtain our second numerical scheme, we use the approximation . This yields the following scheme, called accelerated SETD0 * Eigenfunctions of the operator A 1 in our case.
The numerical method (19) can be written in the following equivalent form, efficient for simulation: The linear operators 0 and 1 are given, respectively, by Note that ASETD1 and ASETD0 are the analog schemes in Lord and Tambue 10 that we are improving their stability in this paper.
The following theorem is the main result of this paper.
Theorem 2 (Main results). Let X(t m ) be the mild solution given by (11) where is given in Assumption 2.1 and ∈ (0, ) is any positive number, small enough.
Remark 2. Remember that as in Lord and Tambue, 10 to simulate our accelerated schemes, the eigenfunctions of the linear operator A should be the same as that of the covariance operator Q * , if not, the projection of the eigenfunctions of Q onto the eigenfunction of A should be done. This is indeed the drawback of the accelerated schemes as in general the projection is costly.

PROOF OF THE MAIN RESULTS
The proofs of the main results need some preliminaries results.

Preparatory results
Lemma 1. Let 0 ≤ ≤ 1 and u ∈ . Then the following sharp integral estimates hold: Proof. The proof of the first estimate can be found in Mukam Proof. First of all, as in Lord and Tambue, 19, Lemma 2.7 it holds that , where is any positive number. Similarly to Lord and Tambue, 19, Theorem 2.6 we have ||(−A) 2 X(t)|| L p (Ω,) < C for ∈ [0, 2). For = 2, we have Taking the norm in both sides of (23), using Assumption 2.1 and the stability properties of the semi-group yields Using the Burkhölder-Davis-Gundy inequality (Kruse 5, Lemma 5.1 ), Assumption 2.1, and Lemma 1, it follows that Substituting (25) in (24) and using (12), Assumption 2.1, and the estimate Using the second estimate of (22), one can readily prove that

Proof of Theorem 2
We only give the proof for ASETD1, since the case of ASETD0 is similar. Iterating the mild solution (11) yields Iterating the numerical scheme (17) yields Subtracting (27) from (26) yields Using Lemma 7 in Mukam and Tambue 16 with r = = , it holds that We can split II 2 in three terms as follows: We start with the estimates of II 22 and II 23 , since they are easier than that of II 21 . Using triangle inequality, Lemma 7 in Mukam and Tambue 16 with r = min( , 2 − ) and = 0, and Assumption 2.1, it holds that Using Assumption 2.1 and the smoothing properties of the semi-group yields Let us now estimate II 21 . Using Taylor's formula in Banach space yields From the mild solution, we have Substituting (34) in (33) yields where I k, s is given by Note that using Assumption 2.1, one easily checks that Substituting (35) in the expression of II 21 yields Using Lemma 2, (37), Assumption 2.1, and the stability of the semi-group yields In view of Assumption 2.1, one can easily get To estimate II (3) 21 , we recast it as follows: Let us start with the estimate of II (32) 21 . Using the triangle inequality and the Hölder inequality twice yields Using Burkhölder-Davis-Gundy inequality, 5 it holds Using Assumption 2.1 and the stability properties of the semi-group yields ≤ (s − r) min(0, −1) .

Using Assumption 2.1 and the stability properties of the semi-group yields
Using the definition of I h m,k,s , Assumption 2.1, and Lemma 2, we arrive at Substituting (45) in (44) and using Lemma 2 yield Substituting (46) in (42) yields Along the same lines as the estimate of II (32) 21 above, we obtain Substituting (47) and (48) Using triangle inequality, we split II 3 as follows: Note that the term II 31 can be written as follows: Let us recall that for any ∈ [0, 1], the following estimate holds; see, for example, Fujita H, Suzuki: 3 Using the Itô isometry, (54), Assumption 2.1, and Lemma 1 yields Using the Itô isometry, Assumption 2.1, the stability properties of the semi-group, and Mukam and Tambue 16, Lemma 7 with r = − and = max(0, − 1) yields (56) Substituting (56) and (55) in (53) yields Using the Itô isometry and the stability properties of the semi-group, it holds that ds. (58) Using the perturbation theorem of Miyadera-Voigt (Engel and Nagel 2, Theorem 3.14, Chapter III, (3.22) ), Note that one can easily check that conditions on applying Miyadera-Voigt theorem are fulfilled. This is due to Engel and Nagel. 2, Chapter III, Corollary 3. 16 Inserting Q 1 2 in (59) yields Taking the norm in both sides of (60) and using Assumption 2.1 and the stabilities properties of the semi-group, it follows that Substituting (61) in (58) yields Using the Itô isometry and splitting the sum in two parts yield We start with the estimate of ||II 33 || 2 L 2 (Ω,) . Using the stability properties of the semi-group yields ds.
Using the fact that A 1 , S 1 , and Q are self-adjoint and Assumption 2.1, it follows that Along the same lines as the estimate of ||II (2) 33 || 2 L 2 (Ω,) , we obtain Substituting (65) and (64) in (63) yields Using the Itô isometry and splitting the sum in two parts, it holds that It is well-known (see, e.g., Larsson 12 ) that ||A 2h v|| ≤ C ′ ||v|| 1 ≤ C||(− A 1 ) 1 2 v|| for all v ∈ V h . Hence, by interpolation theory, it holds that Using the latter estimate together with the smooth properties of the semi-group and Mukam and Tambue 16, Lemma 1 yields ds.

NUMERICAL EXPERIMENTS
As in Lord and Tambue, 10 we consider the following stochastic transport equation: where D > 0 is the diffusion coefficient, q is the Darcy's velocity field as in Lord and Tambue, 10 and D = 10 −2 . We consider two types of boundary conditions: For boundary condition (a), our function f used in (1) to define the Nemystkii operator F is given by One can easily check that (8) holds. In fact, simple estimates yields Therefore, from Proposition 1, it follows that the estimates regarding F in Assumption 2.1 are satisfied. For mixed boundary condition, the Nemystkii operator F also included the trace operator as we can observe in Lord and Tambue. 10 In this case, Assumption 2.1 is not satisfied as the domain of the trace operator is H 2 (Λ). Remember that as in Lord and Tambue, 10 to simulate our accelerated schemes, the eigenfunctions of the linear operator A should be the same as that of the covariance operator Q, if not the projection of the eigenfunctions of Q onto the eigenfunction of A should be done. As in Lord and Tambue, 10 our linear operator in all our simulations is the Laplace operator with Neumann boundary everywhere in the domain as the eigenvalues and eigenfunctions are well-known in rectangular domain. In the decomposition (5), we have used for some small > 0 and > 0. Since the eigenvalues of the Laplace operator with Neumann boundary are given by thus, Assumption 7 is satisfied. Details on simulation/implementation of the accelerated schemes can be found in Lord and Tambue. 10,11 In the legends of all of our graphs, we use the following notation: • ASETD1 is used for graphs from scheme (18) where the boundary condition (b) is used, while ASETD1 = N is used when Neumann boundary condition (a) is used. • ASETD0 is used for graphs from scheme (19) where the boundary condition (b) is used, while ASETD0 = N is used when Neumann boundary condition (a) is used. • SETD1 is used for graphs from the analog of ASETD1 scheme in Lord and Tambue 10 where the boundary condition (b) is used, while ASETD1 = N is used when Neumann boundary condition (a) is used. • ASETD0 is used for graphs from the analog of ASETD0 scheme in Lord and Tambue 10 where the boundary condition (b) is used, while ASETD0 = N is used when Neumann boundary condition (a) is used. • In Figure 3, ASETD1 = N1 and ASETD1 = N2 are used for graphs from scheme (18) where Neumann boundary condition (a) is used with the diffusion coefficient D = 10 −1 and D = 10 −2 , respectively.
We study the convergence for the both the small time steps and large time steps in order to show the weak stability of the schemes SETD1 and SETD0 presented in Lord and Tambue 10 and the good stability properties of our accelerated schemes ASETD1 and ASETD0 that we have proposed in this work. In all our graphs, the exact sample solutions are unknown, and the reference sample solutions in each scheme in our error computations are taken to be the numerical solution samples with that scheme with smallest time step size (1∕15360 for graphs with small time steps in Figure 1 and 1∕240 for graphs with large time steps in Figure 2). Figure 1 shows the convergence graphs with very small time steps with = 1 and = 2. The Peclet number which measures the rate of advection over the diffusion is 24. Although we have used very small time steps, we can observe that the well-known SETD1 and SETD0 schemes developed in Lord and Tambue 10 are unstable at the biggest time step in the graphs.
The convergence rate of all the schemes is close to 1 for = 1 (1.04 in Figure 1A and 1.01 in Figure 1B) and = 2 (1.08 in Figure 1C and 1.05 in Figure 1D). This is in agreement to our theoretical result in Theorem 2. Since for boundary condition (b), Assumption 2.1 is not satisfied for F, we can also conclude that Theorem 2 even holds for large class of nonlinear function F than what we have considered.
In some graphs, we can also observe that the schemes SETD0 and ASETD0 are less accurate comparing to the schemes SETD1 and ASETD1. Indeed, this is normal as the deterministic term on F is approximated accurately in SETD1 (ASETD1) than in ASETD0 (ASETD0) scheme. Figure 2 shows the convergence graphs with large time steps with = 1 and = 2. We can observe that the schemes ASETD1 and SETD0 are still stable for large time steps and that the convergence rate is still called to the theoretical result in Theorem 2 as we have about 0.95 in Figure 2 for = 1 and 1.05 for = 2. In Figure 3, we compare the errors and the orders of convergence for = 2 with two diffusion coefficients D = 10 −1 and D = 10 −2 . As we can observe, the errors for D = 10 −2 are smaller than the errors for D = 10 −1 . We have also observed that the order of convergence has reduced but still close to 1 as the orders of convergence are 0.93 for D = 10 −1 and 1.05 for D = 10 −2 . But we believe that using many samples or small time steps can yield orders of convergence close in the two cases. So, the diffusion coefficient may play an important role in the order of convergence; that is, in Theorem 2 may depend on the diffusion coefficient D.

FIGURE 3
Comparison of the errors and the orders of convergence for = 2 with two diffusion coefficients D = 10 −1 and D = 10 −2 . The errors for D = 10 −2 (ASETD1 = N2) are smaller than the errors for D = 10 −1 (ASETD1 = N1). The order of convergence has also reduced but still close to 1 as the orders of convergence are 0.93 for D = 10 −1 and 1.05 for D = 10 −2 [Colour figure can be viewed at wileyonlinelibrary.com]

CONCLUSION
In this paper, we have proposed two stable explicit exponential integrators to approximate semilinear parabolic partial differential equations driven by additive noise, with a linear operator not necessary self-adjoint. Such equations are also called stochastic convection-recation-diffusion equations. This generalizes the known results in the literature for reaction-diffusion equation. Moreover, our analysis is done under less regularity assumptions of the nonlinear drift function. Our schemes are accelerated and achieve higher convergence rate. For instance, for additive trace class noise, converge rate approximately 1 is recovered. We have provided numerical experiment to illustrate our findings.