## Abstract

We consider the numerical approximation of parabolic stochastic partial differential equations driven by additive space–time white noise. We introduce a new numerical scheme for the time discretization of the finite-dimensional Galerkin stochastic differential equations, which we call the exponential Euler scheme, and show that it converges (in the strong sense) faster than the classical numerical schemes, such as the linear-implicit Euler scheme or the Crank–Nicholson scheme, for this equation with the general noise. In particular, we prove that our scheme applied to a semilinear stochastic heat equation converges with an overall computational order 1/3 which exceeds the barrier order 1/6 for numerical schemes using only basic increments of the noise process reported previously. By contrast, our scheme takes advantage of the smoothening effect of the Laplace operator and of a linear functional of the noise and, therefore overcomes this order barrier.

## 1. Introduction

The numerical approximation of stochastic partial differential equations (SPDEs) encounters all the difficulties that arise in the numerical solution of both deterministic partial differential equations (PDEs) and finite-dimensional stochastic differential equations (SDEs) plus more due to the infinite dimensional nature of the driving noise processes. See, for example, Gyöngy & Nualart (1995, 1997), Grecksch & Kloeden (1996), Gyöngy (1998, 1999), Shardlow (1999), Davie & Gaines (2000), Kloeden & Shott (2001), Hausenblas (2003), Lord & Rougemont (2004) and Müller-Gronbach & Ritter (2007).

In this paper, we consider the numerical approximation of a parabolic SPDE with additive noise, i.e. of the form(1.1)where *A* is in general an unbounded operator (e.g. *A*=*Δ*), *f* is a nonlinear continuous function and *W*_{t} is, for example, a space–time white noise (see §2 for a precise description of equation (1.1)).

Various approximation schemes have been introduced in the literature for this problem (see the papers cited above). All have the property that their rate of convergence for the one-dimensional semilinear stochastic heat equation with additive space–time white noise is 1/6. Indeed, Davie & Gaines (2000; see also Müller-Gronbach & Ritter 2007) showed that no numerical schemes using only equidistant evaluations of the noise can converge with an overall rate faster than 1/6 (with strong convergence in time).

Here, we present a new scheme that overcomes the rate barrier of 1/6 and converges with the rate of 1/3. It uses linear functionals of the noise, which are Gaussian distributed and are thus not difficult to generate. Of course, the choice of a semilinear heat equation is just for illustration and the situation is analogous for SPDE with the more general structure (1.1). The smoothening effect of the semigroup generated by the operator *A* in the SPDE (1.1) and various semigroup estimates play an important role in the proof.

In §2, we give a precise description of the equations under consideration and the assumptions that we need. Then, in §3, we will introduce our numerical scheme and state the convergence theorem, while in §4 we present an example, briefly review the above-cited literature and give some numerical results. Proofs are given in §5.

## 2. Mathematical setting and assumptions

Let *T*>0 and let (*Ω*, , ) be a probability space with a normal filtration , *t*∈[0,*T*]. In addition, let (*H*, 〈.,.〉) be a separable Hilbert space with norm denoted by |.|.

We will interpret the SPDE (1.1) in such a space *H* in the mild sense, i.e. as satisfying the integral equation(2.1)for all *t*∈[0,*T*]. See, for example, Da Prato & Zabczyk (1992) for details. The objects *A*, *u*_{0}, *f* and *W*_{t} here are specified through the following assumptions.

There exist sequences of real eigenvalues 0<*λ*_{1}≤*λ*_{2}≤⋯ and eigenfunctions {*e*_{n}}_{n≥1} of *A* such that the linear operator *A* : *D*(*A*)⊂*H*→*H* is given byfor all *v*∈*D*(*A*) with .

Furthermore, let *D*((−*A*)^{r}) with *r*∈ denote the interpolation spaces of the operator −*A* (e.g. Sell & You 2002).

There exist a sequence *q*_{n}≥0, *n*≥1, of positive real numbers, a real number *γ*∈(0,1) such thatand independent real valued -Brownian motions , *t*≥0, for *n*≥1, i.e. each is -adapted and the increments , *h*>0, are independent of . Then, the cylindrical Brownian motion *W*_{t} is given by

The above series may not converge in *H*, but in some space *U*_{1} into which *H* can be embedded, see Da Prato & Zabczyk (1992) and Prévot & Röckner (2007). In our example with the Laplace operator in one dimension, we will have *λ*_{n}=−*π*^{2}*n*^{2} and *q*_{n}≡1 for *n*≥1. This is the important case of space–time white noise.

The nonlinearity *f* : *H*→*H* is two times continuously Fréchet differentiable, it and its derivatives satisfyfor all *x*, *y*∈*H*, *v*∈*D*((−*A*)^{r}) and *r*=0, 1/2, 1 and they satisfyfor all *v*, *w*, *x*∈*H*, where *L*>0 is a positive constant.

The function *f* is usually given as a real-valued function of a real variable, but in the SPDE (1.1) it is considered as a function defined on *H* and taking values in some function space such as a subspace of *H*.

The initial value *u*_{0} is a *D*((−*A*)^{γ}) valued random variable, which satisfieswhere *γ*>0 is given in assumption 2.2.

Later, for the definition of the numerical scheme, we will also use the particular knowledge of the eigen decomposition of the linear operator *A* given in assumption 2.1. Under assumption 2.2, it is well known that the stochastic convolution for *t*∈[0,*T*] has a predictable modification with values in *D*((−*A*)^{γ}). Hence, under assumptions 2.1, 2.2, 2.4 and 2.6 it is well known (see Da Prato & Zabczyk 1992) that the SPDE (1.1) has an up to modifications unique mild solution *U*_{t} on the time interval [0,*T*], where *U*_{t} is the predictable stochastic process with values in *D*((−*A*)^{γ}), which fulfils (2.1). The following well-known results will be needed later:(2.2)and(2.3)for all *t*, *s* ∈[0,*T*] and a constant *C*>0 (see Jentzen 2008).

## 3. Main theorem

We will now apply a Galerkin projection to the SPDE (1.1) and then introduce a numerical scheme, which we will call the exponential Euler scheme, to the finite-dimensional Galerkin SDE and, finally, we will state a convergence theorem for this scheme.

We define finite-dimensional subspaces *H*_{N} of *H* byand the projections *P*_{N} : *H*→*H*_{N} byfor *N*∈. Then, we truncate or project the initial random variable *u*_{0}, the nonlinearity *f* and the Brownian motion *W*_{t} byfor all *v*∈*H* and each *N*≥1.

Using this notation, we introduce a finite-dimensional SDE in the space *H*_{N} (or, equivalently, in ) by(3.1)which is the Galerkin projection of the SPDE (1.1) onto *H*_{N}, where *A*_{N}=*P*_{N}*A* is the (matrix) operator *A*_{N} : *H*_{N}→*H*_{H}. Since assumption 2.4 also applies to *f*_{N}, the SDE (3.1) also has a unique solution on [0,*T*], which is given (implicitly) by (see Kloeden & Platen 1992)(3.2)

Now, we introduce a numerical method to approximate in time on the interval [0,*T*], which we will call the *exponential Euler scheme*: let and define(3.3)with time-step *h*=*T*/*M* for some *M*∈ and discretization times *t*_{k}=*kh* for *k*=0, 1, …, *M*. Note that *A*_{N} is a diagonal matrix with diagonal entries *λ*_{1}, …, *λ*_{N}. Therefore, we can easily compute the expressions and directly, which can also be seen below. Of course, we need to know the eigen decomposition of the linear operator *A* in order for the method to be applicable (see assumption 2.1). Then, this scheme is easier to simulate than may seem on first sight. Denoting the components of and *f*_{N} bywe can rewrite the numerical scheme (3.3) aswhere the for *i*=1, …, *N* and *k*=0, 1, …, *M*−1 are independent, standard normally distributed random variables.

With this notation, we are able to present our main theorem, which states the strong convergence of the above numerical scheme and also provides a rate for this strong convergence.

*Suppose that* *assumptions 2.1*, *2.2*, *2.4* *and* *2.6* *are satisfied*. *Then*, *there is a constant C*_{T}>0 *such that*(3.4)*holds for all N*, *M*≥2, *where U*_{t} *is the solution of SDE* (*1.1*), *is the numerical solution given by* (*3.3*), for *k*=0, 1, …, *M and γ*>0 *is the constant given in* *assumption 2.2*.

It follows immediately that the exponential Euler scheme (3.3) converges in time with strong order 1−*ϵ* for an arbitrary small *ϵ*>0, since log(*M*) can be estimated by *M*^{ϵ}. Exponential integrators for SDEs have already been studied intensively in the literature, see, for example, Jimenez *et al*. (1999) or also Carbonell *et al*. (2005). The schemes there applied to the finite-dimensional Galerkin SDE (3.1) are similar to the scheme proposed in equation (3.3). However, in Jimenez *et al*. (1999; see eqns (12) and (13) there), they only used evaluations of the noise process instead of linear functionals and so this would not lead to a higher order than the linear-implicit Euler scheme (see §4*c*). Furthermore, in Carbonell *et al*. (2005), they used linear functionals of the noise but they compute the derivative , which is an *N*×*N* matrix, in each computation step, which is very inefficient with respect to the computational cost. They considered a general SDE with additive noise and did not take into account the special structure coming from the SPDE (1.1). In the theorem above, we estimated the mean square error, which is in a way standard for strong convergence (see theorem 10.6.3 in Kloeden & Platen 1992). Alternatively, one could estimate the *p*th moment for every *p*≥1 and so also obtain pathwise convergence (see Jentzen *et al*. accepted), which is a task for further research. For it, one needs a Burkholder–Davis–Gundy inequality in Hilbert spaces, which is a non-trivial tool. Finally, note that for deterministic parabolic PDEs, exponential integrators, i.e. numerical schemes with exponential terms, have for example been considered in Hochbruck & Ostermann (2005*a*,*b*).

## 4. An example

To present more clearly the consequences of theorem 3.1, we consider a much studied example of the SPDE (1.1). After introducing the example and applying theorem 3.1 to it we will review results in the literature concerning the numerical approximation of this example and, more generally, for the SPDE (1.1). Finally, we will illustrate our results with some numerical calculations.

### (a) The semilinear stochastic heat equation

As a simple example of the SPDE (1.1), we consider a semilinear stochastic heat equation with additive space–time white noise on the one-dimensional domain [0,1] over the time interval [0,*T*] with *T*=1, i.e.(4.1)with the Dirichlet boundary condition and the initial value *u*(*x*,0)=*u*_{0}(*x*) for *x*∈(0,1).

Here, *H*=*L*^{2}[0,1] and the linear operator *A* is the one-dimensional Laplace operator with the Dirichlet boundary condition *A*=*Δ*≔*D*(*Δ*)→*H*, wheregiven byfor all *u*∈*D*(*Δ*), where *u*″ denotes the second weak derivative of the function *u*.

The eigenfunctions *e*_{n} in *H* and eigenvalues −*λ*_{n} of *Δ* are given byfor all *n*≥1. Hencewithfor all *u*∈*D*(*Δ*) with .

The noise *W*_{t} here is the space–time Wiener processwith *q*_{n}≡1 for all *n*≥1 in view of assumption 2.2. (The summation here is just formal. In particular, it does not converge in *H*.) Therefore, we have with an arbitrary small *ϵ*>0 in our situation.

### (b) Consequences of theorem 3.1 applied to the SPDE (4.1)

From theorem 3.1, the exponential Euler scheme (3.1) applied to the SPDE (4.1) has the error estimate(4.2)Since *M*.*N* arithmetical operations, random number and function evaluations are needed to calculate , the computational cost of the scheme (3.3) is *K*=*N*.*M*. In view of the above error bound, it is optimal to choose and and we have the optimal error bound(4.3)Hence, the convergence rate of the scheme with respect to the computational cost is (1/3)−*ϵ* for an arbitrary small *ϵ*>0.

Note that, in general, if a numerical solution for some scheme with computational cost *N*.*M* applied to this SPDE converges in the sense(4.4)for *α*, *β*>0, then it converges with a rate of with respect to the computational cost. In our situation, we have and *β*=1−*ϵ*, so we obtain the overall rate (1/3)−*ϵ*.

### (c) Articles in the literature on the SPDE (4.1)

We briefly review some articles in the literature concerning the numerical approximation of the SPDE (4.1) and the more general SPDE (1.1).

Gyöngy & Nualart (1995) introduced an implicit numerical scheme for the SPDE (4.1) and showed that it converges strongly to the exact solution without giving a rate (see their theorem 2.1).

Grecksch & Kloeden (1996) considered a parabolic SPDE driven by a (possibly multiplicative) scalar Wiener process and applied explicit Itô–Taylor schemes (see Kloeden & Platen 1992) to the Galerkin SDEs, while Kloeden & Shott (2001) applied linear-implicit Itô–Taylor schemes to the Galerkin SDEs. In both cases, a higher order of strong convergence was obtained, but this was due to the special nature of the noise. In the general situation with space–time white noise, these schemes only converge with the order of 1/6.

Gyöngy & Nualart (1997) considered an SPDE with multiplicative space–time white noise and applied a temporal numerical scheme, which is based on the simulation of random variables of the formwhere *σ* is the diffusion coefficient and *X*_{k} is the numerical scheme at time *t*_{k}. They showed that this scheme converges with a rate of 1/8 in time (see their theorem 3.1).

Shardlow (1999) applied finite differences to the SPDE (4.1) to obtain a spatial discretization that he then discretized in time with a *θ*-method. This had an overall convergence rate 1/6 with respect to the computational cost (see his theorem 3.7).

Gyöngy (1998, 1999) also applied finite differences for an SPDE driven by space–time white noise and then used several temporal implicit and explicit schemes, in particular, the linear-implicit Euler scheme. He showed that these schemes converge with order 1/2 in space and with order 1/4 in time (assuming a smooth initial value). Hence, he obtained an overall convergence rate of 1/6 with respect to the computational cost in space and time (see theorem 3.1(iii) in Gyöngy 1999).

In a seminal paper, Davie & Gaines (2000) showed that any numerical scheme applied to the SPDE (4.1) with *f*=0 which uses only equidistant values of the noise *W*_{t} cannot converge faster than the rate of 1/6 with respect to the computational cost (see §2.1 in Davie & Gaines 2000). Moreover, for the special situation with *f*=0, they simulated the solution exactly (see also Müller-Gronbach *et al*. 2007). Müller-Gronbach & Ritter (2007) also showed that this is a lower bound for the convergence rate for the SPDE above, but with multiplicative noise. (They even showed that in one dimension (see equation (4.1)), one cannot improve this rate of convergence by choosing non-uniform evaluations of the noise.) However, Davie & Gaines (2000, p. 129) remarked that it may be possible to improve the convergence rate by using suitable linear functionals of the noise, which is essentially what we have done in this paper.

Hausenblas (2003) applied the linear-implicit and explicit Euler scheme and the Crank–Nicholson scheme to an SPDE (4.1) driven by an infinite dimensional noise. In the case of a smoother noise, i.e. trace-class noise, she obtained the order 1/4 with respect to the computational cost. However, in the general case of space–time white noise, the convergence rate is no better than 1/6.

Finally, Lord & Rougemont (2004) also considered the SPDE (4.1) with a smoother noise. They discretized the Galerkin SDE in time with a numerical scheme that also uses the factor with *A*_{N}=*P*_{N}*A* for *A*=*Δ*. The *Lord*–*Rougemont scheme* is given by(4.5)where is the numerical solution at time *t*_{k}. They showed that this scheme is useful when the noise is very smooth in space, in particular, with Gevrey regularity. However, in the general case of space–time white noise, the scheme (4.5) converges also only with the rate of 1/6 with respect to the computational cost, which is a consequence of the work of Davie & Gaines (2000) and can also be seen in the numerical simulations below.

### (d) Numerical results

We illustrate theorem 3.1 with a numerical example. In particular, we consider the semilinear stochastic heat equation (4.1) on the one-dimensional domain (0,1) with , i.e.(4.6)with the Dirichlet boundary condition and the initial value , *T*=1. Of course, in general, *f* will be nonlinear, but we choose *f* to be linear here since we have an exact solution for comparison with the numerical solution.

We applied the linear-implicit Euler scheme (with just the Laplacian part implicit) to the SPDE (4.6), which is given byand the Lord–Rougemont scheme (4.5) as well as the exponential Euler scheme (3.3). The spatial error of the Galerkin projection, i.e. the difference of equation (2.1) and (3.2), is of the order of 1/2, so we first focus on the convergence rate in the time of the numerical schemes above. We fix *N*=100 (space discretization) and then apply the above schemes with different *M*=10, 20, 25, 40, 50, 80, 100, 200, 500 and 1000 (time discretization). Figure 1 shows the error as log–log plot, where the expectation in the error is approximated by the mean of 100 independent realizations.

The linear-implicit Euler and Lord–Rougemont schemes clearly converge with the rate 1/4 in time, while the exponential Euler scheme converges with the temporal rate 1. Now, we focus on their convergence rate with respect to their computational cost. Here, figure 2 shows the error as log–log plot, where the expectation in the error is approximated again by the mean of 100 independent realizations.

Here, the linear-implicit Euler and Lord–Rougemont schemes clearly converge with the rate 1/6, while the exponential Euler scheme converges with the rate 1/3. All the three schemes thus converge with their theoretically predicted order. Finally in figure 3, we present the approximation of the solution of equation (4.6) computed by the numerical scheme introduced in that paper at the four different times *t*=0, 1/1000, 1/100 and *t*=*T*=1.

## 5. Proof of theorem 3.1

In what follows, *C* is a constant which changes from line to line. Let *N*, *M*∈ be arbitrary with *M*≥2. We divide the proof into two parts. Our main goal is to show(5.1)where denotes the *L*^{2}(*Ω*)-norm, given byfor an *H*-valued random variable *X*. First, we will show in §5*a* that the spatial error satisfies(5.2)and then, in §5*b*, that the temporal error satisfies(5.3)

Combining both results will then give the desired estimate (5.1).

### (a) Spatial error

Before we begin, we note the inequalityfor all *v*∈*D*((−*A*)^{γ}), where |*v*|_{γ}=|(−*A*)^{γ}*v*|. We havefrom which we obtainwhich impliesdue to assumption 2.4. From this, it follows that:where we have used the Hölder inequality. The Gronwall lemma then yieldsand thereforedue to equation (2.2). This proves the estimate (5.2).

### (b) Time discretization

Now let *k* ∈ {0, 1, …, *M*−1} be arbitrary. We haveand

Thus, we need to estimate(5.4)whereOf course, we haveWe will now show thatwhich by inequality (5.4), implies thatThus, by Gronwall's Lemma, we havewhich implies the estimate (5.3).

Hence, it remains to estimate the quantity . For this, we use the Taylor ExpansionwithWe obtainwithIn §5*b*(i)–(iv), we show that for *i*=1, …, 4, from which will finally follow that:and this implies inequality (5.3).

#### (i) First term

First, we consider the *E*_{1} term.To begin, we havedue to assumption 2.4. Hencedue to assumption 2.4. Finally, we obtainwhich is the claim for *E*_{1}.

In the literature, the term *E*_{1} is usually estimated in the following way:which yieldsIn this way, one can only obtain a convergence rate of *γ* in time, which for our example, would be , *ϵ*>0 small. To obtain a higher order, we need to use the smoothening effect of the term as we did earlier, which is based on the estimateThe situation is similar for the estimates of the terms *E*_{3} and *E*_{4}.

#### (ii) Second term

We need the estimate(5.5)for all *v*∈*H*, for which this is true, since by assumption 2.4,for all *v*∈*H*. Thus for *E*_{2}, we obtainTherefore, we havewhich implies

#### (iii) Third term

Here, we havedue to assumption 2.4 and inequality (2.3). Therefore, we obtaindue to assumption 2.4 again. Finally, we have

#### (iv) Fourth term

Note thatHence, assumption 2.4 implies that(5.6)Now, we estimate *E*_{4}. Sinceand sinceholds for *l*=0, 1,…, *k*−1 due to the inequalities (2.2) and (2.3), we obtainin view of inequality (5.6). ▪

## Acknowledgments

We thank the referees for their valuable comments. This work was supported by the DFG project ‘Pathwise numerical analysis of stochastic evolution equations’.

## Footnotes

- Received August 7, 2008.
- Accepted October 20, 2008.

- © 2008 The Royal Society