Strong Convergence of a Stochastic Rosenbrock-type Scheme for the Finite Element Discretization of Semilinear SPDEs Driven by Multiplicative and Additive Noise

Strong Convergence of a Stochastic Rosenbrock-type Scheme for the Finite Element Discretization of Semilinear SPDEs Driven by Multiplicative and Additive Noise

Abstract

This paper aims to investigate the numerical approximation of a general second order parabolic stochastic partial differential equation(SPDE) driven by multiplicative and additive noise. Our main interest is on such SPDEs where the nonlinear part is stronger than the linear part, usually called stochastic dominated transport equations. Most standard numerical schemes lose their good stability properties on such equations, including the current linear implicit Euler method. We discretise the SPDE in space by the finite element method and propose a new scheme in time appropriate for such equations, called stochastic Rosenbrock-Type scheme, which is based on the local linearisation of the semi-discrete problem obtained after space discretisation. We provide a strong convergence of the new fully discrete scheme toward the exact solution for multiplicative and additive noise. Our convergence rates are in agreement with results in the literature. Numerical experiments to sustain our theoretical results are provided.

Keywords:
Rosenbrock-type scheme Stochastic partial differential equations Multiplicative & Additive noise Strong convergence Finite element method.
Msc:
MSC 65C30 MSC 74S05 MSC 74S60

1 Introduction

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

(1)

for all , on the Hilbert space , the final time, and are nonlinear functions, is the initial data which is random, is a linear operator, unbounded, not necessary self-adjoint, and is assumed to be a generator of a contraction semigroup , The noise is a Wiener process defined in a filtered probability space . The filtration is assumed to fulfill the usual conditions (see ((34), Definition 2.1.11)). We assume that the noise can be represented as

(2)

where , , are respectively the eigenvalues and the eigenfunctions of the covariance operator , and are independent and identically distributed standard Brownian motions. Precise assumptions on , , and will be given in the next section to ensure the existence of the unique mild solution of (1), which has the following representation (see (32); (34))

(3)

for all . Equations of type (1) are used to model different real world phenomena in different fields such as biology, chemistry, physics etc (35); (36); (5). In many cases explicit solutions of SPDEs are unknown, therefore numerical approximations are the only tool appropriate to approach them. Numerical approximation of SPDE of type (1) is therefore an active research area and has attracted a lot of attentions since two decades, see e.g. (22); (43); (44); (18); (11); (12); (46); (45); (35); (33); (16) and references therein. Due to the time step restriction of the explicit Euler method, linear implicit Euler method is used in many situations. Linear implicit Euler method has been investigated in the literature, see e.g. (18); (44); (23); (40) and the references therein. The resolvent operator plays a key role to stabilise the linear implicit Euler method, where is the discrete form of after space discretization. Such approach is justified when the linear operator is strong. Indeed, when is stronger than , the operator drives the SPDE (1) and the good stability properties of the linear implicit Euler method and exponential integrators are guaranteed. In more concrete applications, the nonlinear function can be stronger. Typical examples are stochastic reaction equations with stiff reaction term. For such equations, both linear implicit Euler method (18); (44); (23); (40) and exponential integrators (22); (43); (11) behave like the standard explicit Euler method (see Section 2.3) and therefore lose their good stabilities properties. For such problems in the deterministic context, Exponential Rosenbrock-type (39); (10) methods and Rosenbrock-type (39); (28); (29) methods were proved to be efficient. Recently, the exponential Rosenbrock method was extended to the case of stochastic partial differential equations (26) and was proved to be very stable for stochastic reactive dominated transport equations. Since solving linear systems are more straightforward than computing the exponential of a matrix, it is important to develop alternative methods based on the resolution of linear systems, which may be more efficient if the appropriate preconditionners are used. In this paper, we propose a new scheme based on the combination of the Rosenbrock-type method and the linear implicit Euler method. The resulting numerical scheme is called Stochatic Rosenbrock-type scheme(SROS). The space discretization is done with the finite element method and the new scheme is based on the local linearization of the nonlinear part of the semi-discrete problem in space. The local linearization therefore weakens the nonlinear part of the drift function such that the linearized semi-disrete problem is driven by its new linear part, which changes at each time step. The standard linear implicit Euler method (18); (44) is applied at the end to the linearized semi-discrete problem. This combination gives our new scheme SROS. We analyze the strong convergence of the new fully discrete scheme toward the exact solution in the -norm. The challenge here is that the resolvent of the operator appearing in the numerical scheme (41) is not constant as it changes at each time step. Furthermore the operator is a random operator. We use some tools from (26) and furthermore, we provide in Section 3.1 some stability estimates to handle the composition of the perturbed random resolvent operators, useful in our convergence analysis. The results indicate how the convergence orders depend on the regularity of the initial data and the noise. More precisely, we achieve the optimal convergence orders for multiplicative noise and the optimal convergence orders for additive noise, where is the regularity’s parameter of the noise (see Assumption 2.2) and is an arbitrary number small enough.

The rest of this paper is organized as follows. Section 2 deals with the well posedness problem, the numerical scheme and the main results. In section 3 we provide some error estimates for the deterministic homogeneous problem as preparatory results and then we provide the proof of the main results. Section 4 provides some numerical experiments to sustain the theoretical findings.

2 Mathematical setting and main results

2.1 Main assumptions and well posedness problem

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

(4)

where is an orthonormal basis of . Note that (4) is independent of the orthonormal basis of . For simplicity we use the notations and . It is well known that for all and , and

(5)

We assume that the covariance operator is positive and self-adjoint. Throughout this paper is a -wiener process. The space of Hilbert-Schmidt operators from to is denoted by with the corresponding norm defined by

(6)

where is an orthonormal basis of . Note that (6) is independent of the orthonormal basis of . In the rest of this paper, we take .

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

Assumption 2.1

[Linear operator ] is a generator of an analytic semigroup .

Assumption 2.2

[Initial value ] Let . We assume that , .

As in the current literature for deterministic Rosenbrock-type methods, (29); (28), deterministic exponential Rosemnbrock-type method (10); (27) and stochastic exponential Rosenbrock-type methods (26), we make the following assumption on the nonlinear drift term.

Assumption 2.3

[Nonlinear term ] We assume the nonlinear map to be Fréchet differentiable with bounded derivative, i.e. there exists a constant such that

(7)

Moreover, as in ((20), Page 6) for deterministic Rosenbrock-type method, we assume that the resolvent set of contains for all .

Remark 1

Inequality (7) together with the mean value theorem show that there exists a constant such that

(8)

As a consequence of (8), there exists a positive constant such that

Remark 2

An illustrative example for which the resolvent set of contains is obtained if we assume to be the generator of a contraction semigroup and the derivative of the nonlinear drift term to satisfy the following coercivity condition

(9)

In fact, it follows from (9) that is an relatively -bounded ((1), Chapter III, Definition 2.1) and dissipative operator with -bound . Therefore, from ((1), Chapter III, Theorem 2.7), it follows that is a generator of a contraction semigroup. Hence , for all .

Remark 3

The assumption that can be relaxed, but the drawback is that the resolvent set of the perturbed semigroup is smaller than that of the initial semigroup.

Following ((32), Chapter 7) or (45); (22); (13); (18) we make the following assumption on the diffusion term.

Assumption 2.4

[Diffusion term ] We assume that the operator satisfies the global Lipschitz condition, i.e. there exists a positive constant such that

(10)

As a consequence, it holds that

(11)

We equip , with the norm , for all . It is well known that is a Banach space (9).

To establish our strong convergence result when dealing with multiplicative noise, we will also need the following further assumption on the diffusion term when , which was also used in (13); (19) to achieve optimal regularity order and in (22); (18); (26) to achieve optimal convergence order in space and time.

Assumption 2.5

We assume that there exists a positive constant such that and for all , where comes from Assumption 2.2.

Typical examples which fulfill Assumption 2.5 are stochastic reaction diffusion equations (see ((13), Section 4)).

When dealing with additive noise (i.e. when ), the strong convergence proof will make use of the following assumption, also used in (44); (43); (26).

Assumption 2.6

We assume that the covariance operator satisfies the following estimate

(12)

where comes from Assumption 2.2.

When dealing with additive noise, to achieve higher convergence order in time, we assume that the nonlinear function satisfies the following assumption, also used in (43); (44); (26).

Assumption 2.7

The deterministic mapping is twice differentiable and there exist two constants and such that

(13)
(14)

Let us recall the following proposition which provides some smooth properties of the semigroup generated by , that will be useful in the rest of the paper.

Proposition 1

[Smoothing properties of the semigroup](9) Let , and , then there exists a constant such that

where , and .
If then .

Theorem 2.8

((32), Theorem 7.2)
Let Assumption 2.1, Assumption 2.3 and Assumption 2.4 be satisfied. If is a - measurable valued random variable, then there exists a unique mild solution of the problem (1) represented by (3) and satisfying

and for any , there exists a constant such that

2.2 Finite element discretization

In the rest of the paper, to simplify the presentation, we consider the linear operator to be of second-order. More precisely, we consider the SPDE (1) to be a second-order semilinear parabolic which takes the following form

(15)

where the functions and are continuously differentiable with globally bounded derivatives. In the abstract framework (1), the linear operator takes the form

(16)
(17)

where , . We assume that there exists a positive constant such that

The functions and are defined by

(18)

for all , , , with . For an appropriate family of eigenfunctions such that , it is well known ((13), Section 4) that the Nemystskii operator related to and the multiplication operator associated to defined in (18) satisfy Assumption 2.3, Assumption 2.4 and Assumption 2.5. As in (22); (4) we introduce two spaces and , such that ; the two spaces depend on the boundary conditions and the domain of the operator . 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 is the normal derivative of and is the exterior pointing normal to the boundary of , given by

Using the Green’s formula and the boundary conditions, the corresponding bilinear form associated to is given by

for Dirichlet and Neumann boundary conditions, and

for Robin boundary conditions. Using the Gårding’s inequality (see e.g. (36)), it holds that there exist two constants and such that

(19)

By adding and substracting in both sides of (1), we have a new linear operator still denoted by , and the corresponding bilinear form is also still denoted by . Therefore, the following coercivity property holds

(20)

Note that the expression of the nonlinear term has changed as we included the term in a new nonlinear term that we still denote by . The coercivity property (20) implies that is sectorial on , i.e. there exist such that

(21)

where (see (9)). The coercivity property (20) implies that is the infinitesimal generator of a contraction semigroup on such that

(22)

where denotes a path that surrounds the spectrum of . The coercivity property (20) also implies that is a positive operator and its fractional powers are well defined for any by

(23)

where is the Gamma function (see (9)). Let us now turn to the space discretization of our problem (1). We start by splitting the domain in finite triangles. Let be the triangulation with maximal length satisfying the usual regularity assumptions, and be the space of continuous functions that are piecewise linear over the triangulation . We consider the projection from to defined for every by

(24)

The discrete operator is defined by

(25)

Like , is also a generator of a bounded analytic semigroup . Let be a constant satisfying

(26)

As any semigroup and its generator, and satisfy the smoothing properties of Proposition 1 with a uniform constant (i.e. independent of ). Following (2); (21); (22); (4), we characterize the domain of the operator as follows:

The semi-discrete version of problem (1) consists to find , such that and

(27)

The following lemma can be found in ((26), Lemma 4 & Lemma 5). Its provides the space and time regularities of the mild solution of (38).

Lemma 1
  • For multiplicative noise, let Assumption 2.1 with , Assumption 2.2, Assumption 2.3 and Assumption 2.4 be fulfilled, then the following estimates hold

    (28)
    (29)

    for all .

    Moreover, if and if Assumption 2.5 is fulfilled, then

    (30)

    for all .

  • For additive noise, let Assumption 2.1, Assumption 2.2, Assumption 2.3 and Assumption 2.7 be fulfilled with , then the following estimates hold

    (31)
    (32)

    for all .

Here is a positive constant, independent of , , and .

Corollary 1

As a consequence of Lemma 1, it holds that

for all .

2.3 Standard linear implicit Euler method and stability properties

Let us recall that the linear implicit Euler scheme applied to the semi-discrete problem (27) leads to

(33)
(34)

If the linear operator tends to the null5 operator, its corresponding discrete operator tends to the null operator and tends to the identity operator . In this case, the numerical scheme (33) and the standard exponential integrator method (22) become the unstable Euler-Maruyama scheme. See also ((26), Section 2.3) for more details. For a simple illustration of the stability properties of such problems, let us consider the following deterministic linear differential equation

(35)

The linear implicit Euler method applied to (35) by considering as the nonlinear part reads as

(36)

The numerical scheme (36) is stable (39); (30) if and only if . Note that when is small enough and large enough, the numerical scheme (36) becomes the explicit Euler method and the stability region becomes very small. The Rosenbrock-type methods were proved to be efficient and were studied in (8); (7); (30) for ordinary differential equations. Applying the Rosenbrock-Euler method to the linear problem (35) reads as

(37)

Note that (37) coincides with the full implicit method with . The Rosenbrock-Euler method (37) is unconditionally stable (A-stable). This demonstrates the strong stability property of the Rosenbrock-type methods for stiff problems. Authors of (28); (29) extended the Rosenbrock-type methods to parabolic partial differential equations and the methods were proved to be efficient for solving transport equations in porous media (39). The case of stiff stochastic partial differential equations is not yet studied to the best of our knowledge and will be the aim of this paper.

2.4 Novel fully discrete scheme and main results

Let us build a more stable scheme, robust when the operator tends to null. For the time discretization, we consider the one-step method which provides the numerical approximated solution of at discrete time , . The method is based on the continuous linearization of (27). More precisely we linearize (27) at each time step as follows

(38)

for all , where is the Fréchet derivative of at and is the remainder at . Both and are random functions and are defined for all by

(39)
(40)

Applying the linear implict Euler method to (38) gives the following fully discrete scheme, called linear implicit Rosenbrock Euler method

(41)

where and are defined respectively by

(42)

and

(43)

In the numerical scheme (41), the resolvent operator in (42) is random and changes at each time step. Having the numerical method (41) in hand, our goal is to analyze its strong convergence toward the exact solution in the norm for multiplicative and additive noise.

Throughout this paper we take , where for , , is fixed, is a generic constant that may change from one place to another but is independent of both and . The main results of this paper are formulated in the following theorems.

Theorem 2.9

Let and be respectively the mild solution given by (3) and the numerical approximation given by (41) at . Let Assumption 2.1, Assumption 2.2 with , Assumption 2.3 and Assumption 2.4 be fulfilled.

  • If , then the following error estimate holds

  • If and if Assumption 2.5 is fulfilled, then the following error estimate holds

Theorem 2.10

When dealing with additive noise (i.e. when ), let Assumption 2.1, Assumption 2.2 with , Assumption 2.3, Assumption 2.6 and Assumption 2.7 be fulfilled. Then the following error estimate holds for the mild solution of (1) and the numerical approximation (41)

(44)

3 Proof of the main results

The proofs of the main results require some preparatory results.

3.1 Preparatory results

For non commutative operators on a Banach space, we introduce the following notation, which will be used in the rest of the paper.

(45)
Lemma 2

Let Assumption 2.2 be fulfilled. Then for all the following estimate holds

(46)
Proof

See ((26), Lemma 10).

Lemma 3

For all and all , the random linear operator is the generator of an analytic semigroup called random (or stochastic) perturbed semigroup and uniformly bounded on , i.e there exists a positive constant independent of , , and the sample such that

Proof

See ((26), Lemma 5).

The following lemma is an analogous of ((24), (3.31)), but here our semigroup is not constant. In fact it is random and further its changes at each time step.

Lemma 4

Let Assumption 2.1 and Assumption 2.3 be fulfilled.

  • For all , and all , it holds that

    (47)
  • For all and , it holds that

    (48)
  • For all , it holds that

    (49)
Proof

Note that for all , . Therefore . It remains to prove the first part of the estimate (47). Using the interpolation theory, we only need to prove that (47) holds for and . Since and the resolvent set of contains , it follows from ((31), (5.23)) that

(50)

Taking the norm in both sides of (50) and using the uniformly boundedness of (see Lemma 3) yields

(51)

Using the change of variable yields

(52)

This shows that (47) holds for . Pre-multiplying both sides of (50) by yields