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>\frac {1}{2}\) and Poisson random measure. Such equations are more realistic in modelling real world phenomena. To the best of our knowledge, numerical schemes for such SPDE have been lacked in the scientific literature. The approximation is done with the standard finite element method in space and three Euler-type timestepping methods in time. More precisely the well-known linear implicit method, an exponential integrator and the exponential Rosenbrock scheme are used for time discretization. In contract to the current literature in the field, our linear operator is not necessary self-adjoint and we have achieved optimal strong convergence rates for 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 the optimal convergence rates of order \(\mathcal {O}(h^{2}+\varDelta t)\) for the exponential integrator and implicit schemes. Numerical experiments are provided to illustrate our theoretical results for the case of SPDE driven by the fBm noise.
Similar content being viewed by others
1 Introduction
We analyse the strong numerical approximation of an SPDE defines in \(\varLambda \subset \mathbb {R}^{d}\), d ∈{1,2,3} with initial value and boundary conditions (Dirichlet, Neumann, Robin boundary conditions or mixed Dirichlet and Neumann). In Hilbert space, our model equation can be formulated as the following parabolic SPDE
in Hilbert space \({\mathscr{H}}=L^{2}(\varLambda )\), with z0 ∈ χ, where χ is the mark set defined by \(\chi :={\mathscr{H}}\setminus \{0\}\). Let \({\mathscr{B}}(\varGamma )\) be the smallest σ-algebra containing all open sets of Γ. Let \((\chi , {\mathscr{B}}(\chi ), \nu )\) be a σ-finite measurable space and ν (with ν≢0) a Lévy measure on \({\mathscr{B}}(\chi )\) such that
Let N(dz,dt) be the \({\mathscr{H}}\)-valued Poisson distributed σ-finite measure on the product σ-algebra \({\mathscr{B}}(\chi )\) and \({\mathscr{B}}(\mathbb {R}_{+})\) with intensity ν(dz)dt, where dt is the Lebesgue measure on \({\mathscr{B}}(\mathbb {R}_{+})\). In our model problem (1), \(\widetilde {N}(dz, dt)\) stands for the compensated Poisson random measure defined by
We denote by T > 0, the final time, \(F:{\mathscr{H}}\rightarrow {\mathscr{H}}\), ϕ are deterministic mappings that will be specified precisely later, X0 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 BH(t) is a \({\mathscr{H}}\)-valued Q-cylindrical fractional Brownian motion of Hurst parameter \(H\in (\frac {1}{2},1]\) in a filtered probability space \((\varOmega , \mathcal {F}, \mathbb {P},\{\mathcal {F}_{t}\}_{t\geq 0})\) with the covariance operator \(Q: {\mathscr{H}}\rightarrow {\mathscr{H}}\), which is positive definite and self-adjoint. The filtered probability space \((\varOmega , \mathcal {F}, \mathbb {P},\{\mathcal {F}_{t}\}_{t\geq 0})\) is assumed to fulfil the usual condition (see [31, Definition 2.2.11]). It is well known [3] that the noise can be represented as
where \(q_{i}, e_{i}, i\in \mathbb {N}^{d}\) are respectively the eigenvalues and eigenfunctions of the covariance operator Q, and \({\beta _{i}^{H}}\) are mutually independent and identically distributed fractional Brownian motions (fBm).
In our study, we first study in details the following particular case where z0 = 0, i.e. the SPDE is driven only by fBm
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, 11, 22]) and traffic networks (see, e.g., [17, 39]). 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\neq \frac {1}{2}\) the process BH 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\in (\frac {1}{2}, 1]\) (see, for example, [2, 3, 24, 25]) and theory of SPDEs driven by fractional Brownian motion has been also studied. For example, stochastic partial differential equations in a Hilbert space with an infinite dimensional fractional Brownian motion are considered in [3, 6, 7, 18,19,20]. 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 [13, 14, 38]). The works in [13, 38] 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 [13, 38] 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 BH 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 will first update the implicit linear for finite element method and not necessarily self-adjoint. We also provide the strong convergence of the exponential scheme [21, 34] for (\(H\in (\frac {1}{2}, 1]\)). Note that this scheme is an explicit stable scheme, where the implementation is based on the computation of matrix exponential functions [21]. As the linear implicit and exponential scheme are stable only when the linear operator A is stronger than the nonlinear function F,Footnote 1 we also provide the strong convergence of the stochastic exponential Rosenbrock scheme (SERS) [26] for (\(H \in (\frac {1}{2}, 1]\)), which is very stable when (5) is driven both by its linear or nonlinear parts.
However, the model (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 [5, 30] 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 optimal convergence rates of order \(\mathcal {O}(h^{2}+\varDelta t)\) for the exponential integrator and implicit schemes.Footnote 2
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.
2 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.
Definition 1
[13, 23, 25, 38] The fractional Brownian motion (fBm) of Hurst parameter H ∈ (0,1) is a centered Gaussian process βH = {βH(t),t ≥ 0} with the covariance function
Notice that if \(H=\frac {1}{2}\), the process is the standard Brownian motion.
Remark 1
[25, Remark 1.2.3] For H = 1, we set βH(t) = β1(t) = tξ, where ξ is a standard normal random variable.
Throughout this paper the Hurst parameter H is assumed to be in the interval (1/2,1]. Let \(\left (K,\langle .,.\rangle _{K},\|.\|\right )\) be a separable Hilbert space. For p ≥ 2 and for a Banach space U, we denote by Lp(Ω,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 \({\mathscr{L}}_{2}(U,K)=HS(U,K)\) the space of Hilbert-Schmidt operators from U to K with
where \((e_{i})_{i\in \mathbb {N}^{d}}\) is an orthonormal basis on U. The sum in (6) is independent of the choice of the orthonormal basis in U. For simplicity, we use the notation L(U,U) =: L(U) and \({\mathscr{L}}_{2}(U,U)=:{\mathscr{L}}_{2}(U)\). It is well known that for all l ∈ L(U,K) and \(l_{1}\in {\mathscr{L}}_{2}(U)\), \(ll_{1}\in {\mathscr{L}}_{2}(U,K)\) and
We denote by \({L^{0}_{2}}\), the space of Hilbert-Schmidt operators from \(Q^{\frac {1}{2}}({\mathscr{H}})\) to \({\mathscr{H}}\) by \({L^{0}_{2}}:=HS(Q^{\frac {1}{2}}({\mathscr{H}}),{\mathscr{H}})\) with corresponding norm \(\|.\|_{{L^{0}_{2}}}\) by
Now, let us introduce the Wiener integral with respect to the one-dimensional fBm βH. Let T > 0 and the linear space \(\mathcal {E}\) denotes the collection of all \(\mathbb {R}\)-valued step functions defined on [0,T], that is, \(\phi \in \mathcal {E}\), if
where t ∈ [0,T], \(\phi _{i}\in \mathbb {R}\) and 0 = t1 < t2 < ⋅⋅⋅ < tn = T. For \(\phi \in \mathcal {E}\), we define its Wiener integral with respect to βH as
Let \(\mathcal {Z}\) be the Hilbert space defined as the closure of \(\mathcal {E}\) with respect to the scalar product
Then, the mapping
is an isometry between \(\mathcal {E}\) and the linear space span{βH,t ∈ [0,T]}, which can be extended to an isometry between \(\mathcal {Z}\) and the first Wiener chaos of the fBm \(\overline {\text {span}}^{L^{2}(\varOmega )}\{\beta _{H}, t\in [0,T]\}\) (see [3, 18]). The image of an element \(\varphi \in \mathcal {Z}\) by this isometry is called the Wiener integral of φ with respect to βH. Next introduce the square integrable Kernel
and its derivative
where \(c_{H} = \left (\frac {H(2H-1)}{B(2H-2,H-1/2)}\right )^{1/2}\), with B denoting the Beta function and t > s. Let us define the operator \(K_{H}^{*}\) from \(\mathcal {E}\) to L2([0,T]) by
We easily have
and
Hence, the operator \(K_{H}^{*}\) is an isometry between \(\mathcal {E}\) and L2([0,T]) which can be extended to \(\mathcal {Z}\). By definition
is a Brownian sheet, and in turn the fractional noise has a representation
In addition, for any \(\varphi \in \mathcal {Z}\)
if and only if \(K_{H}^{*}\varphi \in L^{2}([0,T])\).
Also denoting \(L^{2}_{\mathcal {Z}}([0,T]) = \{\varphi \in \mathcal {Z}, K_{H}^{*}\varphi \in L^{2}([0,T])\}\), since H > 1/2 we have the following embedding property.
Proposition 1
[7, 25] \(L^{2}([0,T])\subset L^{1/H}([0,T])\subset L^{2}_{\mathcal {Z}}([0,T])\).
The next proposition shows the first and second moments of stochastic integral defined for function in \(L^{2}_{\mathcal {Z}}([0,T])\).
Proposition 2
[6, 7] For \(f\in L^{2}_{\mathcal {Z}}([0,T])\) we have
and
Moreover, we have the following lemma very important throughout this paper.
Lemma 1
[3, Lemma 1] For any φ ∈ L1/H([0,T]), the following inequality holds
where CH > 0 is a constant depending only on H.
Next, we are interested in considering an fBm with values in a Hilbert space and giving the definition of the corresponding stochastic integral.
Definition 2
Let \(f\in L^{2}\left ([0,T],{L^{0}_{2}}\right )\). Then, its stochastic integral with respect to the fBm BH is defined, for t ≥ 0, as follows
As the sequence of random variables \(\left ({{\int \limits }_{0}^{T}}f(s)Q^{\frac {1}{2}}e_{i}d{\beta _{i}^{H}}(s), i\in \mathbb {N}^{d}\right )\) are mutually independent Gaussian random variable, by Propositions 1, 2 and Lemma 1, the mean of random variable (12) is also zero and we have the following result for its second moment.
Lemma 2
[3, 18] Let \(f\in L^{2}([0,T];{L^{0}_{2}})\) and \(H\in (\frac {1}{2},1]\), then the following holds
In that follows, we will make some assumptions on F, ϕ, X0 and A, which will allow us to ensure the existence and uniqueness of the mild solution X of (5) represented by (see, e.g., [38])
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] and \(\delta \in \left [\frac {2H+\beta -1}{2},1\right ]\), the deterministic mapping \(\phi :[0,T]\times \varLambda \rightarrow {L_{2}^{0}}\) satisfies
Assumption 2 (Non linearity)
For the deterministic mapping \(F:{\mathscr{H}}\rightarrow {\mathscr{H}}\), we assume that there exists constant \(L\in (0,\infty )\) such that
As a consequence of (18) it holds that
Assumption 3 (Initial value)
We assume that \(X_{0}:\varOmega \rightarrow {\mathscr{H}}\) is a \(\mathcal {F}_{0}/{\mathscr{B}}({\mathscr{H}})\)-measurable mapping and \(X_{0}\in L^{2}\left (\varOmega , D\left (A^{\frac {2H+\beta -1}{2}}\right )\right )\).
In the Banach space \(\mathcal {D}\left (A^{\frac {\alpha }{2}}\right )\), \(\alpha \in \mathbb {R}\), we use notation \(\left \|A^{\frac {\alpha }{2}}\cdot \right \|=\|\cdot \|_{\alpha }\) and we recall the following properties of the semigroup S(t) generated by − A, which will be useful throughout this paper.
Proposition 3 (Smoothing properties of the semigroup)
[29] Let α > 0, δ ≥ 0 and 0 ≤ γ ≤ 1, then there exists a constant C > 0 such that
where l = 0,1 and \(D^{l}=\frac {d^{l}}{dt^{l}}\). If δ > γ then \(D(A^{\delta })\supset D(A^{\gamma })\). Moreover, AδS(t) = S(t)AδonD(Aδ).
The next lemma (specially (24) and (25)) is an important result which plays a crucial role to obtain regularity results, very useful in this work.
Lemma 3
For any 0 ≤ ρ ≤ 1, 0 ≤ γ ≤ 2 and 0 ≤ υ ≤ H with \(H\in \left (\frac {1}{2},1\right ]\), if the linear operator is given by (34), there exists a positive constant C such that for all 0 ≤ t1 ≤ t2 ≤ T,
Proof
See [27, Lemma 2.1] for the proof of (22) and (23). Concerning the proof of (24), the border case \(H=\frac {1}{2}\) if obtained using (22) with ρ = 1 and the order border case H = 1 is also obtained using (23) with γ = 2. Hence the proof of (24) is then completed by interpolation theory. The proof of (25) for 0 ≤ υ ≤ H is an immediate consequence of Proposition 3. The border case υ = H is proved by (24). This completes the proof of Lemma 3. □
Remark 2
Proposition 3 and Lemma 3 also hold with a uniform constant C (independent of h) when A and S(t) are replaced respectively by their discrete versions Ah and Sh(t) defined in Section 3, see, e.g., [16, 21].
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 (15) such that for all t ∈ [0,T], \(X(t)\in L^{2}\left (\varOmega ,D\left (A^{\frac {2H+\beta -1}{2}}\right )\right )\) with
Moreover, if the linear operator is given by (34), the following optimal regularity results in space and time hold
and for 0 ≤ t1 < t2 ≤ T;
where C = C(β,L,T,H) is a positive constant and β is the regularity parameter of Assumption 1.
Proof
[38, Theorem 3.5] gives the result of existence and uniqueness of the mild solution X. For regularity in space, we adapt from [26, 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 ≤ 2a2 + 2b2, (13) and (14), Assumption 1, Proposition 3, Lemma 3 (24) to have
For the proof of (29), triangle inequality yields
Using the stability property of the semigroup S(t) (20) with \(\gamma =\frac {2H+\beta -1}2\) and (28) allows to have
because \(D\left (A^{\frac {2H+\beta -1}{2}}\right )\) is continuously embedded in \(L^{2}\left (\varOmega \right )\) and
For the estimate of II, using triangle inequality, the estimate (a + b)2 ≤ 2a2 + 2b2, (13) and (14), inserting an appropriate power of A, Proposition 3, Assumption 1, Lemma 3 (25) with \(\upsilon =\frac {1-\beta }2\in [0,\frac {1}{2})\) (hence 0 ≤ υ < H), we obtain
Substituting (32) in (31) completes the proof of (29) and therefore that of Theorem 1. □
3 Numerical schemes
Throughout this section, we assume that Λ is bounded and has smooth boundary or is a convex polygon of \(\mathbb {R}^{d}\), d ∈{1,2,3}. In the rest of this paper we consider the SPDE (5) to be of the following form
x ∈Λ, t ∈ [0,T], where the function \(f : \varLambda \times \mathbb {R}\longrightarrow \mathbb {R}\) is continuously twice differentiable and the function \(b: \varLambda \times \mathbb {R}\longrightarrow \mathbb {R}\) is globally Lipschitz with respect to the second variable. In the abstract framework (5), the linear operator A takes the form
\(\mathbf {q}=\left (q_{i} \right )_{1 \leq i \leq d}\), where \(D_{ij}\in L^{\infty }(\varLambda )\), \(q_{i}\in L^{\infty }(\varLambda )\). We assume that there exists a positive constant c1 > 0 such that
The functions \(F : {\mathscr{H}} \longrightarrow {\mathscr{H}}\) and \( \phi : \mathbb {R}\longrightarrow HS\left (Q^{1/2}({\mathscr{H}}), {\mathscr{H}}\right )\) are defined by
for all x ∈Λ, \(v\in {\mathscr{H}}\), \(u\in Q^{1/2}({\mathscr{H}})\), with \({\mathscr{H}}=L^{2}(\varLambda )\). For an appropriate family of eigenfunctions (ei) such that \(\underset {i\in \mathbb {N}^{d}}{\sup }\left [\underset {x\in \varLambda }{\sup }\| e_{i}(x)\|\right ]<\infty \), it is well known [12, Section 4] that the Nemytskii operator F related to f and the operator ϕ associated to b defined in (36) satisfy Assumption 1 and Assumption 2. As in [8, 21] we introduce two spaces \(\mathbb {H}\) and V, such that \(\mathbb {H}\subset 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 V = H1(Λ)
where ∂v/∂A is the normal derivative of v and A is the exterior pointing normal n = (ni) 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 c0 and λ0 such that
By adding and substracting c0Xdt 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 c0X in a new nonlinear term that we still denote by F. The coercivity property (38) implies that − A is sectorial in L2(Λ), i.e. there exist \(C_{1}, \theta \in (\frac {1}{2}\pi ,\pi )\) such that
where \(S_{\theta }=\left \lbrace \lambda \in \mathbb {C} : \lambda =\rho e^{i \phi }, \rho >0, 0\leq \vert \phi \vert \leq \theta \right \rbrace \) (see [10]). Then, − A is the infinitesimal generator of a bounded analytic semigroup S(t) = e−tA on L2(Λ) such that
where \(\mathcal {C}\) denotes a path that surrounds the spectrum of − A. The coercivity property (38) 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 [10]). Under condition (35), it is well known (see, e.g., [8]) that the linear operator − A given by (34) generates an analytic semigroup S(t) ≡ e−tA. Following [8, 21], we characterize the domain of the operator Ar/2 denoted by \(\mathcal {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 [34, 37]. Let \(\mathcal {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}\subset {\mathscr{H}}\) denote the space of continuous functions that are piecewise linear over triangulation \(\mathcal {T}_{h}\). To discretise in space, we introduce the projection Ph from L2(Ω) to Vh define for u ∈ L2(Ω) by
The discrete operator \(A_{h}:V_{h}\rightarrow V_{h}\) is defined by
where a is the corresponding bilinear form of A.
Like the operator A, the discrete operator − Ah is also the generator of an analytic semigroup \(S_{h}(t):=e^{-tA_{h}}\). The semidiscrete space version of problem (5) is to find Xh(t) = Xh(⋅,t) such that for t ∈ [0,T]
The mild solution of (44) can be represented as follows
and we have the following regularity results.
Lemma 4
Assume that Assumptions 1–3 are satisfied, then the unique mild solution Xh(t) given by (45) satisfied
and for 0 ≤ t1 < t2 ≤ T;
Proof
Since the operators Ah and Sh(t) satisfy the same properties as A and S(t) (see Remark 2), then by using [34, (83)] and the boundedness of Ph in the proof of (28) and (29), we obtain the proof of the expression (46) and (47). The proof of Lemma 4 is thus completed. □
Now applying the linear implicit Euler method [13, 37] to (44) gives the following fully discrete scheme
Furthermore applying the stochastic exponential integrator ([21], SETD1) and Rosenbrock scheme ([26], SERS) to (44) yields
and
where \(\varphi _{1}(\varDelta tA_{h})=(\varDelta tA_{h})^{-1}\left (e^{\varDelta tA_{h}}-I\right )=\frac {1}{\varDelta t}{\int \limits }_{0}^{\varDelta t}e^{(\varDelta t-s)A_{h}} ds\), \({J^{h}_{m}}\) is the Frechet derivative of PhF at \({Z^{h}_{m}}\) and Sh,Δt := (I +Δ tAh)− 1. The term \({G^{h}_{m}}\) is the remainder at \({Z^{h}_{m}}\) and defines for all ω ∈Ω by
and
Note that the exponential integrator scheme (49) 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) (50) 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 [26, 27].
Assumption 4
For the deterministic mapping \(F:{\mathscr{H}}\rightarrow {\mathscr{H}}\), we also assume that there exists constant \(L\in (0,\infty )\) such that
3.1 Main result for SPDE driven by fBm
Theorem 2
Let X(tm) be the mild solution of (5) at time tm = mΔt, Δt ≥ 0 represented by (15). Let \({\zeta ^{h}_{m}}\) be the numerical approximations through (48) and (50) (\({\zeta ^{h}_{m}}={X^{h}_{m}}\) for implicit scheme, \({\zeta ^{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.
4 Proofs of the main result for SPDE with fBm
We introduce the Riesz representation \(R_{h}:V\rightarrow V_{h}\) defined by
under the regularity assumptions on the triangulation and in view of the V-ellipticity, it is well known ([8, 16]) that the following error bound holds:
for r ∈ [1,2]. Let us consider the following deterministic linear problem: Find u ∈ V such that
The corresponding semidiscrete problem in space consists to finding uh ∈ Vh such that
Let us define the following operator
Then we have the following lemma
Lemma 5
The following estimates hold for the semidiscrete approximation of (44). There exists a constant C > 0 such that
-
(i)
For v ∈ D(Aγ/2)
$$ \begin{array}{@{}rcl@{}} \|u(t)-u^{h}(t)\|=\|G_{h}(t)v\|\leq Ch^{r}t^{-(r-\gamma)/2}\|v\|_{\gamma}, r\in[0,2], \gamma\leq r \end{array} $$(61)for any t ∈ (0,T].
-
(ii)
For \(v\in D(A^{\frac {\gamma -1}2})\)
$$ \begin{array}{@{}rcl@{}} \left( {{\int}_{0}^{t}}\|G_{h}(s)v\|^{2}ds\right)^{\frac{1}{2}}\leq Ch^{\gamma}\|v\|_{\gamma-1},\qquad 0\leq \gamma\leq 2,\qquad t>0. \end{array} $$(62) -
(iii)
For \(v\in D(A^{-\frac {\rho }2})\)
$$ \begin{array}{@{}rcl@{}} \left\|{{\int}_{0}^{t}}G_{h}(s)vds\right\|\leq Ch^{2-\rho}\|v\|_{-\rho},\qquad 0\leq \rho\leq 1,\qquad t>0. \end{array} $$(63) -
(iv)
For \(v\in D(A^{\frac {\delta -1}2})\)
$$ \begin{array}{@{}rcl@{}} \left( {{\int}_{0}^{t}}\|G_{h}(s)v\|^{\frac{1}{H}}ds\right)^{2H}\leq Ch^{2(2H+\delta-1)}\|v\|^{2}_{\delta-1}, 0\leq \delta\leq 1, t>0. \end{array} $$(64)
Proof
See [34, Lemmas 3.1 and 3.2 (iv) and (v)] for the proofs of (i)–(iii). Let us prove (iv).
-
For \(H=\frac {1}{2}\), using (62) with γ = δ, we obtain
$$ \begin{array}{@{}rcl@{}} \left( {{\int}_{0}^{t}}\|G_{h}(s)v\|^{\frac{1}{H}}ds\right)^{2H}&=&{{\int}_{0}^{t}}\|G_{h}(s)v\|^{2}ds\\ &\leq & Ch^{2\delta}\|v\|^{2}_{\delta-1}=Ch^{2(2H+\delta-1)}\|v\|^{2}_{\delta-1}. \end{array} $$(65) -
For H = 1, using (63) with ρ = 1 − δ, we obtain
$$ \begin{array}{@{}rcl@{}} \left( {{\int}_{0}^{t}}\|G_{h}(s)v\|^{\frac{1}{H}}ds\right)^{2H}&=& \left( {{\int}_{0}^{t}}\|G_{h}(s)v\|ds\right)^{2}\\ &\leq & C\left( h^{2-(1-\delta)}\|v\|_{\delta-1}\right)^{2}\\ &=&Ch^{2+2\delta}\|v\|^{2}_{\delta-1}=Ch^{2(2H+\delta-1)}\|v\|^{2}_{\delta-1}. \end{array} $$(66)
Hence, the proof of (64) is thus completed by interpolation theory. □
Lemma 6 (Space error)
Let Assumptions 1–3 be fulfilled, then the following error estimate holds for the mild solution (15) and the discrete problem (45) holds
Proof
Using triangle inequality, we have
with
Note that the deterministic error is already estimated, so we will mostly concentrate our study on the stochastic error. Indeed, Lemma 5 with r = γ = 2H + β − 1 yields
Using triangle inequality, the boundedness of Sh(t − s) and Ph, Assumption 2 (more precisely (18)), we estimate the error e1(t) as follows
Applying Lemma 5 (i) with r = 2H + β − 1 and γ = 0, Assumption 2 (18), Theorem 1 (more precisely (29)) to the first term and Lemma 5 (iii) with ρ = 0, Theorem 1 (more precisely (26)) and Lemma 4 to the second term yields
For the estimation of e2(t), triangle inequality, the estimate (a + b)2 ≤ 2a2 + 2b2, (13) and (14), Lemma 5 ((i) with r = 2H + β − 1, γ = β − 1, (iv) with δ = β) and Assumption 1 yields
Combining the estimates (68), (70), (71) and applying Gronwall inequality ends the proof. □
4.1 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 [37, (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 7
For any m, h and Δt the following estimates holds
-
(i)
$$ \begin{array}{@{}rcl@{}} \left\|(I+\varDelta t A_{h})^{-m}\right\|\leq 1. \end{array} $$(72)
-
(ii)
For all \(u\in D(A^{\frac {\gamma -1}2})\), 0 ≤ γ ≤ 1,
$$ \begin{array}{@{}rcl@{}} \left\|(S^{m}_{h,\varDelta t}-S_{h}(t_{m}))P_{h}u\right\|\leq C\varDelta t^{\frac{2H+\gamma-1}2}t_{m}^{-H}\|u\|_{\gamma-1}. \end{array} $$(73) -
(iii)
For all \(u\in D(A^{\frac {\gamma -1}2})\), 0 ≤ γ ≤ 1,
$$ \begin{array}{@{}rcl@{}} \left\|(S^{m}_{h,\varDelta t}-S_{h}(t))P_{h}u\right\|\leq C\varDelta t^{\frac{2H+\gamma-1}2}t^{-H}\|u\|_{\gamma-1},t\in[t_{j-1},t_{j}), \end{array} $$(74)for any j = 1,2,⋅⋅⋅,M.
-
(iv)
If u ∈ D(Aμ/2), 0 ≤ μ ≤ 2, then
$$ \begin{array}{@{}rcl@{}} \|(S^{m}_{h,\varDelta t}-S_{h}(t_{m}))P_{h}u\|\leq C\varDelta t^{\mu/2}\|u\|_{\mu}. \end{array} $$(75) -
(v)
For all non smooth data \(u\in {\mathscr{H}}\),
$$ \begin{array}{@{}rcl@{}} \|(S^{m}_{h,\varDelta t}-S_{h}(t_{m}))P_{h}u\|\leq C\varDelta tt_{m}^{-1}\|u\|. \end{array} $$(76)
Proof
See [34, Lemma 3.3] for the proof of (i), (iv) and (v). For the proof of (ii), we use[34, Lemma 3.3 (iii) (88)] as follows
Hence substituting (76) in [34, (84)] completes the proof of (ii). Now for the proof of (iii), triangle inequality, Lemma 7 (ii), the property of discrete semigroup and [34, (83)] yields
The proof of Lemma 7 is thus completed. □
Lemma 8
-
(i)
For any ρ ∈ [0,1] and u ∈ D(A−ρ/2) there exists a positive constant C such that
$$ \begin{array}{@{}rcl@{}} \left\|\sum\limits_{j=1}^{m}{\int}_{t_{j-1}}^{t_{j}}(S^{j}_{h,\varDelta t}-S_{h}(s))P_{h}uds\right\|\leq C\varDelta t^{\frac{2-\rho}2-\epsilon}\|u\|_{-\rho}. \end{array} $$(79) -
(ii)
For any μ ∈ [0,1] and \(u\in D(A^{\frac {\mu -1}2})\) the following estimate holds
$$ \begin{array}{@{}rcl@{}} \left( \sum\limits_{j=1}^{m}{\int}_{t_{j-1}}^{t_{j}}\|(S^{j}_{h,\varDelta t}-S_{h}(s))P_{h}u\|^{\frac{1}{H}}ds\right)^{2H}\leq C\varDelta t^{2H+\mu-1-\epsilon}\|u\|_{\mu-1}^{2}. \end{array} $$(80)
Where 𝜖 is an arbitrary small number.
Proof
See [34, Lemma 3.5] for the proof of (i). For the proof of (ii), we have
Using triangle inequality and the estimate (a + b)k ≤ 2k− 1ak + 2k− 1bk (with k ≥ 1 and a,b ≥ 0) we obtain
By Lemma 7 (ii) with γ = μ we obtain
By inserting an appropriate power of Ah, [34, (81)] and Remark 2 yields
Substituting (83) and (84) in (82) yields
Concerning the estimate of K2, let 𝜖 > 0 small enough, Lemma 7 (iii) with γ = μ yields
hence
This completes the proof of Lemma 8. □
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 err0 is estimated by Lemma 6. It remains to estimate the time error err1. We recall that the exact solution at tm of the semidiscrete problem (44) is given by
We also recall that the numerical solution at tm given by (48) can be rewritten as
where the notation [t], [t]m are defined by
It follows from (89) and (90) that
As we said at the beginning of this section, following closely the work done in [37, (4.7)–(4.14)] and replacing its preparatory results with Lemma 7 (i), (iv) with μ = 2H + β − 1, (v), Lemma 8 (i) with ρ = 0, Remark 2 (20) with δ = γ = 1, Assumptions 2-3, boundedness of Sh(tm − s) and Ph, the stability properties of a discrete semigroup Sh(t), (18) and (47), we have
Note that in this work, we do not need to impose an assumption on \(F^{\prime \prime }\) to increase the convergence rate as it is done in [37]. Indeed, thanks to (47) the following estimate is largely sufficient to reach a higher rate.
Let us focus now on the estimate I2, using triangle inequality and the estimate (a + b)2 ≤ 2a2 + 2b2 we split it in three terms
Firstly using (14), inserting an appropriate power of Ah, [34, (81)], Assumption 1 (more precisely (17)) and Remark 2 ((22) with ρ = 1 − β) we obtain
Secondly (14), the change of variable j = m − k and ς = tm − s, Lemma 7 (iii) with γ = β and Assumption 1 (more precisely (17)) yields
Thirdly using (13), the change of variable j = m − k and ς = tm − s, Lemma 8 (ii) with μ = β and Assumption 1 (more precisely (16)) we obtain
Hence, inserting (95)–(97) in (94) and taking the square-root give
Adding (93) and (99) we obtain
Applying the discrete version of the Gronwall inequality yields
Adding (67) and (100) completes the proof.
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 \(\mathcal {O}(\varDelta t)\) is reached.
Corollary 1
Let X(tm) be the mild solution of (5) (A self-adjoint) at time tm = mΔt, Δt ≥ 0 represented by (15). Let \({X^{h}_{m}}\) be the numerical approximation through (48). 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 8 in the self-adjoint case. The result is presented in the following lemma:
Lemma 9
-
(i)
For any ρ ∈ [0,1] and u ∈ D(A−ρ/2) there exists a positive constant C such that
$$ \begin{array}{@{}rcl@{}} \left\|\sum\limits_{j=1}^{m}{\int}_{t_{j-1}}^{t_{j}}(S^{j}_{h,\varDelta t}-S_{h}(s))P_{h}uds\right\|\leq C\varDelta t^{\frac{2-\rho}2}\|u\|_{-\rho}. \end{array} $$(102) -
(ii)
For any μ ∈ [0,1] and \(x\in D(A^{\frac {\mu -1}2})\) the following estimate holds
$$ \begin{array}{@{}rcl@{}} &&\sum\limits_{i,j=0}^{m-1}{\int}_{t_{j}}^{t_{j+1}}{\int}_{t_{i}}^{t_{i+1}}\left\langle\mathcal{S}_{h}(u,t_{i})P_{h}x,\mathcal{S}_{h}(v,t_{j})P_{h}x\right\rangle \kappa(u,v)dudv\\ &\leq & C\varDelta t^{2H+\mu-1}\|x\|_{\mu-1}^{2}, \end{array} $$(103)and
$$ \begin{array}{@{}rcl@{}} &&\sum\limits_{i,j=0}^{m-1}{\int}_{t_{j}}^{t_{j+1}}{\int}_{t_{i}}^{t_{i+1}}\left\langle\mathcal{T}_{h}(i)P_{h}x,\mathcal{T}_{h}(j)P_{h}x\right\rangle \kappa(u,v)dudv\\ &\leq & C\varDelta t^{2H+\mu-1}\|x\|_{\mu-1}^{2}, \end{array} $$(104)
where
Proof
See [15, Proof of Lemma 4.4 (i)] for the proof of (i) and [38, Lemmas 4.8 and 4.9], [34, (83)] for the proof of (ii). □
With this new lemma, we are now in position to prove our Corollary 1.
Proof Proof of Corollary 1
Recall that the time error err1 is defined as
Following closely the work done in [37, (4.7)–(4.14)] and replacing its preparatory results with Lemma 7 (i), (iv) with μ = 2H + β − 1, (v) with σ = 0, Lemma 9 (i) with ρ = 0, Remark 2 (20) with δ = γ = 1, Assumptions 2–3, boundedness of Sh(tm − s) and Ph, the stability properties of a discrete semigroup Sh(t), (18) and (47), we have
Concerning the estimate I2, we also split it in three terms as in (95). The estimates \(I_{21}^{2}\) and \(I_{22}^{2}\) are still the same but we need to re-estimate \(I_{23}^{2}\). In this fact, since the sequence of random variables
\(\left ({\int \limits }_{0}^{t_{m}}\left (S_{h}(t_{m}-s)-S_{h,\varDelta t}^{m-[s]^{m}}\right ) P_{h}\phi (t_{m-1}) Q^{\frac {1}{2}}e_{i}d{\beta _{i}^{H}}(s), i\in \mathbb {N}^{d}\right )\) is mutually independent Gaussian random variable, using the estimate (a + b)2 ≤ 2a2 + 2b2, notation (105), Assumption 1 (more precisely (16)), Lemma 9 (ii) with μ = β, we obtain
hence inserting (96), (97) and (108) in (95) and taking the square-root give
Adding (107) and (109) we obtain
Applying the discrete version of the Gronwall inequality yields
Adding (67) and (110) completes the proof. □
4.2 Proof of Theorem 2 for SETD1
As usual, splitting the fully discrete error in two terms yields
Since the space error err0 has been estimated by Lemma 6, we only need to estimate the time error err2. Remember that the exact solution at tm is given by
and we recall that the numerical solution at tm given by (49) can be rewritten as
where the notations [t] and [t]m are given by (91). By (111) and (112), we have
Applying the triangle inequality yields
Using the boundedness of Ph and Sh(tm − s), Lemma 3 and (47), we easily have
and
Adding (115) and (116), we obtain
We estimate at now \(I_{2}^{\prime }\). Using triangle inequality and the estimate (a + b)2 ≤ 2a2 + 2b2, we split it in three terms
Thanks to (96) we have
Thereafter, (14), the change of variable j = m − k and ς = tm − s, inserting an appropriate power of Ah, [34, (81)] and Remark 2 (more precisely (20) with \(\gamma =\frac {2H+\beta -1}2\), (21) and Assumption 1 (more precisely (17)) yields
Finally, using (13), inserting an appropriate power of Ah, [34, (81)] and Remark 2 (more precisely (20) with \(\gamma =\frac {2H+\beta -1}2\), (21) and (24)), Assumption 1 (more precisely (16)) we obtain
hence inserting (119)–(121) in (118) and taking the square-root give
Adding (117) and (122) we obtain
Using the discrete version of the Gronwall inequality yields
Combining (67) and (124) completes the proof.
4.3 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 [26] we obtain
Lemma 10
[26, Lemma 5] For all \(M\in \mathbb {N}\) and all ω ∈Ω, there is a positive constant \(C^{\prime }_{1}\) independent of h, m, Δt and the sample ω such that
Lemma 11
[26, Lemma 6] The function \({G^{h}_{m}}(\omega )\) defined by (52) satisfies the global Lipschitz condition with a uniform constant \(C^{\prime }>0\), independent of h, m and ω such that
Lemma 12
[26, Lemma 9] For all ω ∈Ω, the stochastic perturbed semigroup \({S^{h}_{m}}(\omega )(t):=e^{(-A_{h}+{J^{h}_{m}}(w))t}\) satisfies the following properties
-
(i)
For γ1,γ2 ≤ 1 such that 0 ≤ γ1 + γ2 ≤ 1,
$$ \begin{array}{@{}rcl@{}} \left\|A_{h}^{-\gamma_{1}}({S^{h}_{m}}(\omega)(t)-I)A_{h}^{-\gamma_{2}}\right\|_{L(\mathcal{H})}\leq C t^{\gamma_{1}+\gamma_{2}},t\in(0,T]. \end{array} $$ -
(ii)
For γ1 ≥ 0 we have
$$ \begin{array}{@{}rcl@{}} \left\|{S^{h}_{m}}(\omega)(t)A_{h}^{\gamma_{1}}\right\|_{L(\mathcal{H})}\leq C t^{-\gamma_{1}},t\in(0,T]. \end{array} $$ -
(iii)
For γ1 ≥ 0 and 0 ≤ γ2 < 1such that γ2 − γ1 ≥ 0, we have
$$ \begin{array}{@{}rcl@{}} \|A_{h}^{-\gamma_{1}}{S^{h}_{m}}(\omega)(t)A_{h}^{\gamma_{2}}\|_{L(\mathcal{H})}\leq C t^{\gamma_{1}-\gamma_{2}},t\in(0,T]. \end{array} $$ -
(iv)
For γ1,γ2 > 0 such that 0 ≤ γ1 − γ2 ≤ 1, then the following estimate holds
$$ \begin{array}{@{}rcl@{}} \|A_{h}^{-\gamma_{1}}({S^{h}_{m}}(\omega)(t)-I)A_{h}^{\gamma_{2}}\|_{L(\mathcal{H})}\leq C t^{\gamma_{1}-\gamma_{2}},t\in(0,T]. \end{array} $$
Lemma 13
[26, Lemma 10] The stochastic perturbed semigroup \({S^{h}_{m}}(\omega )\) 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 6 we have the estimate of the space error
Consider now the estimate of the time error err3. We recall that the semidiscrete problem (44) can be rewritten as
for all tm ≤ t ≤ tm+ 1 where \({J^{h}_{m}}\) and \({G^{h}_{m}}\) is given by (51) and (52). Hence, the exact solution at time tm of the semidiscrete problem (126) is given by
and the numerical solution (50) can be rewritten as
If m = 1 then from (127) and (128) we obtain
By [26, (93)] we have the estimate
Using triangle inequality and the estimate (a + b)2 ≤ 2a2 + 2b2 we split it into three terms
For the estimate \(II_{1}^{\prime \prime 2}\), using (14), inserting an appropriate power of Ah, [34, (81)], Assumption 1 (more precisely (17)) and Lemma 12 (ii) with \(\gamma _{1}=\frac {1-\beta }2\) we obtain
We denote by 𝜖 a positive constant small enough, using (13), inserting an appropriate power of Ah, Assumption 1, [34, (83)], Lemma 12 (ii) with \(\gamma _{1}=H-\frac {\epsilon }2\), (iv) with \(\gamma _{1}=H-\frac {\epsilon }2\) and \(\gamma _{2}=\frac {1-\beta }2\) if 0 ≤ β < 1 (or (i) with \(\gamma _{1}=H-\frac {\epsilon }2\) and γ2 = 0 if β = 1) we have
Hence, putting (133) and (133) in (131) and taking the square-root give
For m ≥ 2, we recall that the solution at tm of the semidiscrete problem (126) is given by
We recall also that the numerical solution at tm given by (128) can be rewritten as
Using (136), (137) and the triangle inequality, we have
where
and
By Lemma 10, triangle inequality and Lemma 4 (47), we easily have
In a similar way, using Lemma 13 with ν = 0, Lemma 10, Lemma 11, triangle inequality and Lemma 4 (47), we obtain
For the estimate \(IV^{\prime \prime }\), using triangle inequality and the estimate (a + b)2 ≤ 2a2 + 2b2 we split it in two terms
Firstly, using (14), inserting an appropriate power of Ah, [34, (81)], Assumption 1 (more precisely (17)) and Lemma 12 (ii) with \(\gamma _{1}=\frac {1-\beta }{2}\) we obtain
Secondly by (13), inserting an appropriate power of Ah, Assumption 1, [34, (81)], Lemma 12 (ii) with \(\gamma _{1}=H-\frac {\epsilon }{2}\), (iv) with \(\gamma _{1}=H-\frac {\epsilon }{2}\) and \(\gamma _{2}=\frac {1-\beta }{2}\), we obtain
Hence, substituting (142) and (143) in (141) and taking the square-root give
For estimate \(VI^{\prime \prime }\), using triangle inequality and the estimate (a + b + c)2 ≤ 3a2 + 3b2 + 3c2 we split it in two terms as
Let 𝜖 be a sufficient small number. At first, using (14), inserting an appropriate power of Ah, [34, (81)], Assumption 1 (more precisely (17)), Lemma 13 (ii) with \(\nu =H-\frac {\epsilon }{2}\), Lemma 12 (iii) with \(\gamma _{1}=\gamma _{2}=\frac {1-\beta }{2}\), the variable change j = m − k − 1 and [26, (169)] we have
Using (14), inserting an appropriate power of Ah, Lemma 13 with ν = H, Lemma 12 (iii) with γ1 = γ2 = H, (iv) with γ1 = H and \(\gamma _{2}=\frac {1-\beta }{2}\) if 0 < β < 1 (or (i) with γ1 = H and γ2 = 0 if β = 1), [34, (81)], Assumption 1 (more precisely (17)) the variable change j = m − k − 1 and [26, (169)] yields
Afterwards, (13), inserting an appropriate power of Ah, Lemma 13 with \(\nu =H-\frac {\epsilon }2\), Lemma 12 (iii) with \(\gamma _{1}=\gamma _{2}=H-\frac {\epsilon }2\), (iv) with \(\gamma _{1}=H-\frac {\epsilon }2\) and \(\gamma _{2}=\frac {1-\beta }{2}\) if 0 < β < 1 (or (i) with \(\gamma _{1}=H-\frac {\epsilon }2\) and γ2 = 0 if β = 1), [34, (81)], Assumption 1 (more precisely (16)) the variable change j = m − k − 1 and [26, (169)] yields
Hence, inserting (146), (147) and (148) in (145) and taking the square-root give
Adding (135), (139), (140), (144), (149) and applying Gronwall’s lemma yields
5 Extension to SPDE driven simultaneously by fBm and Poisson random measure
5.1 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 \({\mathscr{H}}=L^{2}(\varLambda )\), (151) is equivalent to (1) where the linear operator A and the nonlinear function F are defined as in (34) and (36). 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 [38]. 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
with \({Z^{h}_{0}}=P_{h}X_{0}\), where
To obtain the optimal order in time, as in [28], we need the following assumption in Poisson measure noise.
Assumption 5
The covariance operator \(Q:{\mathscr{H}}\longrightarrow {\mathscr{H}}\) and the jump coefficient satisfy the following estimate
where η = 2H + β − 1 with β ∈ (0,1] as in Assumptions 3 and 1.
Remark 3
All the regularity results in space and time, both for continuous (1) (or semidiscrete equation) important to achieve optimal convergence orders can easily be extended from our results on fBm in Theorem 1 and Lemma 4 by just following [28, Proposition 3.1].
5.2 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(tm) be the mild solution of (1) at time tm = mΔt, Δt ≥ 0. Let \({\zeta ^{h}_{m}}\) be the numerical approximations through (153) and (128) (\({\zeta ^{h}_{m}}={X^{h}_{m}}\) for implicit scheme, \({\zeta ^{h}_{m}}={Z^{h}_{m}}\) for SERS) and \({Y^{h}_{m}}\) the numerical approximation through the SETD1 given in (152). If Assumptions 1–4 and Assumption 5 hold with β ∈ (0,1], then
and
where 𝜖 is a positive constant small enough.
Corollary 2
Let X(tm) be the mild solution of (1) (A self-adjoint) at time tm = mΔt, Δt ≥ 0. Let \({X^{h}_{m}}\) be the numerical approximations through (153). If Assumptions 1–3 and Assumption 5 hold with β ∈ (0,1], then
5.3 Proof of convergence results for SPDE with fBm and Poisson measure noise
As in [28], the proofs are based on Burkholder-Davis-Gundy Inequality where the fBm version is given in [32, Theorem 1.2]. Under Assumptions 1–3 and 5, The regulatity result in time is
Where C = C(β,L,T,H) is a positive constant and β is the regularity parameter of Assumption 1. In all our schemes, the error can be splitted in space error err0 and the time error err1. The space error err0 can be estimated as in Lemma 6 using results from the proof of [28, Theorem 4.1]. More precisely in the estimation of the poisson term in their case, we replace β ∈ [0,2] by η = 2H + β − 1 defined in Assumption 5. The time error err1 will be here splitted in three terms. More precisely the deterministic I1 related to the nonlinear function F, the fBm term I2 and the Poisson term I3. The estimation of I1 is done with the aid of Assumption 4 and (160) similarly as the work done in [28, (77)–(88)] for implicit and exponential schemes. As in the proof of [26, Theorem 10], we use the Taylor expansion in Banach space (see [28, (77)]) to estimate I1 for SERS. The fBm term I2 is done exactly as in the previous section for the scheme without Poisson. By replacing β ∈ [0,2] by η = 2H + β − 1 defined in Assumption 5, the estimation of the Poisson term I3 is done as in [28, Theorem 5.2] for implicit and SETD1 schemes using Burkholder-Davis-Gundy Inequality, the work in [26, Theorem 10] and preparatory results Lemmas 10-13 for SERS.
6 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 [13] and the circulant method [32], which will be used in this work to generate the fBm. Here we consider the stochastic advection diffusion reaction SPDE (5)–(33) with constant diagonal diffusion tensor \( \mathbf {D}= 10^{-2} \mathbf {I}_{2}=(D_{i,j}) \) in (34), and mixed Neumann-Dirichlet boundary conditions on Λ = [0,L1] × [0,L2]. The Dirichlet boundary condition is X = 1 at Γ = {(x,y) : x = 0} and we use the homogeneous Neumann boundary conditions elsewhere. This example can be seen as an engineering application where our SPDE (5)–(33) models the uncertainties in chemical reaction rates for stochastic transport of pollutants in porous media due by diffusion-advection-reaction with long-range dependence. If we include the Poisson noise term, our (1)–(151) should model the uncertainties in chemical reaction rates for stochastic transport with long long-range dependence and some rare high chemical reaction rate, which can occur in the fracture networks. If we include the dispersion tensor as in [9], our (5)–(33) models also the sub-grid fluctuations in transport velocities and uncertainties in the reaction rates with long long-range dependence. In all these cases, the unknown function X is the concentration of the pollutants.
The eigenfunctions \( \{e_{i,j} \} =\{e_{i}^{(1)}\otimes e_{j}^{(2)}\}_{i,j\geq 0} \) of the covariance operator Q are taken to be the same as for Laplace operator −Δ with homogeneous boundary condition, so we have
where \(l \in \left \lbrace 1, 2 \right \rbrace , x\in \varLambda \). In the noise representation (4), we have used
for some small δ > 0. We have used b(x,t) = 2 in (33), so ϕ in Assumption 1 is obviously satisfied for β = (0,1]. In our simulations, we have used δ = 0.001. The function f used in (36) to be \(f(x,z)= \frac {z}{1+z}\) for all \((x,z)\in \varLambda \times \mathbb {R}\). Therefore the corresponding Nemytskii operator F defined by (36) obviously satisfies Assumption 2. We obtain the Darcy velocity field q = (qi) by solving the following system
with Dirichlet boundary conditions on \({\varGamma _{D}^{1}}=\left \lbrace 0,L_{1} \right \rbrace \times \left [ 0,L_{2}\right ] \) and Neumann boundary conditions on \({\varGamma _{N}^{1}}=\left (0,L_{1}\right )\times \left \lbrace 0,L_{2}\right \rbrace \) such that
and −k∇p(x,t) ⋅n = 0 in \({\varGamma _{N}^{1}}\). Note that k is the permeability tensor and p the presure. We use a random permeability field as in [36, Figure 6]. The streamline of the velocity field q are given in Fig. 1d. To deal with high Péclet number, we discretise in space using finite volume method, viewed as a finite element method (see [33]). We take L1 = 3 and L2 = 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 X0 = 0, so we can therefore expect high-order convergence, which depends only on the noise term and H.
Convergence in the root mean square L2 norm at T = 1 as a function of Δt for implicit scheme (a), exponential scheme (b) and exponential Rosenbrock scheme (c). We have used here 50 realizations. The streamline of the velocity field q is given in (d)
Figure 1a is the errors graph for the implicit scheme with different values of H. We have observed that the order of convergence is 0.48 in time for H = 0.51 and β = 1, 0.6476 for H = 0.65 and β = 1.
Figure 1b is the errors graph for the exponential scheme with two values of H. We have observed the order of convergence is 0.5012 in time for H = 0.51 and β = 1, 0.6653 for H = 0.65 and β = 1.
Figure 1c 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.65 and H = 0.51. Here we have fixed β = 1 and same Gaussian random 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.
Numerical samples solution using exponential scheme with H = 0.65 in (a) and H = 0.51 in (b)
Notes
In this case the SPDE (5) is said to be driven by its linear part.
Linear operator A self-adjoint for implicit scheme
References
Alberverio, S., Gawarecki, L., mandrekar, V., Rüdiger, B., Sarkar, B.: Ito formula for mild solutions of SPDEs with Gaussian and non-Gaussian noise and applications to stability propertieŝ. Random Oper. Stoch. Equ. 25(2), 79–105 (2017)
Alós, E., Nualart, D.: Stochastic integration with respect to the fractional Brownian motion. Stoch. Stoch. Rep. 75, 129–152 (2003)
Caraballo, T., Garrido-Atienza, M. J., Taniguchi, T.: The existence and exponential behavior of solutions to stochastic delay evolution equations with a fractional Brownian motion. Nonlinear Anal. 74, 3671–3684 (2011)
Cheridito, P.: Arbitrage in fractional Brownian motion models. Finance Stoch. 7, 533–553 (2003)
Cont, R., Tankov, P.: Financial modelling with jump process. In: Financial Mathematics Series. CRC Press, Boca Raton (2000)
Duncan, T. E., Maslowski, B., Pasik-Duncan, B.: Fractional Brownian motion and stochastic equations in Hilbert spaces. Stoch. Dyn. 2, 225–250 (2002)
Duncan, T. E., Maslowski, B., Pasik-Duncan, B.: Semilinear stochastic equations in a Hilbert spaces with a fractional Brownian motion. SIAM J. Math. Anal. 40, 2286–2315 (2009)
Fujita, F., Suzuki, T.: Evolution problems (Part 1). Handbook of Numerical Analysis. In: Ciarlet, P.G., Lions, J.L. (eds.) , vol. 2, pp 789–928. North-Holland, Amsterdam (1991)
Geiger, S., Lord, G., Tambue, A.: Exponential time integrators for stochastic partial differential equations in 3D reservoir simulation. Comput. Geosci. 16, 323–334 (2012)
Henry, D.: Geometric Theory of Semilinear Parabolic Equations. Lecture Notes in Mathematics, vol. 840. Springer, Berlin (1981)
Hu, Y., Oksendal, B.: Fractional white noise calculus and applications to finance. Infin. Dimens. Quantum Probab. Relat. Top. 6, 1–32 (2003)
Jentzen, A., Röckner, M.: Regularity analysis for stochastic partial differential equations with nonlinear multiplicative trace class noise. J. Diff. Equ. 252(1), 114–136 (2012)
Kamrani, M., Jamshidi, N.: Implicit Euler approximation of stochastic evolution equations with fractional Brownian motion. Commun. Nonlinear Sci. Numer. Simul. 44, 1–10 (2017)
Kim, Y. T., Rhee, J. H.: Approximation of the solution of stochastic evolution equation with fractional Brownian motion. J. Korean Stat. Soc. 33, 459–470 (2004)
Kruse, R.: Optimal error estimates of Galerkin finite element methods for stochastic partial differential equations with multiplicative noise. IMA J. Numer. Anal. 34, 217–251 (2014)
Larsson, S.: Nonsmooth data error estimates with applications to the study of the long-time behavior of finite element solutions of semilinear parabolic problems. Preprint 36, Department of Mathematics, Chalmers University of Technology. Available online at http://citeseerx.ist.psu.edu/viewdoc/summary-doi=10.1.1.28.1250 (1992)
Leland, W., Taqqu, M., Willinger, W., Wilson, D.: On the self-similar nature of ethernet traffic (Extended version). IEEE/ACM Trans. Netw. 2, 1–15 (1994)
Li, Z.: Global attractiveness and quasi-invariant sets of impulsive neutral stochastic functional differential equations driven by fBm. Neurocomputing 177, 620–627 (2016)
Li, Z., Yan, L.: Ergodicity and stationary solution for stochastic neutral retarded partial differential equations driven by fractional Brownian motion. J. Theor. Probab. 32(3), 1399–1419 (2019)
Li, Z., Yan, L.: Stochastic averaging for two-time-scale stochastic partial differential equations with fractional Brownian motion. Nonlinear Anal.: Hybrid Syst. 31, 317–333 (2019)
Lord, G. J., Tambue, A.: Stochastic exponential integrators for the finite element discretization of SPDEs for multiplicative & additive noise. IMA J. Numer. Anal. 2, 515–543 (2013)
Mandelbrot, B. B.: The variation of certain speculative prices. J. Bus. 36, 394–419 (1963)
Maslowsk, B., Nualart, D.: Evolution equations driven by a fractional Brownian motion. J. Funct. Anal. 202, 277–305 (2003)
Mémina, J., Mishura, Y., Valkeila, E.: Inequalities for the moments of Wiener integrals with respect to a fractional Brownian motion. Stat. Probab. Lett. 51, 197–206 (2001)
Mishura, Y.: Stochastic Calculus for Fractional Brownian Motion and Related Processes. Lecture Notes in Mathematics, vol. 1909. Springer, Berlin (2008)
Mukam, J. D., Tambue, A.: Strong convergence analysis of the stochastic exponential Rosenbrock scheme for the finite element discretization of semilinear SPDEs driven by multiplicative and additive noise. J. Sci. Comput. 74, 937–978 (2018)
Mukam, J. D., Tambue, A.: A note on exponential Rosenbrock-Euler method for the finite element discretization of a semilinear parabolic partial differential equation. Comput. Math. Appl. 76, 1719–1738 (2018)
Mukam, J. D., Tambue, A.: Optimal strong convergence rates of numerical methods for semilinear parabolic SPDE driven by Gaussian noise and Poisson random measure. Comput. Math. Appl. 77(10), 2786–2803 (2019)
Pazy, A.: Semigroups of Linear Operators and Applications to Partial Differential Equations. Applied Mathematical Sciences, p 44. Springer, New York (1983)
Platen, E., Bruti-Liberati, N.: Numerical solution of stochastic differential equations with jumps in finance. In: Stochastic Modelling and Applied Probability, vol. 64. Springer, Berlin (2010)
Prévôt, C., Röckner, M.: A Concise Course on Stochastic Partial Differential Equations. Lecture Notes in Mathematics, vol. 1905. Springer, Berlin (2007)
Shevchenko, G.: Fractional Brownian Motion in a Nutshell. Int. J. Mod. Phys.: Conf. Ser. 36, 1560002 (2015). (16 pages)
Tambue, A.: An exponential integrator for finite volume discretization of a reaction-advection-diffusion equation. Comput. Math. Appl. 71(9), 1875–1897 (2016)
Tambue, A., Mukam, J. D.: Strong convergence of the linear implicit Euler method for the finite element discretization of semilinear SPDEs driven by multiplicative or additive noise. Appl. Math. Comput. 346, 23–40 (2019)
Tambue, A., Ngnotchouye, J. M. T.: Weak convergence for a stochastic exponential integrator and finite element discretization of stochastic partial differential equation with multiplicative & additive noise. Appl. Num. Math. 108, 57–86 (2016)
Tambue, A., Lord, G. J., Geiger, S.: An exponential integrator for advection-dominated reactive transport in heterogeneous porous media. J. Comput. Phys. 229(10), 3957–3969 (2010)
Wang, X.: Strong convergence rates of the linear implicit Euler method for the finite element discretization of SPDEs with additive noise. IMA J. Numer. Anal. 37(2), 965–984 (2017)
Wang, X., Qi, R., Jiang, F.: Sharp mean-square regularity results for SPDEs with fractional noise and optimal convergence rates for the numerical approximations. BIT Numer. Math. 57(2), 557–585 (2017)
Willinger, W., Taqqu, M. S., Leland, W. E., Wilson, D. V.: Self-similarity in high-speed packet traffic: analysis and modelling of Ethernet traffic measurements. Stat. Sci. 10, 67–85 (1995)
Acknowledgements
Aurelien Junior Noupelah thanks Prof Louis Aime Fono and Prof Jean Louis Woukeng for their constant supports. We would like to thank Jean Daniel Mukam for very useful discussions.
Funding
Open Access funding provided by Western Norway University Of Applied Sciences
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Noupelah, A.J., Tambue, A. 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. Numer Algor 88, 315–363 (2021). https://doi.org/10.1007/s11075-020-01041-1
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11075-020-01041-1
Keywords
- Stochastic parabolic partial differential equations
- Fractional Brownian motion
- Finite element method
- Errors estimate
- Finite element methods
- Timestepping methods
- Poisson random measure