Weak backward error analysis for SDEs

Arnaud Debussche, Erwan Faou

Introduction

In the last decades, backward error analysis has become one of the most powerful tool to analyze the long time behavior of numerical schemes applied to evolution equations. The main idea can be described as follows: Let us consider an ordinary differential equation of the form

The idea of backward error analysis is to show that $\Phi_{\tau}$ can be interpreted as the exact flow $\varphi_{\tau}^{f_{\tau}}$ or a modified vector field defined as a series in powers of $\tau$

where $f_{\tau}^{N}$ is the truncated series:

Under some analyticity assumptions, the constant $C_{N}\tau^{N}$ can be optimized in $N$ , so that the error term in the previous equation can be made exponentially small with respect to $\tau$ .

Such a result is very important and has many applications in the case where $f$ has some strong geometric properties, such as Hamiltonian or reversible structure. In this situation, and under some compatibility conditions on the numerical method $\Phi_{\tau}$ , the modified vector field $f_{\tau}$ inherits the structure of $f$ . For example if $\Phi_{\tau}$ is symplectic and $f$ Hamiltonian, then $f_{\tau}$ remains Hamiltonian. This has major consequences such as the preservation of a modified Hamiltonian over very long time (of order $\tau^{-N}$ ) for the numerical solution, from which we can deduce long time stability results, existence of numerical invariant tori in the integrable case, etc…

In the Hamiltonian case, this idea goes back to Moser , but was applied later to symplectic integrator by Benettin & Giorgilli , Hairer & Lubich and Reich . Such results now form the core of the modern geometric numerical integration theory for which we refer to the classical textbooks and .

More recently, these ideas have been extended in some situations to Hamiltonian PDEs: First in the linear case , and then in the semilinear case (nonlinear Schrödinger or wave equations), see .

An attempt has been made by Shardlow to extend the backward error analysis to this context. He has shown that the construction of a modified SDE associated with the Euler scheme can be performed, but only at the first step, ie for $N=2$ , and only for additive noise, ie when $\sigma(X)$ does not depend on $X$ . In this case, he is able to write down a modified SDE:

He explains that for multiplicative noise or higher order, there are too many conditions to be satisfied by the coefficients of the modified equations.

where $L$ is the order $2$ Kolmogorov operator associated with the SDE.

In the case of the Euler method applied to a SDE, we show that with the numerical solution, we can associate modified Kolmogorov operator of the form

Note that in contrast with the classical case, we do not have a modified SDE and cannot straightforwardly define a solution to the modified equation

However, in the case where the SDE is elliptic or hypoelliptic, we can build an approximated solution $v^{(N)}$ such that

Furthermore, using the exponential convergence to equilibrium, we prove that in fact the constant $c_{4}$ does not depend on $T$ so that we have an approximation result valid on very long times. We also show that there exist a modified invariant measure for $L^{(N)}(\tau,x,\partial_{x})$ .

We can then use this weak backward error analysis to prove that the numerical solution $X_{p}$ , $p\geq 1$ obtained by the Euler scheme is exponentially mixing up to some very small error, and for all times. This is typically a geometric numerical integration result in the sense that we prove the persistence of a qualitative property of the exact flow (exponential mixing) to the numerical approximation, over long times.

Note that $v^{(N)}$ is in fact constructed as a truncated series: $v^{(N)}=\sum_{n=0}^{N}\tau^{n}v_{n}$ and that $v_{0}=u$ is the solution of the Kolmogorov equation. Therefore, our result provides an expansion of the error as in (see also ). However, the expansion is different here.

Error estimates on long times for elliptic and hypoelliptic SDEs have already been proved. In , it is shown that for a sufficiently small time step the Euler scheme defines an ergodic process and that the invariant measure of the Euler scheme is close to the invariant measure of the SDE. In , the first term of an expansion of the invariant measure of the Euler scheme with respect to $\tau$ is also given. In our work, we provide the expansion at any order.

We emphasize that in our result, there is no particular smallness assumption on the stepsize $\tau$ used to define the numerical solution. In particular, the discrete process is not supposed to have a unique invariant measure, as in or .

This is also the case in the recent work . There, it is shown that given an elliptic or hypoelliptic SDE, the ergodic averages provided by the Euler scheme are asymptotically close to the average of the invariant measure of the SDE. Higher order schemes also considered. The main tool in is the ellipticity or hypoellipticity of the Poisson equation, ie the equation $L(x,\partial_{x})u=g$ .

Preliminaries

We consider the stochastic differential equation

We denote by $L(x,\partial_{x})$ the Kolmogorov generator associated with the stochastic equation:

where we use the summation convention for repeated indices and $\partial_{i}=\partial_{x_{i}}$ and

It is well known that the Kolmogorov equation:

We wish to investigate the approximation properties of the Euler scheme for long times. We need assumptions on the long time behavior of the law of the solutions of (2.1), ie of the law of the Markov process. We assume the following mixing properties:

These hypothesis are usually satisfied under elliptic or hypoelliptic assumptions on the operator $L(x,\partial_{x})$ . The reader may find in conditions to ensure [H1]. We also refer to for a general definition of hypoelliticity and applications to numerical schemes. Combining kernel estimates for hypoelliptic diffusion () and exponential convergence to equilibrium, [H3] can be proved. Note that similar estimates are used in , where specific examples are considered. Finally, we mention that these hypothesis can be proved to be fulfilled using partial differential equations techniques (see ).

for $n\geq 0$ . Our main result can be stated as follows:

Let $N$ and $\tau_{0}>0$ be fixed. Then there exists a modified smooth density

Using this result, we can also recover the weak convergence result

In the next sections, the constants appearing in the estimate depend in general on bounds on derivatives of $f$ and $g$ defining the SDE. They will also depend in general on $\tau_{0}$ and $N$ , but not on $\varphi$ .

Asymptotic expansion of the weak error

We have the formal expansion for small $t$ :

Since the solution $u$ of the Kolmogorov equation is smooth and has its derivatives bounded in terms of the initial data $\varphi$ , the above formal expansion can be justified and we have the following proposition whose proof is easy and left to the reader:

With the Euler scheme defined in (2.6), we associate the continuous process

In this work, we are only interested in the distributions of the solutions and of their approximation. We now examine in detail the first time step and its approximation properties in terms of the law. By Markov property, it is sufficient to then obtain information at all steps. Next result gives an expansion similar to Proposition 3.1 for the Euler process.

Then for all $n\geq 1$ , there exist operators $A_{n}(x,\partial_{x})$ of order $2n$ , such that for all $N\geq 1$ , there exist a constant $C_{N}$ satisfying

Proof. Using (3.3) and the Itô formula, we get for $t\leq\tau$ ,

Note that the last term is a martingale. We define the operator

Hence applying (3.5) to $\partial_{i}\varphi$ and $\partial_{ij}\varphi$ , we obtain

We set $A_{0}=I$ , $A_{1}=L$ and plugg this in (3.5) to obtain

are operators of order $3$ . Taking the expectation so that the last two term disappear, we easily deduce the result for $N=1$ .

Note the expectation of the last term vanishes so that (3.6) easily implies (3.4).

and obtain (3.6) with $N+1$ replaced by $N+2$ .

Modified generator

Let us now consider $\tau$ as fixed. We want to construct a formal series

where the operators $A_{n}(x,\partial_{x})$ are defined in Theorem 3.2.

where $\widetilde{A}(\tau)=\sum_{n\geq 1}\tau^{n-1}A_{n}$ .

Note that the (formal) inverse of the series is given by

where the $B_{n}$ are the Bernoulli numbers: see for instance and for a similar analysis involving operators. Hence equations (4.1), (4.2) are equivalent in the sense of formal series to

Identifying the right hand sides of (4.1) and (4.3), we get the following recursion formula

Each of the terms of the above sum is an operator of order $2n+2$ with smooth coefficients and therefore $L_{n}$ is also an operator of order $2n+2$ with smooth coefficients.

Note that (4.2) gives immediately the inverse relation of this formal series equation:

where $\mathds{1}$ denote the constant function equal to $1$ .

2 Approximate solution of the modified flow

For a given $N$ , we have constructed in the previous section an operator

In order to perform weak backward error analysis and estimate recursively the modified invariant law of the numerical process, we should be able to define a solution $v^{N}(t,x)$ of the modified flow

However, in our situation we do not know whether this equation has a solution. This is in contrast with standard backward error analysis where the modified flow can always be defined.

The goal of the next proposition is to give a proper definition of the modified flow (4.7).

with initial conditions $v_{0}(0,x)=\varphi(x)$ and $v_{n}(0,x)=0$ for $n\geq 1$ . For all $N\geq 0$ , setting

(i) There exists a constant $C_{N}$ such that for all time $t\geq 0$ , and all $\tau\geq 0$ ,

(ii) For fixed $\tau_{0}>0$ , there exists a constant $C_{N}$ such that for all $\tau\leq\tau_{0}$ ,

be the right-hand side in (4.8). Then $v_{n}$ is uniquely defined and given by the formula

where the constant $C(t)$ depends on $t$ , $k$ , $n$ and on the coefficient of the equation. Clearly, this type of estimate can be improved in the elliptic case. This proves the first part of the Theorem.

Let us consider the successive time derivatives of the functions $w_{n}(s,x)$ . We have using the definition of $w_{n}$ , for all $s\geq 0$ ,

and we see by induction that for all $m\geq 1$ and $s\geq 0$

Now let us consider the Taylor expansion of $w_{n}(\tau)$ , for $\tau\leq\tau_{0}$ and $n=0,\ldots,N$ ,

Using the bounds on the time derivatives of $w_{n}(s,x)$ , we obtain that for all $\tau\geq 0$ and all $n=0,\ldots,N$ ,

for some constant depending on $N$ , $m$ . After summation in $n$ , and using the expression (4.5) of the operators $A_{n}$ and the definition of $w_{n}$ , we get

The second estimate (ii) is then a consequence of (i) with $t=0$ and (4.13).

Note that in the previous theorem, we have constructed a function $v^{(N)}(t,x)$ which is an approximate solution of (4.7). More precisely, we can easily show that we have for all time $t\geq 0$ ,

Asymptotic expansion of the invariant measure and long time behavior

We now analyze the long time behavior of the solution of the modified equation (4.7). In the following, for a given operator $B(x,\partial_{x})$ , we denote by $B(x,\partial_{x})^{*}$ its adjoint with respect to the $L^{2}$ product. We start by an asymptotic expansion of a formal invariant measure for the numerical scheme.

Let $N\geq 0$ be fixed and $L^{(N)}(\tau;x,\partial_{x})$ the operator defined by (4.6). Then the function

with for all $k$ and all $\tau\in[0,\tau_{0}]$ ,

Proof. Assume that $\mu_{0}=\rho$ and $\mu_{j}(x)$ are known, for $j=0,\ldots,n-1$ with $n\geq 1$ . Consider the equation (5.1) given by

Note that the right-hand side $G_{n}(x)$ is a smooth function satisfying

For all $n$ and $k$ there exists a polynomial $P_{k,n}(t)$ such that for all $t\geq 0$ ,

We claim that $c_{n}(t)$ does not depend on time. Indeed:

by definition of the coefficients $\mu_{n}$ , see (5.1). Note that, thanks to the smoothness properties of all the functions, the computation above is easily justified.

Next, we compute the average of $F_{n}$ . By (4.8), (4.11) and (5.3), we have

Using the previous expression obtained for $\langle F_{n}(s)\rangle$ and recalling the initial data for $v_{n}$ , we deduce that

We give now our main result concerning the long time behavior of the numerical solution:

Let $\tau_{0}$ and $N$ be fixed. Then there exists $C_{N}$ and a polynomial $P_{N}(t)$ such that the following holds: Let $X_{p}$ be the discrete process defined by (2.6), then we have for $p\geq 0$ , $\tau\geq 0$ and smooth function $\varphi(x)$

where for all $p$ , $t_{p}=p\tau$ . Moreover, we have

Proof. For all $p$ , with $t_{j}=j\tau$ , we have

Using (4.10) with $t=t_{j}$ , and Proposition 5.2, we deduce that