Stochastic Exponential Integrators for a Finite Element Discretization of SPDEs

Stochastic Exponential Integrators for a Finite Element Discretization of SPDEs


We consider the numerical approximation of general semilinear parabolic stochastic partial differential equations (SPDEs) driven by additive space-time noise. In contrast to the standard time stepping methods which uses basic increments of the noise and the approximation of the exponential function by a rational fraction, we introduce a new scheme, designed for finite elements, finite volumes or finite differences space discretization, similar to the schemes in [4] for spectral methods and [1] for finite element methods. We use the projection operator, the smoothing effect of the positive definite self-adjoint operator and linear functionals of the noise in Fourier space to obtain higher order approximations. We consider noise that is white in time and either in or in space and give convergence proofs in the mean square norm for a diffusion reaction equation and in mean square norm in the presence of an advection term. For the exponential integrator we rely on computing the exponential of a non-diagonal matrix. In our numerical results we use two different efficient techniques: the real fast Léja points and Krylov subspace techniques. We present results for a linear reaction diffusion equation in two dimensions as well as a nonlinear example of two-dimensional stochastic advection diffusion reaction equation motivated from realistic porous media flow.


We consider the strong numerical approximation of Ito stochastic partial differential equations

in a Hilbert space , where and , and initial data is given. The linear operator is the generator of an analytic semigroup with eigenfunctions and is a nonlinear function. The noise can be represented as a series in the eigenfunctions of the covariance operator and we assume for convenience that and have the same eigenfunctions (without loss the generality as the orthogonal projection can be used). In which case [24] we have

where , are the eigenvalues of the covariance operator . The are independent and identically distributed standard Brownian motions.

The study of numerical solutions of SPDEs is an active research area and there is an extensive literature on numerical methods for SPDEs. Recent work by Jentzen and co-workers [2] uses the Taylor expansion and linear functionals of the noise for Fourier–Galerkin discretizations of . In these schemes the diagonalization of the operator through the discretization plays a key role. Using a linear functional of the noise overcomes the order barrier encountered using a standard increment of Wiener process [4]. In [1] the use of linear functionals of the noise is extended to finite–element discretizations (where the operator does not diagonalize) with a semi-implicit Euler–Maruyama method. In contrast to [1] here we consider two exponential based methods for time-stepping as in [21]. We prove a strong convergence result for two versions of the scheme with noise that is white in time and in and in space that shows that the exponential integrators are more accurate that the semi-implicit Euler-Maruyama method. Furthermore we have weaker restrictions on the regularity of the initial data and high accuracy for linear problems comparing to the scheme in [1]. The cost of the extra accuracy though is that to implement these methods we need to compute the exponential of a non–diagonal matrix, which is a notorious problem in numerical analysis [35]. However, new developments for both Léja points and Krylov subspace techniques [33] have led to efficient methods for computing matrix exponentials. Compared to the Fourier-Galerkin methods of [2] we gain the flexibility of finite element (or finite volume methods) to deal with complex boundary conditions and we can apply well developed techniques such as upwinding to deal with advection.

We consider two examples of where is the second order operator and is the diffusion coefficient. For the first example

where is a constant, we can construct an exact solution. Our second example is a stochastic advection diffusion reaction equation in a heterogeneous porous media

where is the Darcy velocity field [22]. In the linear example we take Neumann boundary conditions and for the example from porous media we take a mixed Neumann–Dirichlet boundary condition.

The paper is organised as follows. In Section we present the two numerical schemes based on the exponential integrators and our assumptions on . We present and comment on our convergence results. In Section we present the proofs of our convergence theorems. We conclude in Section by presenting some simulations and discuss implementation of these methods.

2Numerical scheme and main results

We start by introducing our notation. We denote by the norm associated to the standard inner product of the Hilbert space and the norm of the Sobolev space , for . For a Banach space we denote by the set of bounded linear mapping from to and the set of bounded bilinear mapping from to . We introduce further spaces and notation below as required.

Consider the stochastic partial differential equation (Equation 1), under some technical assumptions it is well known (see [24] and references therein) that the unique mild solution is given by

with the stochastic process given by the stochastic convolution

We consider discretization of the spatial domain by a finite element triangulation. Let be a set of disjoint intervals of (for ), a triangulation of (for ) or a set of tetrahedra (for ). Let denote the space of continuous functions that are piecewise linear over . To discretize in space we introduce two projections. Our first projection operator is the projection onto defined for by

We can then define the operator , the discrete analogue of , by

We denote by the semigroup generated by the operator . The second projection , is the projection onto a finite number of spectral modes defined for by

Furthermore we can project the operator

We discretize in space using finite elements and project the noise first onto a finite number of modes and then onto the finite element space. The semi-discretized version of is to find the process such that

The mild solution of (Equation 4) at time , is given by

Given the mild solution at the time , we can construct the corresponding solution at as

For our first numerical scheme (SETD1), we use the following approximations


Then we approximate of by


For efficiency to avoid computing two matrix exponentials we can rewrite the scheme (Equation 7) as

We call this scheme (SETD1).

Our second numerical method (SETD0) is similar to the one in [21]. It is based on approximating the deterministic integral in at the left–hand endpoint of each partition and the stochastic integral as follows

With this we can define the (SETD0) approximation of by


If we project the eigenfunctions of onto the eigenfunctions of the linear operator then by a Fourier spectral method the process

is reduced to an Ornstein–Uhlenbeck process in each Fourier mode as in [4] and we therefore know the exact variance in each mode. We comment further on the implementation in Section . We describe now in detail the assumptions that we make on the linear operator , on our finite element discretization, the nonlinear term and the noise .

Note that for convenience of presentation we take to be a second order operator as this simplifies notation for the norm equivalence below. Similar result hold, however, for higher order operators. We recall some basic properties of the semigroup generated by that may be found for example in [13].

We introduce two spaces and where that depend on the choice of the boundary conditions. For Dirichlet boundary conditions we let

and for Robin boundary conditions, for which Neumann boundary conditions are a particular case, and

see [16] for details. Functions in satisfy the boundary conditions and with in hand we can characterize the domain of the operator and have the following norm equivalence [14] for

We now introduce our assumptions on the spatial domain and finite element space . We consider the space of continuous functions that are piecewise linear over the triangulation with .

This inequality is sometimes called the Bramble and Hilbert inequality, see [16] or [15]. It follows that

We assume that the function is defined in , although in general may be defined in any Hilbert space. The possible choice of can be , or the Banach space of continuous functions from to denoted by if .

We now turn our attention to the noise term and introduce spaces and notation that we need to define the -Wiener process. Denoting by the Banach algebra of bounded linear operators on with the usual norm. We recall that an operator is Hilbert-Schmidt if

If we denote the space of Hilbert-Schmidt operator from to by i.e

the corresponding norm by

Let be a process such that for every sample , . Then we have the following equality

using Ito’s isometry [24]. We assume sufficient regularity of the noise for the existence of a mild solution and to project the noise into the finite element space . To be specific we assume the noise is in either in or in space.

Using the equivalence of norms, we have that

2.1Main results

Throughout the paper we let be the number of terms of truncated noise and and take , where for . We take to be a constant that may depend on and other parameters but not on , or . We also assume that when initial data then , and .

Our first result is a strong convergence result in when the non-linearity satisfies the Lipschitz condition of Assumption (a) with scheme (SETD1). This is, for example, the case of reaction–diffusion SPDEs.

Our first result for scheme (SETD0) is a strong convergence result in when the non-linearity satisfies the Lipschitz condition of Assumption (a).

For convergence in the mean square norm where the non-linearity satisfies the Lipschitz condition from norm to ( Assumption (b)) we can state results for (SETD1) and (SETD0) together.

We note that this theorem covers the case of advection-diffusion-reaction SPDEs, such as that arising in our example from porous media.

We remark that if we denote by the numbers of vertices in the finite elements mesh then it is well known (see for example [27]) that if \(N \geq N_{h}\\) then

As a consequence the estimates in Theorem ?, Theorem ? and Theorem ? can be expressed in function of and only and it is the error from the finite element approximation that dominates. If then it is the error from the projection of the noise onto a finite number of modes that dominates.

From Theorem we also get an estimate in the root mean square norm in the case that the nonlinear function satisfies Assumption (b). We cannot do the proof directly in due to the Lipschitz condition in Assumption (b). Simulations for Theorem will be do in since the discrete norm is more easy to use in all type of boundary conditions.

Finally if we compare these theorems to those in [1] for a modified semi-implicit Euler-Maruyama method then we see that using the exponential based integrators we have weaker conditions on the initial data and in particular the scheme (SETD1) has better convergence properties.

3Proofs of main results

3.1Preparatory results

We start by examining the deterministic linear problem. Find such that such that

The corresponding semi-discretization in space is : find such that

where . Define the operator

so that .

We now consider the SPDE

The first claim of part (i) of the Lemma can be found in [1] and so we prove the second part of (i). Consider the difference

so that

We estimate each of the terms . For , using Proposition yields

Then . For the term , we have

We now estimate each term and . For boundedness of gives

For we have

Using the Lipschitz condition in Assumption (a) with the first claim of (i) yields

Assumption (a) gives

Using the fact that is bounded we find

Combining the previous estimates ends the proof of the second claim of (i).

We now prove part (ii) of the lemma. Consider the difference

and then

Let us estimate the terms and and we start with . If using Proposition yields


For the term , we have

We now estimate each term above. Using the fact that in we have the equivalency of norm , we have

where is the norm of the semigroup viewed as a bounded operator in . We also have the similar relationship for the operator with , .

Using similar inequality as (Equation 12) yields

For small enough, we have

We also have using Proposition

Hence, if with , we have

We also have for the term

The Ito isometry property yields

Using Proposition , the fact that is bounded and yields

with small enough. Let us estimate . The fact that yields