Optimal strong convergence rates of some Euler-type timestepping schemes for the finite element discretization SPDEs driven by additive fractional Brownian motion and Poisson random measure

In this paper, we study the numerical approximation of a general second order semilinear stochastic partial differential equation (SPDE) driven by a additive fractional Brownian motion (fBm) with Hurst parameter $H>\frac 12$ and Poisson random measure, more realistic in modelling real world phenomena. To the best of our knowledge, numerical schemes for such SPDE have been lacked in scientific literature. The approximation is done with the standard finite element method in space and three Euler-type timestepping methods in time, more precisely linear implicit method, exponential integrator and exponential Rosenbrock scheme are used for time discretisation. In contract to the current literature in the field for SPDE driven only by fBm, our linear operator is not necessary self-adjoint and optimal strong convergence rates have been achieved for SPDE driven only by fBm and SPDE driven by fBm and Poisson measure. The results examine how the convergence orders depend on the regularity of the noise and the initial data and reveal that the full discretization attains an optimal convergence rate of order $\mathcal{O}(h^2+\Delta t)$ for the exponential integrator and implicit schemes (linear operator $A$ self-adjoint for implicit). Numerical experiments are provided to illustrate our theoretical results for the case of SPDE driven with fBm noise.

Optimal strong convergence rates of some Euler-type timestepping schemes for the finite element discretization SPDEs driven by additive fractional Brownian motion and Poisson random measure Aurelien Junior Noupelah, Antoine Tambue.
Received: date / Accepted: date Abstract In this paper, we study the numerical approximation of a general second order semilinear stochastic partial differential equation (SPDE) driven by a additive fractional Brownian motion (fBm) with Hurst parameter H > 1 2 and Poisson random measure, more realistic in modelling real world phenomena. To the best of our knowledge, numerical schemes for such SPDE have been lacked in scientific literature. The approximation is done with the standard finite element method in space and three Euler-type timestepping methods in time, more precisely linear implicit method, exponential integrator and exponential Rosenbrock scheme are used for time discretisation. In contract to the current literature in the field for SPDE driven only by fBm, our linear operator is not necessary self-adjoint and optimal strong convergence rates have been achieved for SPDE driven only by fBm and SPDE driven by fBm and Poisson measure. The results examine how the convergence orders depend on the regularity of the noise and the initial data and reveal that the full discretization attains an optimal convergence rate of order O(h 2 + ∆t) for the exponential integrator and implicit schemes (linear operator A self-adjoint for implicit).
Let N (dz, dt) be the H-valued Poisson distributed σ-finite measure on the product σ-algebra B(χ) and B(R + ) with intensity ν(dz)dt, where dt is the Lebesgue measure on B(R + ). In our model problem (1), N (dz, dt) stands for the compensated Poisson random measure defined by N (dz, dt) := N (dz, dt) − ν(dz)dt.
We denote by T > 0, the final time, F : H → H, φ are deterministic mappings that will be specified precisely later, X 0 is the initial data which is random, −A is a linear operator, not necessary self-adjoint, unbounded and generator of an analytic semigroup S(t) := e −tA , t ≥ 0. Note that B H (t) is a H-valued Q-cylindrical fractional Brownian motion of Hurst parameter H ∈ ( 1 2 , 1] in a filtered probability space (Ω, F, P, {F t } t≥0 ) with the covariance operator Q : H → H, which is positive definite and self-adjoint. The filtered probability space (Ω, F, P, {F t } t≥0 ) is assumed to fulfil the usual condition (see [25,Def 2.2.11]). It is well known [3] that the noise can be represented as where q i , e i , i ∈ N d are respectively the eigenvalues and eigenfunctions of the covariance operator Q, and β H i are mutually independent and identically distributed fractional Brownian motions (fBm).
In our study , we first study in details the following particular case where z 0 = 0, i.e, the SPDE is driven only by fBm dX(t) + AX(t)dt = F (X(t))dt + φ(t)dB H (t), The self-similar and long-range dependence properties of the fBm make this process a suitable candidate to model many phenomena like financial markets (see, e.g., [4,10,17]) and traffic networks (see, e.g., [5,35]). In most cases, SPDEs of type (5) do not have explicit solutions and therefore numerical algorithms are required for their approximations. It is important to mention that if H = 1 2 the process B H is not a semi-martingale and the standard stochastic calculus techniques are therefore obsolete while studying SPDEs of type (5). Alternative approaches to the standard Itô calculus are therefore required in order to build a stochastic calculus framework for such fBm. In recent years, there have been various developments of stochastic calculus and stochastic differential equations with respect to the fBm especially for H ∈ ( 1 2 , 1] (see, for example [2,3,19]) and theory of SPDEs driven by fractional Brownian motion has been also studied. For example, linear and semilinear stochastic equations in a Hilbert space with an infinite dimensional fractional Brownian motion are considered in [6,7]. In contrast to standard Brownian (H = 1/2) where there are numerous literature on numerical algorithms for SPDEs, few works have been done for numerical methods for fBm for SPDEs of type (5). Indeed, standard explicit and linear implicit schemes have been investigated in the literature for SPDEs of type (5) (see [12,13,34]). The works in [12,34] deal with self-adjoint operator and use the spectral Galerkin method for the spatial discretization. This is very restrictive as many concrete applications use non self-adjoint operators. Beside numerical algorithms used for spatial discretization and time discretization in [12,34] are limited to few applications. Our goal in this work is to extend keys time stepping methods, which have been built for standard Brownian motion (H = 1/2). These extensions are extremely complicated due to the fact that the process B H is not a semi-martingale. Our results will be based on many novel intermediate lemmas. Indeed our schemes here are based on finite element method (or finite volume method) for spatial discretization so that we gain the flexibility of these methods to deal with complex boundary conditions and we can apply well-developed techniques such as upwinding to deal with advection. For time discretization, we first updated implicit linear for finite element method and not necessarily self-adjoint. We also provide the strong convergence of the exponential scheme [16] for (H ∈ ( 1 2 , 1]). Note that this scheme is an explicit stable scheme, where the implementation is based on the computation of matrix exponential functions [16]. As the linear implicit and exponential scheme are stable only when the linear operator A is stronger than the nonlinear function F 1 , we also provide the strong con-vergence of the Stochastic Exponential Rosenbrock Scheme (SERS) [20] for (H ∈ ( 1 2 , 1]), which is very stable when (5) is driven both by its linear or nonlinear part.
However the model equation (5) can be unsatisfactory and less realistic. For instance, in finance, the unpredictable nature of many events such as market crashes, announcements made by the central banks, changing credit risk, insurance in a changing risk, changing face of operational risk [24,27] might have sudden and significant impacts on the stock price. As for standard Brownian motion, we can incorporate a non-Gaussian noise such as Lévy process or Poisson random measure to model such events. The corresponding equation is our model equation given in (1). In contrast to SPDE driven by fBm in (5) where at least few numerical schemes exist, numerical schemes for such SPDE of type (1) driven by fBm and Poisson measure have been lacked in scientific literature, to the best of our knowledge. In this work, we will also fill the gap by extending the implicit scheme, the exponential scheme and the Stochastic Exponential Rosenbrock Scheme to SPDE of type (1). For SPDE of type (5) and SPDE of type (1), our strong convergence results examine how the convergence orders depend on the regularity of the noise and the initial data and reveal that the full discretization attains an optimal convergence rate of order O(h 2 + ∆t) for the exponential integrator and implicit schemes (linear operator A self-adjoint for implicit).
The paper is structured as follows. In Section 2, Mathematical setting for fBm is presented, along with the well posedness and regularities results of the mild solution of SPDE (5) driven by fBm. In Section 3, numerical schemes based on implicit scheme, stochastic exponential integrator and stochastic exponential Rosenbrock scheme for SPDE (5) driven by fBm are presented. In Section 4, the strong convergence proofs of schemes presented in Section 3 are provided. In Section 5, numerical schemes based on semi implicit scheme, stochastic exponential integrator scheme and stochastic exponential Rosenbrock scheme are presented for SPDE (1) driven by fBm and Poisson measure, along with the extension of their strong convergence proofs. We end the paper in Section 6 with numerical experiments illustrating our theoretical results for SPDE (5) driven by fBm noise.

Mathematical setting
In this section, we review some standard results on fractional calculus and introduce notations, definitions and preliminaries results which will be needed throughout this paper.
Throughout this paper the Hurst parameter H is assumed to be in the interval (1/2, 1]. Let (K, ., . K , . ) be a separable Hilbert space. For p ≥ 2 and for a Banach space U, we denote by L p (Ω, U ) the Banach space of p-integrable U -valued random variables. We denote by L(U, K) the space of bounded linear mapping from U to K endowed with the usual operator norm . L(U,K) and L 2 (U, K) = HS(U, K) the space of Hilbert-Schmidt operators from U to K with In that follows, we will make some assumptions on F , φ, X 0 and A, which will allow us to ensure the existence and uniqueness of the mild solution X of (5) represented by (see e.g [34]) for t ∈ [0, T ]. To ensure the existence and the uniqueness of solution for SPDE (5) and for the purpose of convergence analysis, we make the following assumptions.
Assumption 1 (Noise term) We assume that for some constant β ∈ (0, 1] Assumption 2 (Non linearity) For the deterministic mapping F : H → H, we assume that there exists constant L ∈ (0, ∞) such that As a consequence of (17) it holds that In the Banach space D A α 2 , α ∈ R, we use notation A α 2 · = · α and we recall the following properties of the semigroup S(t) generated by −A, that will be useful throughout this paper.
Proposition 1 (Smoothing properties of the semigroup) [23] Let α > 0, δ ≥ 0 and 0 ≤ γ ≤ 1, then there exist a constant C > 0 such that where l = 0, 1 and The next lemma (specially (23) and (24)) is an important result which plays a crucial role to obtain regularity results, very useful in this work.
Lemma 2 For any 0 ≤ ρ ≤ 1, 0 ≤ γ ≤ 2 and 0 ≤ υ ≤ H with H ∈ 1 2 , 1 , if the linear operator is given by (33), there exists a positive constant C such that for all Proof. See [21, Lemma 2.1] for the proof of (21) and (22). Concerning the proof of (23), the border case H = 1 2 if obtained using (21) with ρ = 1 and the order border case H = 1 is also obtained using (22) with γ = 2. Hence the proof of (23) is thus completed by interpolation theory. The proof of (24) for 0 ≤ υ ≤ H is an immediate consequence of Proposition 1. The border case υ = H is proved by (23). This completes the proof of Lemma 2.
Remark 2 Proposition 1 and Lemma 2 also hold with a uniform constant C (independent of h) when A and S(t) are replaced respectively by their discrete versions A h and S h (t) defined in Section 3, see e.g. [15,16].
The well posedness result is given in the following theorem along with optimal regularity results in both space and time.
Theorem 1 Assume that Assumptions 1-3 are satisfied, then there exists a unique mild solution given by (14) such that for all t ∈ [0, T ], X(t) ∈ Moreover, if the linear operator is given by (33), the following optimal regularity results in space and time hold (27) and for 0 ≤ t 1 < t 2 ≤ T ; Where C = C(β, L, T, H) is a positive constant and β is the regularity parameter of Assumption 1.
Proof [34,Theorem 3.5] gives the result of existence and uniqueness of the mild solution X. For regularity in space, we adapt from [20, Theorem 2.1 (23), (24)] by just replacing β in their case by 2H + β − 1. The difference will therefore be made at the level of the estimate of the stochastic integral To reach our goal, we use triangle inequality, the estimate (a + b) 2 ≤ 2a 2 + 2b 2 , (12) and (13), Assumption 1, Proposition 1, Lemma 2 (23) to have For the proof of (28), triangle inequality yields .

Numerical schemes
Throughout this section, we assume that Λ is bounded and has smooth boundary or is a convex polygon of R d , d ∈ {1, 2, 3}. In the rest of this paper we consider the SPDE (5) to be of the following form where the function f : Λ × R −→ R is continuously twice differentiable and the function b : Λ×R −→ R is globally Lipschitz with respect to the second variable. In the abstract framework (5), the linear operator A takes the form The functions F : Section 4] that the Nemytskii operator F related to f and the operator φ associated to b defined in (35) satisfy Assumption 1 and Assumption 2. As in [8,16] we introduce two spaces H and V , such that H ⊂ V ; the two spaces depend on the boundary conditions and the domain of the operator A. For Dirichlet (or first-type) boundary conditions we take For Robin (third-type) boundary condition and Neumann (second-type) boundary condition, which is a special case of Robin boundary condition, we take where ∂v/∂v A is the normal derivative of v and v A is the exterior pointing normal n = (n i ) to the boundary of A, given by Using the Green's formula and the boundary conditions, the corresponding bilinear form associated to A is given by for Dirichlet and Neumann boundary conditions, and for Robin boundary conditions. Using the Gårding's inequality, it holds that there exist two constants c 0 and λ 0 such that By adding and substracting c 0 Xdt in both sides of (5), we have a new linear operator still denoted by A, and the corresponding bilinear form is also still denoted by a. Therefore, the following coercivity property holds Note that the expression of the nonlinear term F has changed as we included the term c 0 X in a new nonlinear term that we still denote by F . The coercivity where S θ = λ ∈ C : λ = ρe iφ , ρ > 0, 0 ≤ |φ| ≤ θ (see [9]). Then −A is the infinitesimal generator of a bounded analytic semigroup S(t) = e −tA on L 2 (Λ) such that where C denotes a path that surrounds the spectrum of −A. The coercivity property (37) also implies that A is a positive operator and its fractional powers are well defined for any α > 0, by where Γ (α) is the Gamma function (see [9]). Under condition (34), it is well known (see e.g. [8]) that the linear operator −A given by (33) generates an analytic semigroup S(t) ≡ e −tA . Following [8,16], we characterize the domain of the operator A r/2 denoted by D(A r/2 ), r ∈ {1, 2} with the following equivalence of norms, useful in our convergence proofs We consider the discretization of the spatial domain by a finite element triangulation [30,33]. Let T 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) with maximal length h satisfying the usual regularity assumptions. Let V h ⊂ H denote the space of continuous functions that are piecewise linear over triangulation T h . To discretise in space, we introduce the projection P h from L 2 (Ω) to V h define for u ∈ L 2 (Ω) by The discrete operator where a is the corresponding bilinear form of A. Like the operator A, the discrete operator −A h is also the generator of an The mild solution of (43) can be represented as follows and we have the following regularity results.
Lemma 3 Assume that Assumptions 1-3 are satisfied, then the unique mild solution X h (t) given by (44) satisfied and for 0 ≤ t 1 < t 2 ≤ T ; Proof Since the operators A h and S h (t) satisfy the same properties as A and S(t) ( see Remark 2), then by using [30, (83)] and the boundedness of P h in the proof of (27) and (28), we obtain the proof of the expression (45) and (46). The proof of Lemma 3 is thus completed. Now applying the linear implicit Euler method [12,33] to (43) gives the following fully discrete scheme Furthermore applying the stochastic exponential integrator ( [16], SETD1) and Rosenbrock scheme ( [20], SERS) to (43) yields and Where and Note that the exponential integrator scheme (48) is an explicit stable scheme when the SPDE (5) is driven by its linear part as the linear implicit method, while the Stochastic Exponential Rosenbrock Scheme (SERS) (49) is very stable when (5) is driven by its linear or nonlinear part. When dealing with SERS, the strong convergence proof will make use of the following assumption, also used in [20,21].
Assumption 4 For the deterministic mapping F : H → H, we also assume that there exists constant L ∈ (0, ∞) such that

Main Result for SPDE driven by fBm
Theorem 2 Let X(t m ) be the mild solution of (5) at time t m = m∆t, ∆t ≥ 0 represented by (14). Let ζ h m be the numerical approximations through (47) and (49)(ζ h m = X h m for implicit scheme, ζ h m = Z h m for SERS). Under Assumptions 1-3 and 4 (essentially for SERS), β ∈ (0, 1], then the following holds and where is a positive constant small enough.

Proofs of the main result for SPDE with fBm
We introduce the Riesz representation under the regularity assumptions on the triangulation and in view of the Vellipticity, it is well known ( [8,15]) that the following error bound holds: for r ∈ [1,2]. Let us consider the following deterministic linear problem: The corresponding semi-discrete problem in space consists to finding u h ∈ V h such that Let us define the following operator Then we have the following lemma The following estimates hold for the semi-discrete approximation of (43). There exists a constant C > 0 such that Proof (see [30, Lemma 3.1 and Lemma 3.2 (iv) and (v)]) for the proof of (i) − (iii). Let us prove (iv).
-For H = 1 2 , using (60) with γ = δ, we obtain (64) Hence the proof of (62) is thus completed by interpolation theory. Lemma 5 (Space error) Let Assumptions 1-3 be fulfilled, then the following error estimate holds for the mild solution (14) and the discrete problem (44) holds Proof Using triangle inequality, we have .

Proof of Theorem 2 for implicit scheme
It is important to mention that the estimates made in this section are inspired by the results in [33, (4.7)-(4.14), (4.25)-(4.29)], when the linear operatorA is self-adjoint. For our case where A is not necessarily self-adjoint, let us present some preparatory results.
Lemma 6 For any m, h and ∆t the following estimates holds (i) for any j = 1, 2, · · ·, M .
Where is an arbitrary small number.
Proof See [30,Lemma 3.5] for the proof of (i). For the proof of (ii), we have m j=1 tj tj−1 Using triangle inequality and the estimate (a + b) k ≤ 2 k−1 a k + 2 k−1 b k (with k ≥ 1 and a, b ≥ 0) we obtain By Lemma 6 (ii) with γ = µ we obtain By inserting an appropriate power of A h , [30, (81)] and Remark 2 yields Substituting (81) and (82) in (80) yields Concerning the estimate of K 2 , let > 0 small enough, Lemma 6 (iii) with γ = µ yields Adding (83) and (84) This completes the proof of Lemma 7. With these two lemmas, we are now ready to prove our theorem for the implicit scheme . In fact, using the standard technique in the error analysis, we split the fully discrete error in two terms as Note that the space error err 0 is estimated by Lemma 5. It remains to estimate the time error err 1 . We recall that the exact solution at t m of the semidiscrete problem (43) is given by We also recall that the numerical solution at t m given by (47) can be rewritten as  (87) and (88) that As we said at the beginning of this section, following closely the work done in [33, (4.7)-(4.14)] and replacing its preparatory results with Lemma 6 (i), (iv) with µ = 2H + β − 1, (v), Lemma 7 (i) with ρ = 0, Remark 2 (19) with δ = γ = 1, Assumptions 2-3, boundedness of S h (t m − s) and P h , the stability properties of a discrete semigroup S h (t), (17) and (46), we have Note that in this work, we do not need to impose an assumption on F to increase the convergence rate as it is done in [33]. Indeed, thanks to (46) the following estimate is largely sufficient to reach a higher rate.
In what follows, we will present a corollary of Theorem 2 for the implicit Euler scheme where the linear operator A is assumed to be self-adjoint. The optimal strong convergence rate in time O(∆t) is reached.
Corollary 1 Let X(t m ) be the mild solution of (5) (A self-adjoint) at time t m = m∆t, ∆t ≥ 0 represented by (14). Let X h m be the numerical approximation through (47). Under Assumptions 1-3, β ∈ (0, 1], then the following holds For the proof of this corollary, we need to update our preparatory results, more precisely Lemma 7 in the self-adjoint case. The result is presented in the following lemma: and where (103) Proof. See [14, Proof of Lemma 4.4 (i)] for the proof of (i) and [34,Lemmas 4.8 and 4.9], [30, (83)] for the proof of (ii).
With this new lemma, we are now in position to prove our Corollary 1.

Proof of Theorem 2 for SETD1
As usual, spliting the fully discrete error in two terms yields Since the space error err 0 has been estimated by Lemma 5, we only need to estimate the time error err 2 . Remember that the exact solution at t m is given by and we recall that the numerical solution at t m given by (48) can be rewritten as where the notations [t] and [t] m are given by (89). By (109) and (110), we have Applying the triangle inequality yields Using the boundedness of P h and S h (t m − s), Lemma 2 and (46), we easily have and Adding (113) and (114), we obtain We estimate at now I 2 . Using triangle inequality and the estimate (a + b) 2 ≤ 2a 2 + 2b 2 , we split it in three terms Thanks to (94) we have Thereafter, (13), the change of variable j = m − k and ς = t m − s, inserting an appropriate power of A h , [30, (81)] and Remark 2 (more precisely (19) with 20) and Assumption 1 (more precisely (16)) yields Finally, using (12), inserting an appropriate power of A h , [30, (81)] and Remark 2 (more precisely (19) with γ = 2H+β−1 2 , (20) and (23)), Assumption 1 ( more precisely (15)) we obtain hence inserting (117)-(119) in (116) and taking the square-root gives Adding (115) and (120) we obtain Using the discrete version of the Gronwall inequality yields Combining (65) and (122) completes the proof.

Proof of Theorem 2 for SERS scheme
Before moving to the proof, we first present some preparatory results. Thanks to Assumption 4 and the works done in [20] we obtain Lemma 9 [20, Lemma 5] For all M ∈ N and all ω ∈ Ω, there is a positive constant C 1 independent of h, m, ∆t and the sample ω such that Lemma 10 [20, Lemma 6] The function G h m (ω) defined by (50) satisfies the global Lipschitz condition with a uniform constant C > 0, independent of h, m and ω such that Lemma 11 [20,Lemma 9] For all ω ∈ Ω, the stochastic perturbed semigroup S h m (ω)(t) := e (−A h +J h m (w))t satisfies the following properties (iv) For γ 1 , γ 2 > 0 such that 0 ≤ γ 1 − γ 2 ≤ 1, then the following estimate holds Lemma 12 [20, Lemma 10] The stochastic perturbed semigroup S h m (ω) satisfies the following property Where C is a positive constant independent of m, k, h, ∆t and the sample ω.
We can now prove our theorem. As in the proof of the previous schemes, we split the fully discrete error in two terms as By Lemma 5 we have the estimate of the space error Consider now the estimate of the time error err 3 . We recall that the semidiscrete problem (43) can be rewritten as for all t m ≤ t ≤ t m+1 where J h m and G h m is given by (49) and (50). Hence the exact solution at time t m of the semidiscrete problem (124) is given by and the numerical solution (49) can be rewritten as If m = 1 then from (125) and (126) we obtain By [20, (93)] we have the estimate Using triangle inequality and the estimate (a + b) 2 ≤ 2a 2 + 2b 2 we split it in three terms For the estimate II 2 1 , using (13), inserting an appropriate power of A h , [30, (81)], Assumption 1 (more precisely (16)) and Lemma 11 (ii) with γ 1 = 1−β 2 we obtain We denote by a positive constant small enough, using (12), inserting an appropriate power of A h , Assumption 1, [30, (83)], Lemma 11 (ii) with Hence putting (131) and (131) in (129) and taking the square-root gives Adding (128) and (132) yields For m ≥ 2, we recall that the solution at t m of the semidiscrete problem (124) is given by We recall also that the numerical solution at t m given by (126) can be rewritten as Using (134), (135) and the triangle inequality, we have , .
5 Extension to SPDE driven simultaneously by fBm and Poisson random measure

Numerical schemes
Here the goal is to show how the previous results can be extended to the following SPDE driven simultaneously by fBm and Poisson random measure. The corresponding model equation is given by In the Hilbert space H = L 2 (Λ), (149) is equivalent to (1) where the linear operator A and the nonlinear function F are defined as in (33) and (35). The well posedness result for H = 1/2 presented in [1] can easily be extended to (1) for H ∈ [1/2, 1] by combining with [34]. The corresponding exponential Euler (SETD1) scheme in integral form is therefore given with Y h 0 = P h X 0 . In the same way, the semi implicit scheme is given by while the stochastic Exponential Rosenbrock Scheme (SERS) is given by To obtain the optimal order in time, as in [22], we need the following assumption in Poisson measure noise.
Assumption 5 The covariance operator Q : H −→ H and the jump coefficient satisfy the following estimate where η = 2H + β − 1 with β ∈ (0, 1] as in Assumption 3 and Assumption 1 . Remark 1 All the regularity results in space and time, both for continuous equation (1) (or semi-discrete equation) important to achieve optimal convergence orders can easily be extended from our results on fBm in Theorem 1 and Lemma 3 by just following [22, Proposition 3.1].

Convergence results for SPDE with fBm and Poisson measure noise
The convergence result is exactly as for fBm when Assumption 5 is used Theorem 3 Let X(t m ) be the mild solution of (1) at time t m = m∆t, ∆t ≥ 0. Let ζ h m be the numerical approximations through (151) and (126)(ζ h m = X h

Numerical simulations
In opposite to the standard Brownian motion where the simulation is obvious, the simulation of fBm is not obvious and is an important research field in numerical analysis. Keys methods for simulations of fBm are Cholesky method [12] and the circulant method [26], which will be used in this work to generate the fBm. Here we consider the stochastic advection diffusion reaction SPDE where l ∈ {1, 2} , x ∈ Λ. In the noise representation (4), we have used for some small δ > 0. We have used b(x, t) = 2 in (32), so φ in Assumption 1 is obviously satisfied for β = (0, 1]. In our simulations, we have used δ = 0.001. The function f used in (35) to be f (x, z) = z 1+z for all (x, z) ∈ Λ × R. Therefore the corresponding Nemytskii operator F defined by (35) obviously satisfies Assumption 2. We obtain the Darcy velocity field q = (q i ) by solving the following system ∇ · q = 0, q = −k∇p, and −k ∇p(x, t) · n = 0 in Γ 1 N . Note that k is the permeability tensor and p the presure. We use a random permeability field as in [28, Figure 6]. The streamline of the velocity field q are given in Figure 1(d). To deal with high Péclet number, we discretise in space using finite volume method, viewed as a finite element method (see [29]). We take L 1 = 3 and L 2 = 2 and our reference solutions samples are numerical solutions using at time step of ∆t = 1/4096. The errors are computed at the final time T = 1. The initial solution is X 0 = 0, so we can therefore expect high orders convergence, which depend only on the noise term and H.   Figure 1(c) is the errors graph for the exponential Rosenbrock scheme with two values of H. We have observed the order of convergence is 0.5562 in time for H = 0.51 and β = 1, 0.6197 for H = 0.65 and β = 1.
As we can observe, our numerical orders in time are close to our theoretical results in Theorem 2 even if we have only used 50 samples in our Monte Carlo simulations. Figure 2 shows two samples of the solution for H = 0.75 and H = 0.51. Here we have fixed β = 1 and same Gaussian randoms numbers have been to generate our fBm samples. As we can observe, the parameter H has significant influence on the sample of the numerical solution. This is independent of our timestepping methods.