Error estimates of the backward Euler-Maruyama method for multi-valued stochastic differential equations

TL;DR

This paper establishes error estimates for the backward Euler-Maruyama method applied to multi-valued stochastic differential equations, including those with non-smooth convex potentials, demonstrating convergence of order at least 1/4.

Contribution

It provides the first rigorous error analysis for the backward Euler-Maruyama method in the context of multi-valued SDEs with non-smooth potentials, extending existing techniques.

Findings

Method is well-defined and converges with order at least 1/4

Applicable to stochastic gradient flows with discontinuous gradients

Verified on overdamped Langevin and stochastic p-Laplace equations

Abstract

In this paper, we derive error estimates of the backward Euler-Maruyama method applied to multi-valued stochastic differential equations. An important example of such an equation is a stochastic gradient flow whose associated potential is not continuously differentiable, but assumed to be convex. We show that the backward Euler-Maruyama method is well-defined and convergent of order at least $1/4$ with respect to the root-mean-square norm. Our error analysis relies on techniques for deterministic problems developed in [Nochetto, Savar\'e, and Verdi, Comm.\ Pure Appl.\ Math., 2000]. We verify that our setting applies to an overdamped Langevin equation with a discontinuous gradient and to a spatially semi-discrete approximation of the stochastic $p$-Laplace equation.

Equations

$$\{\begin{array}{l}\mathrm{d}X(t)=-\nabla \Phi (X(t))\mathrm{d}t+g_0\mathrm{d}W(t),t\in (0,T],\\ X(0)=X_0,\end{array}$$

$$\pi =\{0=t_0<\dots <t_n<\dots <t_N=T\}$$

$$\{\begin{array}{l}{X}^n=X^{n-1}-k\nabla \Phi (X^n)+g_0\Delta W^n,n\in \{1,\dots ,N\},\\ X^0=X_0,\end{array}$$

$$\nabla \Phi (x)=\{\begin{array}{ll}-1,& \text{ if }x<0,\\ c,& \text{ if }x=0,\\ 1,& \text{ if }x>0.\end{array}$$

$$\displaystyle\partial\Phi(x)=\big{\{}v\in{\mathbf{R}}^{d}\,:\,\Phi(x)+\langle v,y-x\rangle\leq\Phi(y)\,\text{for all }y\in{\mathbf{R}}^{d}\big{\}}.$$

$$\partial \Phi (x)=\{\begin{array}{ll}\{-1\},& \text{ if }x<0,\\ [-1,1],& \text{ if }x=0,\\ \{1\},& \text{ if }x>0.\end{array}$$

$$\{\begin{array}{l}\mathrm{d}X(t)+f(X(t))\mathrm{d}t\ni b(X(t))\mathrm{d}t+g(X(t))\mathrm{d}W(t),t\in (0,T],\\ X(0)=X_0.\end{array}$$

$$\{\begin{array}{l}{X}^n\in X^{n-1}-kf(X^n)+kb(X^n)+g(X^{n-1})\Delta W^n,n\in \{1,\dots ,N\},\\ X^0=X_0.\end{array}$$

$$⟨f_v-f_z,z-w⟩\le \gamma ⟨f_v-f_w,v-w⟩$$

$$D(f)=\{x\in M:f(x)\ne \varnothing \}.$$

$$⟨f_u-f_v,u-v⟩\ge 0$$

$$⟨y-f_v,x-v⟩\ge 0\text{ for all }v\in D(f),f_v\in f(v),$$

$$\displaystyle{\mathbf{E}}\Big{[}\Big{|}\int_{s}^{t}g(\tau)\,\mathrm{d}W(\tau)\Big{|}^{p}\Big{]}$$

$$Z(t)\le \xi (t)+{\int }_0^t\phi (s)Z(s)\mathrm{d}s+M(t),\mathbf{P}\text{-almost surely,}$$

$$\displaystyle{\mathbf{E}}\big{[}Z(t)\big{]}\leq\exp\Big{(}\int_{0}^{t}\varphi(s)\,\mathrm{d}s\Big{)}{\mathbf{E}}\big{[}\sup_{s\in[0,t]}\xi(s)\big{]}$$

$$\{\begin{array}{l}\mathrm{d}X(t)=-\nabla \Phi (X(t))\mathrm{d}t+g_0\mathrm{d}W(t),t\in [0,T],\\ X(0)=X_0,\end{array}$$

$$⟨f(v),w-v⟩\le \Phi (w)-\Phi (v),v,w\in {\mathbf{R}}^d,$$

$$⟨f(v)-f(w),v-w⟩\ge 0$$

$$⟨f(v)-f(z),z-w⟩\le ⟨f(v)-f(w),v-w⟩$$

$$⟨f(v)-f(w),v-w⟩$$

$$\displaystyle\geq-\big{(}\Phi(w)-\Phi(v)\big{)}-\big{(}\Phi(v)-\Phi(w)\big{)}=0$$

$$⟨f(v)-f(z),z-w⟩$$

$$=⟨f(v),z⟩-⟨f(v),w⟩+⟨f(z),w-z⟩\pm ⟨f(v),v⟩$$

$$=⟨f(v),z-v⟩+⟨f(v),v-w⟩+⟨f(z),w-z⟩\pm ⟨f(w),v-w⟩$$

$$=⟨f(v),z-v⟩+⟨f(v)-f(w),v-w⟩+⟨f(z),w-z⟩+⟨f(w),v-w⟩.$$

$$⟨f(v)-f(z),z-w⟩$$

$$+\Phi (w)-\Phi (z)-\sigma (z,w)+\Phi (v)-\Phi (w)-\sigma (w,v)$$

$$=⟨f(v)-f(w),v-w⟩-\sigma (v,z)-\sigma (z,w)-\sigma (w,v).$$

$$X(t)=X_0-{\int }_0^{t}f(X(s))\mathrm{d}s+g_{0}W(t).$$

$$\displaystyle\sup_{t\in[0,T]}{\mathbf{E}}\big{[}|X(t)|^{2}\big{]}\leq C\big{(}1+{\mathbf{E}}\big{[}|X_{0}|^{2}\big{]}\big{)}.$$

$$\{\begin{array}{l}{X}^n=X^{n-1}-kf(X^n)+g_0\Delta W^n,n\in \{1,\dots ,N\},\\ X^0=X_0,\end{array}$$

$$\displaystyle\max_{n\in\{1,\ldots,N\}}{\mathbf{E}}[|X^{n}|^{2}]\leq{\mathbf{E}}\big{[}|X_{0}|^{2}\big{]}+2T\big{(}\Phi(0)+|g_{0}|^{2}\big{)}$$

$$\displaystyle\sum_{n=1}^{N}{\mathbf{E}}\big{[}|X^{n}-X^{n-1}|^{2}\big{]}+4k\sum_{n=1}^{N}{\mathbf{E}}\big{[}\Phi(X^{n})\big{]}\leq 2{\mathbf{E}}\big{[}|X_{0}|^{2}\big{]}+4T\big{(}\Phi(0)+|g_{0}|^{2}\big{)}.$$

Full text

1. Introduction

In this paper, we investigate the numerical approximation of multi-valued stochastic differential equations (MSDE). An important example of such equations is provided by stochastic gradient flows with a convex potential. More precisely, let $T\in(0,\infty)$ and $(\Omega,{\mathcal{F}},({\mathcal{F}}_{t})_{t\in[0,T]},{\mathbf{P}})$ be a filtered probability space satisfying the usual conditions. By $W\colon[0,T]\times\Omega\to{\mathbf{R}}^{m}$, $m\in{\mathbf{N}}$, we denote a standard $({\mathcal{F}}_{t})_{t\in[0,T]}$-adapted Wiener process. For instance, let us consider the numerical treatment of nonlinear, overdamped Langevin-type equations of the form

[TABLE]

where $X_{0}\in L^{p}(\Omega,{\mathcal{F}}_{0},{\mathbf{P}};{\mathbf{R}}^{d})$, $p\in[2,\infty)$, $g_{0}\in{\mathbf{R}}^{d,m}$, and $\Phi\colon{\mathbf{R}}^{d}\to{\mathbf{R}}$ are given. These equations have many important applications, for example, in Bayesian statistics and molecular dynamics. We refer to [10, 22, 23, 44, 49] and the references therein.

We recall that, if the gradient $\nabla\Phi$ is of superlinear growth, then the classical forward Euler–Maruyama method is known to be divergent in the strong and weak sense, see [18]. This problem can be circumvented by using modified versions of the explicit Euler–Maruyama method based on techniques such as taming, truncating, stopping, projecting, or adaptive strategies, cf. [4, 6, 17, 19, 29, 48].

In this paper, we take an alternative approach by considering the backward Euler–Maruyama method. Our main motivation for considering this method lies in its good stability properties, which allow its application to stiff problems arising, for instance, from the spatial semi-discretization of stochastic partial differential equations. Implicit methods have also been studied extensively in the context of stochastic differential equations with superlinearly growing coefficients. For example, see [1, 15, 16, 30, 31].

The error analysis in the above mentioned papers on explicit and implicit methods typically requires a certain degree of smoothness of $\nabla\Phi$ such as local Lipschitz continuity. The purpose of this paper is to derive error estimates of the backward Euler–Maruyama method for equations of the form (1.1), where the associated potential $\Phi\colon{\mathbf{R}}^{d}\to{\mathbf{R}}$ is not necessarily continuously differentiable, but assumed to be convex.

For the formulation of the numerical scheme, let $N\in{\mathbf{N}}$ be the number of temporal steps, let $k=\frac{T}{N}$ be the step size, and let

[TABLE]

be an equidistant partition of the interval $[0,T]$, where $t_{n}=nk$ for $n\in\{0,\dots,N\}$. The backward Euler–Maruyama method for the Langevin equation (1.1) is then given by the recursion

[TABLE]

where $\Delta W^{n}=W(t_{n})-W(t_{n-1})$.

An example of a non-smooth potential is found by setting $d=m=1$ and $\Phi(x)=|x|^{p}$, $x\in{\mathbf{R}}$, for $p\in[1,2)$. Evidently, the gradient of $\Phi$ is not locally Lipschitz continuous at $0\in{\mathbf{R}}$ for $p\in(1,2)$. Moreover, if $p=1$, then the gradient $\nabla\Phi$ has a jump discontinuity of the form

[TABLE]

Here, the value $c\in{\mathbf{R}}$ at $x=0$ is not canonically determined. We have to solve a nonlinear equation of the form $x+k\nabla\Phi(x)=y$ in each step of the backward Euler method (1.3). However, if $y\in(-k,k)$, then the sole candidate for a solution is $x=0$, since otherwise $|x+k\nabla\Phi(x)|\geq k$. But $x=0$ is only a solution if $kc=y$. Therefore, the mapping ${\mathbf{R}}\ni x\mapsto x+k\nabla\Phi(x)\in{\mathbf{R}}$ is not surjective for any single-valued choice of $c$.

This problem can be bypassed by considering the multi-valued subdifferential $\partial\Phi\colon{\mathbf{R}}^{d}\to 2^{{\mathbf{R}}^{d}}$ of a convex potential $\Phi\colon{\mathbf{R}}^{d}\to{\mathbf{R}}$, which is given by

[TABLE]

Recall that $\partial\Phi(x)=\{\nabla\Phi(x)\}$ if the gradient exists at $x\in{\mathbf{R}}^{d}$ in the classical sense. See [45, Section 23] for further details.

In the above example, one easily verifies that

[TABLE]

This allows us to solve the nonlinear inclusion where we want to find $x\in{\mathbf{R}}$ with $x+k\partial\Phi(x)\ni y$ for any $y\in{\mathbf{R}}$.

For this reason, we study the more general problem of the numerical approximation of multi-valued stochastic differential equations (MSDE) of the form

[TABLE]

Here, we assume that the mappings $b\colon{\mathbf{R}}^{d}\to{\mathbf{R}}^{d}$ and $g\colon{\mathbf{R}}^{d}\to{\mathbf{R}}^{d,m}$ are globally Lipschitz continuous. Moreover, the multi-valued drift coefficient function $f\colon{\mathbf{R}}^{d}\to 2^{{\mathbf{R}}^{d}}$ is assumed to be a maximal monotone operator, cf. Definition 2.1 below. See also Section 4 for a complete list of all imposed assumptions on the MSDE (1.5). Let us emphasize that the subdifferential of a proper, lower semi-continuous and convex potential is an important example of a possibly multi-valued and maximal monotone mapping $f$, cf. [45, Corollary 31.5.2].

The backward Euler–Maruyama method for the approximation of the MSDE (1.5) on the partition $\pi$ is then given by the recursion

[TABLE]

We discuss the well-posedness of this method (1.6) under our assumptions on $f$, $b$, and $g$ in Section 5. In particular, it will turn out that both problems, (1.5) and (1.6), admit single-valued solutions $(X(t))_{t\in[0,T]}$ and $(X^{n})_{n=0}^{N}$, respectively.

The main result of this paper, Theorem 6.4, then states that the backward Euler–Maruyama method is convergent of order at least $1/4$ with respect to the norm in $L^{2}(\Omega;{\mathbf{R}}^{d})$. For the error analysis we rely on techniques for deterministic problems developed in [37]. An important ingredient is the additional condition on $f$ that there exists $\gamma\in(0,\infty)$ with

[TABLE]

for all $v,w,z\in D(f)\subset{\mathbf{R}}^{d}$ and $f_{v}\in f(v)$, $f_{w}\in f(w)$, $f_{z}\in f(z)$. This assumption is easily verified for a subdifferential of a convex potential, cf. Lemma 3.2. As already noted in [37] for deterministic problems, this inequality allows us to avoid Gronwall-type arguments in the error analysis for terms involving the multi-valued mapping $f$.

Before we give a more detailed outline of the content of this paper let us mention that multi-valued stochastic differential equations have been studied in the literature before. The existence of a uniquely determined solution to the MSDE (1.5) has been investigated, e.g., in [7, 21, 41]. We also refer to the more recent monograph [40] and the references therein. In [14, 51] related results have been derived for multi-valued stochastic evolution equations in infinite dimension. The numerical analysis for MSDEs has also been considered in [3, 26, 42, 53, 55]. However, these papers differ from the present paper in terms of the considered numerical methods, the imposed conditions, or the obtained order of convergence.

Further, we also mention that several authors have developed explicit numerical methods for SDEs with discontinuous drifts in recent years. For instance, we refer to [9, 24, 25, 33, 34, 35, 36]. While these results often apply to more irregular drift coefficients, which are beyond the framework of maximal monotone operators, the authors have to employ more restrictive conditions such as the global boundedness of the drift, which is not required in our framework.

This paper is organized as follows: In Section 2 we fix some notation and recall important terminology for multi-valued mappings. In Section 3 we demonstrate how to apply the techniques from [37] to the simplified setting of the Langevin equation (1.1). In addition, we also show that if the gradient $\nabla\Phi$ is more regular, say Hölder continuous with exponent $\alpha\in(0,1]$, then the order of convergence increases to $\frac{1+\alpha}{4}$. Moreover, it turns out that the error constant does not grow exponentially with the final time $T$. This is an important insight if the backward Euler method is used within an unadjusted Langevin algorithm [44], which typically requires large time intervals. See Theorem 3.7 and Remark 3.8 below.

In Section 4 we turn to the more general multi-valued stochastic differential equation (1.5). We state all assumptions imposed on the appearing drift and diffusion coefficients and collect some properties of the exact solution. In Section 5 we show that the backward Euler–Maruyama method (1.6) is well-posed under the assumptions of Section 4. In Section 6 we prove the already mentioned convergence result with respect to the root-mean-square norm. Finally, in Section 7 we verify that the setting of Section 4 applies to a Langevin equation with the discontinuous gradient (1.4). Further, we also show how to apply our results to the spatial discretization of the stochastic $p$-Laplace equation which indicates their usability for the numerical analysis of stochastic partial differential equations. However, a complete analysis of the latter problem will be deferred to a future work.

2. Preliminaries

In this section, we collect some notation and introduce some background material. First we recall some terminology for set valued mappings and (maximal) monotone operators. For a more detailed introduction we refer, for instance, to [47, Abschn. 3.3] or [39, Chapter 6].

By ${\mathbf{R}}^{d}$, $d\in{\mathbf{N}}$, we denote the Euclidean space with the standard norm $|\cdot|$ and inner product $\langle\cdot,\cdot\rangle$. Let $M\subset{\mathbf{R}}^{d}$ be a set. A set-valued mapping $f\colon M\to 2^{{\mathbf{R}}^{d}}$ maps each $x\in M$ to an element of the power set $2^{{\mathbf{R}}^{d}}$, that is, $f(x)\subseteq{\mathbf{R}}^{d}$. The domain $D(f)$ of $f$ is given by

[TABLE]

Definition 2.1**.**

Let $M\subset{\mathbf{R}}^{d}$ be a non-empty set. A set-valued mapping $f\colon M\to 2^{{\mathbf{R}}^{d}}$ is called monotone if

[TABLE]

for all $u,v\in D(f)$, $f_{u}\in f(u)$, and $f_{v}\in f(v)$.

Moreover, a set-valued mapping $f\colon M\to 2^{{\mathbf{R}}^{d}}$ is called maximal monotone if $f$ is monotone and for all $x\in M$ and $y\in{\mathbf{R}}^{d}$ satisfying

[TABLE]

it follows that $x\in D(f)$ and $y\in f(x)$.

Next, we recall a Burkholder–Davis–Gundy-type inequality. For a proof we refer to [28, Chapter 1, Theorem 7.1]. For its formulation we take note that the Frobenius or Hilbert–Schmidt norm of a matrix $g\in{\mathbf{R}}^{d,m}$ is also denoted by $|g|$.

Lemma 2.2**.**

Let $p\in[2,\infty)$ and $g\in L^{p}(\Omega;L^{p}(0,T;{\mathbf{R}}^{d,m}))$ be stochastically integrable. Then, for every $s,t\in[0,T]$ with $s<t$, the inequality

[TABLE]

holds.

Let us also recall a stochastic variant of the Gronwall inequality. A proof that can be modified to this setting can be found in [54]. Compare also with [50].

Lemma 2.3**.**

Let $Z,M,\xi\colon[0,T]\times\Omega\to{\mathbf{R}}$ be $({\mathcal{F}}_{t})_{t\in[0,T]}$-adapted and almost surely continuous stochastic processes such that $M$ is a local $({\mathcal{F}}_{t})_{t\in[0,T]}$-martingale with $M(0)=0$. Moreover, suppose that $Z$ and $\xi$ are nonnegative. In addition, let $\varphi\colon[0,T]\to{\mathbf{R}}$ be integrable and nonnegative. If, for all $t\in[0,T]$, we have

[TABLE]

then, for every $t\in[0,T]$, the inequality

[TABLE]

holds.

Moreover, we often make use of generic constants. More precisely, by $C$ we denote a finite and positive quantity that may vary from occurrence to occurrence but is always independent of numerical parameters such as the step size $k=\frac{T}{N}$ and the number of steps $N\in{\mathbf{N}}$.

3. Application to the Langevin equation with a convex potential

In order to illustrate our approach, we first consider a more regular stochastic differential equation with single-valued (Hölder) continuous drift term. More precisely, we consider the overdamped Langevin equation [23, Section 2.2]

[TABLE]

where $X_{0}\in L^{2}(\Omega,{\mathcal{F}}_{0},{\mathbf{P}};{\mathbf{R}}^{d})$, $g_{0}\in{\mathbf{R}}^{d,m}$, and $W\colon[0,T]\times\Omega\to{\mathbf{R}}^{m}$ is a standard ${\mathbf{R}}^{m}$-valued Wiener process. In addition, we impose the following assumption on the potential $\Phi\colon{\mathbf{R}}^{d}\to{\mathbf{R}}$.

Assumption 3.1**.**

Let $\Phi\colon{\mathbf{R}}^{d}\to{\mathbf{R}}$ be a convex, nonnegative, and continuously differentiable function.

In the following, we denote by $f\colon{\mathbf{R}}^{d}\to{\mathbf{R}}^{d}$ the gradient of $\Phi$, that is $f(x)=\nabla\Phi(x)$. It is well-known that the convexity of $\Phi$ implies the variational inequality

[TABLE]

see, for example, [45, § 23].

In the following lemma, we collect some properties of $f$, which are direct consequences of Assumption 3.1. Both inequalities are well-known. The proof of (3.4) is taken from [37].

Lemma 3.2**.**

Under Assumption 3.1 and with $f=\nabla\Phi$, the inequalities

[TABLE]

and

[TABLE]

are fulfilled for all $v,w,z\in{\mathbf{R}}^{d}$.

Proof.

The first inequality follows directly from (3.2) since

[TABLE]

for all $v,w\in{\mathbf{R}}^{d}$. For the proof of the second inequality we start by rewriting its left-hand side. For arbitrary $v,w,z\in{\mathbf{R}}^{d}$ we rearrange the terms to obtain

[TABLE]

Setting $\sigma(v,w):=\Phi(w)-\Phi(v)-\langle f(v),w-v\rangle$ for all $v,w\in{\mathbf{R}}^{d}$, we see that

[TABLE]

But (3.2) says that $\sigma(v,w)\geq 0$ for all $v,w\in{\mathbf{R}}^{d}$, which completes the proof. ∎

It follows from Assumption 3.1 and Lemma 3.2 that the drift $f=\nabla\Phi$ of the stochastic differential equation (3.1) is continuous and monotone. Therefore, by [43, Thm. 3.1.1] the stochastic differential equation (3.1) has a solution in the strong (probabilistic) sense satisfying ${\mathbf{P}}$-a.s. for all $t\in[0,\infty)$

[TABLE]

Moreover, the solution is unique up to ${\mathbf{P}}$-indistinguishability and it is square-integrable with

[TABLE]

Next, we turn to the numerical approximation of the solution of (3.1). Recall that for a single-valued drift the backward Euler–Maruyama method is given by the recursion

[TABLE]

where $\Delta W^{n}=W(t_{n})-W(t_{n-1})$, $t_{n}=nk$ and $k=\frac{T}{N}$.

The next lemma contains some a priori estimates for the backward Euler–Maruyama method (3.6).

Lemma 3.3**.**

Let $g_{0}\in{\mathbf{R}}^{d,m}$ be given and let Assumption 3.1 be satisfied. For an arbitrary step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, let $(X^{n})_{n\in\{0,\dots,N\}}$ be a family of $({\mathcal{F}}_{t_{n}})_{n\in\{0,\ldots,N\}}$-adapted random variables satisfying (3.6). If $X_{0}\in L^{2}(\Omega,{\mathcal{F}}_{0},{\mathbf{P}};{\mathbf{R}}^{d})$, then

[TABLE]

and

[TABLE]

Proof.

First, we recall the identity

[TABLE]

Using also (3.6), we then get

[TABLE]

for every $n\in\{1,\ldots,N\}$. Hence, an application of (3.2) yields

[TABLE]

for every $n\in\{1,\ldots,N\}$. From applications of the Cauchy–Schwarz inequality and the weighted Young inequality we then obtain

[TABLE]

for every $n\in\{1,\ldots,N\}$.

The third term on the right-hand side is absorbed in the third term on the left-hand side. Summation then yields

[TABLE]

An inductive argument over $n\in\{1,\ldots,N\}$ then yields that $X^{n}$ is square-integrable due to the assumption $X_{0}\in L^{2}(\Omega,{\mathcal{F}}_{0},{\mathbf{P}};{\mathbf{R}}^{d})$. Therefore, after taking expectation the last sum vanishes. Moreover, an application of the Itō isometry then gives

[TABLE]

Since this is true for any $n\in\{1,\ldots,N\}$ the assertion follows. ∎

As the next theorem shows, Assumption 3.1 is also sufficient to ensure the well-posedness of the backward Euler–Maruyama method. The result follows directly from the fact that $f$ is continuous and monotone due to (3.3). For a proof we refer, for instance, to [4, Sect. 4], [38, Chap. 6.4], and [52, Theorem C.2]. The assertion also follows from the more general result in Theorem 5.3 below.

Theorem 3.4**.**

Let $X_{0}\in L^{2}(\Omega,{\mathcal{F}}_{0},{\mathbf{P}};{\mathbf{R}}^{d})$ as well as $g_{0}\in{\mathbf{R}}^{d,m}$ be given and let Assumption 3.1 be satisfied. Then, for every equidistant step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, there exists a uniquely determined family of square-integrable and $({\mathcal{F}}_{t_{n}})_{n\in\{0,\ldots,N\}}$-adapted random variables $(X^{n})_{n\in\{0,\dots,N\}}$ satisfying (3.6).

We now turn to an error estimate with respect to the $L^{2}(\Omega;{\mathbf{R}}^{d})$-norm. Since we do not impose any (local) Lipschitz condition on the drift $f$, classical approaches based on discrete Gronwall-type inequalities are not applicable. Instead we rely on an error representation formula, which was introduced for deterministic problems in [37].

For the formulation of this, we introduce the following notation: For a given equidistant partition $\pi=\{0=t_{0}<t_{1}<\ldots<t_{N}=T\}\subset[0,T]$ with step size $k=\frac{T}{N}$, we denote by ${\mathcal{X}}\colon[0,T]\times\Omega\to{\mathbf{R}}^{d}$ the piecewise linear interpolant of the sequence $(X^{n})_{n\in\{0,\ldots,N\}}$ generated by the backward Euler method (3.6). It is defined by ${\mathcal{X}}(0)=X^{0}$ and

[TABLE]

In addition, we introduce the processes $\overline{{\mathcal{X}}},\underline{{\mathcal{X}}}\colon[0,T]\times\Omega\to{\mathbf{R}}^{d}$, which are piecewise constant interpolants of $(X^{n})_{n\in\{0,\ldots,N\}}$ and defined by $\overline{{\mathcal{X}}}(0)=\underline{{\mathcal{X}}}(0)=X^{0}$ and

[TABLE]

Analogously, we define the piecewise linear interpolated process $\mathcal{W}\colon[0,T]\times\Omega\to{\mathbf{R}}^{m}$ by $\mathcal{W}(0)=0$ and

[TABLE]

for all $t\in(t_{n-1},t_{n}]$, $n\in\{1,\ldots,N\}$.

We are now prepared to state Lemma 3.5. The underlying idea of this lemma was introduced in [37], where it is used to derive a posteriori error estimates for the backward Euler method. In fact, in the absence of noise, only the first term on the right-hand side of (3.12) is non-zero. In [37] this term is used as an a posteriori error estimator, since it is explicitly computable by quantities generated by the numerical method.

Lemma 3.5**.**

Let $X_{0}\in L^{2}(\Omega,{\mathcal{F}}_{0},{\mathbf{P}};{\mathbf{R}}^{d})$ as well as $g_{0}\in{\mathbf{R}}^{d,m}$ be given and let Assumption 3.1 be satisfied. Let $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, be an arbitrary equidistant step size and let $t_{n}=nk$, $n\in\{0,\dots,N\}$. Then, for every $n\in\{1,\ldots,N\}$ the estimate

[TABLE]

holds, where $(X(t))_{t\in[0,T]}$ and $(X^{n})_{n\in\{0,\ldots,N\}}$ are the solutions of (3.1) and (3.6), respectively.

Proof.

From (3.6) we directly deduce that for every $n\in\{1,\ldots,N\}$

[TABLE]

Then, one easily verifies for all $t\in(t_{n-1},t_{n}]$, $n\in\{1,\ldots,N\}$, that

[TABLE]

Hence, due to (3.5) the error process $E:=X-{\mathcal{X}}$ fulfills

[TABLE]

for all $t\in[0,T]$. Here, we have $E_{2}(t_{n})=0$, since $\mathcal{W}$ is an interpolant of $W$. Hence, for all $n\in\{0,\ldots,N\}$,

[TABLE]

To estimate the norm of $E_{1}(t_{n})$, we first note that $E_{1}$ has absolutely continuous sample paths with $E_{1}(0)=0$. Hence,

[TABLE]

is fulfilled for almost all $t\in[0,T]$. Therefore, by integration with respect to $t$, we get

[TABLE]

Next, we write

[TABLE]

and use (3.3) and (3.4) to obtain, for almost every $t\in(t_{n-1},t_{n}]$, that

[TABLE]

Furthermore, the expectation of the second integral on the right-hand side of (3.15) is equal to

[TABLE]

Therefore,

[TABLE]

Since $\int_{t_{i-1}}^{t_{i}}(t_{i}-t)\,\mathrm{d}t=\frac{1}{2}k^{2}$ the assertion follows. ∎

The next lemma contains an estimate of the difference between the Wiener process $W$ and its piecewise linear interpolant $\mathcal{W}$.

Lemma 3.6**.**

For every $g_{0}\in{\mathbf{R}}^{d,m}$ and every step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, the equality

[TABLE]

holds.

Proof.

From the definition (3.11) of $\mathcal{W}$ it follows that

[TABLE]

where we used that the two increments of the Wiener process are independent for every $t\in(t_{n-1},t_{n}]$, $n\in\{1,\ldots,N\}$, and we also applied Itō’s isometry. By symmetry of the two terms it then follows that

[TABLE]

and the proof is complete. ∎

The error estimates in Lemma 3.5 and Lemma 3.6 allow us to determine the order of convergence of the backward Euler–Maruyama method without relying on discrete Gronwall-type inequalities. The following theorem imposes the additional assumption that the drift $f$ is Hölder continuous. We include the parameter value $\alpha=0$, which simply means that $f$ is continuous and globally bounded. The case of less regular $f$ is treated in Section 6.

Observe that we recover the standard rate $\frac{1}{2}$ if $\alpha=1$, that is, if the drift $f$ is assumed to be globally Lipschitz continuous. Compare also with the standard literature, for example, [20, Chap. 12] or [32, Sect. 1.3].

For processes $X\colon[0,T]\times\Omega\to{\mathbf{R}}^{d}$ and exponents $\alpha\in[0,1]$, we define the family of Hölder semi-norms by

[TABLE]

Theorem 3.7**.**

Let $X_{0}\in L^{2}(\Omega,{\mathcal{F}}_{0},{\mathbf{P}};{\mathbf{R}}^{d})$ as well as $g_{0}\in{\mathbf{R}}^{d,m}$ be given, let Assumption 3.1 be fulfilled and let $f=\nabla\Phi$ be Hölder continuous with exponent $\alpha\in[0,1]$, i.e., there exists $L_{f}\in(0,\infty)$ such that

[TABLE]

Then there exists $C\in(0,\infty)$ such that for every step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, the estimate

[TABLE]

holds, where $(X(t))_{t\in[0,T]}$ and $(X^{n})_{n\in\{0,\ldots,N\}}$ are the solutions to (3.1) and (3.6), respectively.

Proof.

Since $f$ is assumed to be $\alpha$-Hölder continuous it follows that

[TABLE]

In particular, $f$ grows at most linearly. Therefore, as stated in [28, Chap. 2, Thm 4.3], the solution $(X(t))_{t\in[0,T]}$ of (3.1) satisfies $X\in C^{\frac{1}{2}}([0,T];L^{2}(\Omega;{\mathbf{R}}^{d}))$.

We will use Lemma 3.5 to prove the error bound. To this end, we first show that

[TABLE]

Indeed, we make use of the Hölder continuity of $f$ directly and obtain

[TABLE]

where we also used Hölder’s inequality with $p=\frac{2}{1+\alpha}\in[1,2]$ and $\frac{1}{p}+\frac{1}{q}=1$ as well as Jensen’s inequality. Due to the a priori estimate (3.8) the sum \sum_{i=1}^{N}{\mathbf{E}}\big{[}|X^{i}-X^{i-1}|^{2}\big{]} is bounded independently of the step size $k$. Hence, we arrive at (3.17).

Therefore, it remains to estimate the second error term in Lemma 3.5:

[TABLE]

where we inserted the definition of $\overline{{\mathcal{X}}}$ from (3.10). Moreover, from (3.11) we get

[TABLE]

for $t\in(t_{j-1},t_{j}]$. Hence, the random variable in the second slot of the inner product on the right-hand side (3.18) is centered and is independent of any ${\mathcal{F}}_{t_{j-1}}$-measurable random variable. Thus, we may write

[TABLE]

To estimate $T_{1}$ we first recall the definitions of $\underline{{\mathcal{X}}}$ and $\overline{{\mathcal{X}}}$ from (3.10). Then we apply the Cauchy–Schwarz inequality and obtain

[TABLE]

From the Hölder continuity of $f$ we then deduce that

[TABLE]

where the last inequality is in fact an equality if $\alpha=1$, $\frac{1}{q}=0$ or if $\alpha=0$, $\frac{1}{q}=1$. Otherwise the inequality follows from Hölder’s inequality with $p=\frac{1}{\alpha}\in(1,\infty)$ and $\frac{1}{p}+\frac{1}{q}=1$, followed by an application of Jensen’s inequality. Furthermore, Lemma 3.6 states that

[TABLE]

Therefore, together with (3.8) we arrive at the estimate

[TABLE]

for all $n\in\{1,\ldots,N\}$.

The estimate of $T_{2}$ works similarly by additionally making use of the Hölder continuity of the exact solution. To be more precise, we have that

[TABLE]

Together with the Cauchy–Schwarz inequality and (3.19), we therefore obtain

[TABLE]

Inserting the estimates for $T_{1}$, $T_{2}$ and (3.17) into Lemma 3.5 completes the proof. ∎

Remark 3.8**.**

The precise form of the constant $C$ appearing in Theorem 3.7 is, after taking squares,

[TABLE]

with $C_{0}=2{\mathbf{E}}[|X_{0}|^{2}]+4T(\Phi(0)+|g_{0}|^{2})$.

Observe that, since we avoid the use of Gronwall-type inequalities, the error constant does not grow exponentially with time $T$. This indicates that the backward Euler–Maruyama method is particularly suited for long-time simulations as is often required in Markov-chain Monte Carlo methods, for example, in the unadjusted Langevin algorithm [44].

4. Properties of the exact solution

In this section, we turn our attention to the multi-valued stochastic differential equation (MSDE) in (1.5). We give a complete account of the assumptions imposed on the coefficient functions. In addition, we collect some results on the existence and uniqueness of a strong solution to the MSDE. We also include useful results on higher moment bounds of the exact solution.

Assumption 4.1**.**

The set valued mapping $f\colon{\mathbf{R}}^{d}\to 2^{{\mathbf{R}}^{d}}$ is maximal monotone with $\operatorname{int}D(f)\neq\emptyset$. Moreover, there exist constants $\beta,\lambda\in[0,\infty)$, $\mu\in(0,\infty)$, and $p\in[1,\infty)$ such that

[TABLE]

for every $v\in D(f)$ and $f_{v}\in f(v)$.

Assumption 4.2**.**

The function $b\colon{\mathbf{R}}^{d}\to{\mathbf{R}}^{d}$ is Lipschitz continuous; i.e., there exists a constant $L_{b}\in[0,\infty)$ such that

[TABLE]

for all $v,w\in{\mathbf{R}}^{d}$.

Assumption 4.3**.**

The function $g\colon{\mathbf{R}}^{d}\to{\mathbf{R}}^{d,m}$ is Lipschitz continuous; i.e., there exists a constant $L_{g}\in[0,\infty)$ such that

[TABLE]

for all $v,w\in{\mathbf{R}}^{d}$.

Assumption 4.4**.**

The initial value $X_{0}$ is an ${\mathcal{F}}_{0}$-measurable and $D(f)$-valued random variable. Furthermore,

[TABLE]

where the value of $p$ is the same as in Assumption 4.1.

Observe that Assumptions 4.2 and 4.3 directly imply that $b$ and $g$ grow at most linearly. More precisely, after possibly increasing the values of $L_{b}$ and $L_{g}$, we obtain the bounds

[TABLE]

for all $v\in{\mathbf{R}}^{d}$.

Remark 4.5**.**

Without loss of generality we will assume that $0\in D(f)$. Otherwise, since the graph of $f$ is not empty, we take $v_{0}\in D(f)$ and $f_{v_{0}}\in f(v_{0})$ and replace $f$, $b$, and $g$ by suitably shifted mappings, for instance, $\tilde{f}(v):=f(v+v_{0})$. Then $0\in D(\tilde{f})$ holds. Compare further with [47, Abschn. 3.3.3].

Next, we introduce the notion of a solution of (1.5), which we use for the remainder of this paper.

Definition 4.6**.**

A tuple $(X,\eta)$ is called a solution of the multi-valued stochastic differential equation (1.5), if the following conditions hold.

• (i)

The mapping $X\colon[0,T]\times\Omega\to{\mathbf{R}}^{d}$ is an $({\mathcal{F}}_{t})_{t\in[0,T]}$-adapted, almost surely continuous stochastic process such that $X(t)\in\overline{D(f)}$ for all $t\in(0,T]$ with probability one.

• (ii)

The mapping $\eta\colon[0,T]\times\Omega\to{\mathbf{R}}^{d}$ is an $({\mathcal{F}}_{t})_{t\in[0,T]}$-adapted stochastic process such that

[TABLE]

• (iii)

The equality

[TABLE]

holds for all $t\in[0,T]$ and ${\mathbf{P}}$-almost surely.

• (iv)

For almost all $\omega\in\Omega$ and $t\in[0,T]$, it follows that $\eta(t,\omega)\in f(X(t,\omega))$; in other words, for every $y\in D(f)$ and $f_{y}\in f(y)$ the inequality

[TABLE]

is satisfied for almost every $t\in[0,T]$ and ${\mathbf{P}}$-almost surely, cf. Definition 2.1.

This notion of a solution has been considered in, for example, [7], [21], [41], and [51], where also the existence of a unique solution is shown. Due to their importance for the error analysis, we next prove certain moment estimates.

Theorem 4.7**.**

Let Assumptions 4.1 and 4.4 be satisfied with $p\in[1,\infty)$. Then there exists a unique solution $(X,\eta)$ of (1.5) in the sense of Definition 4.6. There is a constant $C\in(0,\infty)$ such that

[TABLE]

Furthermore, if $p\in(1,\infty)$ and $\frac{1}{p}+\frac{1}{q}=1$, then

[TABLE]

Proof.

Existence and uniqueness is shown, for instance, in [21]. For

[TABLE]

the equality

[TABLE]

holds by an application of Itō’s formula (see [12, Chap. 4.7, Theorem 7.1]). From the coercivity assumption on $f$ we obtain that

[TABLE]

for every $f_{X(s)}\in f(X(s))$ and almost every $s\in[0,T]$. The fact that $\eta(s)\in f(X(s))$ for almost every $s\in[0,T]$ then implies that

[TABLE]

Since $b$ and $g$ satisfies the linear growth bound (4.1), we have

[TABLE]

as well as

[TABLE]

Thus, we get

[TABLE]

We introduce

[TABLE]

Then $Z$, $M$, and $\xi$ are $({\mathcal{F}}_{t})_{t\in[0,T]}$-adapted and almost surely continuous stochastic processes. Furthermore, $M$ is a local $({\mathcal{F}}_{t})_{t\in[0,T]}$-martingale with $M(0)=0$. Thus, an application of Lemma 2.3 yields, for every $t\in[0,T]$, that

[TABLE]

Inserting the definition of $Z$ then proves the first estimate.

Furthermore, if Assumption 4.1 holds with $p\in(1,\infty)$, then we have, for every $f_{x}\in f(x)$, $x\in{\mathbf{R}}^{d}$, that

[TABLE]

with $q=\frac{p}{p-1}$. Therefore, it follows that

[TABLE]

since $\eta(s)\in f(X(s))$ for almost every $s\in[0,T]$. ∎

Remark 4.8**.**

Let us mention that, for instance, in [40, Chapter 4] and the references therein, a weaker notion of a solution to (1.5) is found. More precisely, if $(X,\eta)$ is a solution in the sense of Definition 4.6, then $(X,H)$ is a solution in the sense of [40, Chapter 4] with the definition

[TABLE]

In particular, the process $H$ is a continuous, progressively measurable process with bounded total variation and $H(0)=0$ almost surely. The stronger condition of absolute continuity of the process $H$, which is required in Definition 4.6, is essential in the proof of Theorem 6.4 below. This explains why we work with the stronger notion of a solution in Definition 4.6.

5. Well-posedness of the backward Euler method

In this section, we show that the backward Euler–Maruyama method (1.6) for the MSDE (1.5) is well-posed under the same assumptions as in the previous section.

Lemma 5.1**.**

Let Assumptions 4.1 and 4.2 be satisfied. Furthermore, let $w\in{\mathbf{R}}^{d}$ and $k\in(0,T]$ be given with $L_{b}k\in[0,1)$. Then there exist uniquely determined $x_{0}\in D(f)$ and $\eta_{x_{0}}\in f(x_{0})$, which satisfy the nonlinear equation

[TABLE]

Proof.

We first show that there exists a unique $x_{0}\in D(f)$ such that

[TABLE]

To this end, notice that for all $x,y\in{\mathbf{R}}^{d}$, the inequalities

[TABLE]

hold due to the step-size bound. In addition, it follows from (4.1) that

[TABLE]

for all $x\in{\mathbf{R}}^{d}$. Hence, $(\mathrm{id}+kf-kb)$ is the sum of the maximal monotone operator $kf$ and the mapping $(\mathrm{id}-kb)$, which is single-valued, Lipschitz continuous, monotone and coercive.

Thus, we can apply [2, Theorem 2.1] and obtain the existence of $x_{0}\in D(f)$ such that (5.2) holds. Furthermore, there necessarily exists a corresponding unique element $\eta_{x_{0}}\in f(x_{0})$ with

[TABLE]

It remains to prove the uniqueness of $x_{0}$, which directly implies the uniqueness of $\eta_{x_{0}}$. Assume that there exist $x_{1}\in D(f)$ and $\eta_{x_{1}}\in f(x_{1})$ as well as $x_{2}\in D(f)$ and $\eta_{x_{2}}\in f(x_{2})$ such that

[TABLE]

By considering the difference of these equations tested with $x_{1}-x_{2}$, we obtain

[TABLE]

Since $1-kL_{b}>0$ we must have $x_{1}=x_{2}$ and the proof is complete. ∎

For later use, we note that the solution operator for (5.1) is Lipschitz continuous.

Lemma 5.2**.**

Let Assumptions 4.1 and 4.2 be satisfied. For $k\in(0,T]$ with $L_{b}k\in[0,1)$ let $S_{k}\colon{\mathbf{R}}^{d}\to D(f)$ be the solution operator that maps $w\in{\mathbf{R}}^{d}$ to the unique solution $x_{0}\in D(f)$ of (5.1). Then $S_{k}$ is globally Lipschitz continuous with

[TABLE]

Proof.

Let $w_{1},w_{2}\in{\mathbf{R}}^{d}$ and $k\in(0,T]$ with $L_{b}k\in[0,1)$ be given. Let $x_{i}=S_{k}(w_{i})\in D(f)$ and $\eta_{x_{i}}\in f(x_{i})$, $i\in\{1,2\}$, denote the unique solutions of the equations

[TABLE]

By considering the difference of these equations, tested with $x_{1}-x_{2}$, we obtain

[TABLE]

By using the Cauchy–Schwarz inequality for the right-hand side as well as the monotonicity and the Lipschitz continuity for the left-hand side, we get

[TABLE]

Reinserting $x_{i}=S_{k}(w_{i})$ then shows that

[TABLE]

as claimed. ∎

Theorem 5.3**.**

Let Assumptions 4.1 to 4.4 be satisfied. Then for every step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, with $L_{b}k\in[0,1)$ there exist uniquely determined families of square-integrable, ${\mathbf{R}}^{d}$-valued and $({\mathcal{F}}_{t_{n}})_{n\in\{0,\ldots,N\}}$-adapted random variables $(X^{n})_{n\in\{0,\dots,N\}}$ and $(\eta^{n})_{n\in\{0,\dots,N\}}$ such that $X^{n}\in D(f)$, $\eta^{n}\in f(X^{n})$ for every $n\in\{0,\dots,N\}$ and

[TABLE]

for every $n\in\{1,\dots,N\}$, ${\mathbf{P}}$-almost surely.

Proof.

We prove the existence of $(X^{n})_{n\in\{0,\dots,N\}}$ and $(\eta^{n})_{n\in\{0,\dots,N\}}$ by induction over $n\in\{0,\ldots,N\}$. From the assumptions on $X_{0}$ and $f$ it is clear that $X^{0}=X_{0}$ and $\eta^{0}\in f(X_{0})$ are ${\mathcal{F}}_{t_{0}}$-adapted and square-integrable. In particular, it follows from Assumptions 4.1 and 4.4 that

[TABLE]

Next, we assume that $(X^{j})_{j\in\{0,\ldots,n-1\}}$ and $(\eta^{j})_{j\in\{0,\ldots,n-1\}}$ are $({\mathcal{F}}_{t_{j}})_{j\in\{0,\ldots,n-1\}}$-adapted, square-integrable and satisfy (5.3) for all $j\in\{1,\ldots,n-1\}$. By Lemma 5.1 there exist uniquely determined $X^{n}(\omega)\in D(f)$ and $\eta^{n}(\omega)\in f(X^{n}(\omega))$ for almost every $\omega\in\Omega$ such that

[TABLE]

By Lemma 5.2, the solution operator $S_{k}\colon{\mathbf{R}}^{d}\to D(f)$ that maps $X^{n-1}(\omega)+g(X^{n-1}(\omega))\Delta W^{n}(\omega)$ to $X^{n}(\omega)\in D(f)$ is Lipschitz continuous. As $S_{k}$ is Lipschitz continuous and, hence, of linear growth it follows that $X^{n}$ is an ${\mathcal{F}}_{t_{n}}$-measurable and square-integrable random variable. To be more precise, we have the bound

[TABLE]

This implies, in particular, that

[TABLE]

is also a ${\mathcal{F}}_{t_{n}}$-measurable and square-integrable random variable as $X^{n}$, $X^{n-1}$ and $g(X^{n-1})\Delta W^{n}$ have these properties. This finishes the proof of the induction and hence that of the theorem. ∎

Next we state an a priori estimate for the sequence of random variables satisfying recursion (1.6).

Lemma 5.4**.**

Let Assumptions 4.1 to 4.4 be satisfied. For a step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, with $5L_{b}k\in[0,1)$, let $(X^{n})_{n\in\{0,\dots,N\}}$ and $(\eta^{n})_{n\in\{0,\dots,N\}}$ be two families of $({\mathcal{F}}_{t_{n}})_{n\in\{0,\ldots,N\}}$-adapted random variables as stated in Theorem 5.3. Then there exists $K_{X}\in(0,\infty)$ independent of the step size $k=\frac{T}{N}$ such that

[TABLE]

In addition, if $p\in(1,\infty)$, then there exists $K_{\eta}\in(0,\infty)$ independent of the step size $k=\frac{T}{N}$ such that

[TABLE]

where $q\in(1,\infty)$ is given by $\frac{1}{p}+\frac{1}{q}=1$.

Remark 5.5**.**

If $p=1$ in Assumption 4.1, then $f$ and, hence, $(\eta^{n})_{n\in\{1,\ldots,N\}}$ are bounded. In particular, (5.5) holds for any $q\in(1,\infty)$ and for any step size $k=\frac{T}{N}$ with $L_{b}k\in[0,1)$.

Proof of Lemma 5.4.

First, we recall the identity

[TABLE]

As $\eta^{n}\in f(X^{n})$, using Assumptions 4.1 and 4.2, it follows that

[TABLE]

where we also applied (4.1). Hence,

[TABLE]

for every $n\in\{1,\ldots,N\}$, where we also applied the Cauchy–Schwarz and weighted Young inequalities. After a kick-back, we sum from $1$ to $n\in\{1,\ldots,N\}$ to obtain

[TABLE]

After taking expectations, the last term on the right-hand side vanishes. Then, applications of Itō’s isometry and (4.1) give

[TABLE]

Since the step-size bound $5L_{b}k\in[0,1)$ ensures that

[TABLE]

the discrete Gronwall inequality (see, for example, [8]) is applicable and completes the proof of (5.4). Finally, it follows from the polynomial growth bound on $f$ that

[TABLE]

and an application of (5.4) then yields (5.5). ∎

6. Error estimates in the general case

In this section, we derive an error estimate for the backward Euler method given by (1.6) for the MSDE (1.5).

To prove the convergence of the scheme (1.6), let us fix some notation. Throughout this section, we assume that the equidistant step size $k=\frac{T}{N}$ is small enough so that the a priori estimates in Lemma 5.4 hold. Further, as in (3.9) and (3.10), we denote the piecewise linear interpolants of the discrete values by ${\mathcal{X}}(0)=X^{0}$, $\mathcal{H}(0)=\eta^{0}$ for $\eta^{0}\in f(X^{0})$ and

[TABLE]

for all $t\in(t_{n-1},t_{n}]$ and $n\in\{1,\ldots,N\}$. Similarly, we define the piecewise constant interpolant by $\overline{{\mathcal{X}}}(0)=\underline{{\mathcal{X}}}(0)=X^{0}$ and

[TABLE]

for all $t\in(t_{n-1},t_{n}]$ and $n\in\{1,\ldots,N\}$. Moreover, we introduce the stochastic processes $G\colon[0,T]\times\Omega\to{\mathbf{R}}^{d}$ and ${\mathcal{G}}\colon[0,T]\times\Omega\to{\mathbf{R}}^{d}$ defined by

[TABLE]

as well as by ${\mathcal{G}}(0)=0$ and, for all $n\in\{1,\ldots,N\}$ and $t\in(t_{n-1},t_{n}]$,

[TABLE]

In view of (1.6) and the definition of ${\mathcal{G}}$ for $t\in(t_{n-1},t_{n}]$, $n\in\{1,\dots,N\}$, we obtain the representation

[TABLE]

We begin the derivation of our error estimate by considering the difference between the stochastic integral $G$ and its approximation ${\mathcal{G}}$.

Lemma 6.1**.**

Let Assumptions 4.1 to 4.4 be satisfied. Then there exists $K_{G}\in(0,\infty)$ such that, for every equidistant step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$ with $5L_{b}k\in[0,1)$ and every $t\in[0,T]$, we have

[TABLE]

In addition, for every $\rho\in[2,\infty)$, there exists $K_{\rho}\in(0,\infty)$ such that, for every $n\in\{1,\ldots,N\}$ and $t\in(t_{n-1},t_{n}]$, the following estimates hold:

[TABLE]

and

[TABLE]

Proof.

Recall the definitions of $G$ and ${\mathcal{G}}$ from (6.1) and (6.2). First, we add and subtract a term and then apply the triangle inequality. Then, for every $n\in\{1,\ldots,N\}$ and $t\in(t_{n-1},t_{n}]$ we arrive at

[TABLE]

by an application of Itō’s isometry. Furthermore, due to the Lipschitz continuity of $g$ we obtain

[TABLE]

where the last step follows from the identity

[TABLE]

which holds for every $s\in(t_{i-1},t_{i}]$, $i\in\{1,\ldots,N\}$. Finally, it follows from the same arguments as in the proof of Lemma 3.6 and by (4.1) for every $t\in(t_{n-1},t_{n}]$ that

[TABLE]

Together with the a priori bounds from Lemma 5.4 this shows (6.4).

It remains to prove the estimates (6.5) and (6.6). For (6.5) we first apply the Burkholder–Davis–Gundy-type inequality from Lemma 2.2 with constant $C_{\rho}$ and obtain for every $n\in\{1,\ldots,N\}$ and $t\in(t_{n-1},t_{n}]$ that

[TABLE]

where we also made use of the linear growth bound (4.1) in the last step. This proves (6.5). The bound in (6.6) can be shown by analogous arguments. ∎

The next lemma generalizes an important estimate from the proof of Theorem 3.7 to the multi-valued setting. In particular, we refer to Lemma 3.5 and (3.17).

Lemma 6.2**.**

Let Assumptions 4.1 to 4.4 be satisfied. For every step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, with $5L_{b}k\in[0,1)$, let the families $(X^{n})_{n\in\{0,\dots,N\}}$ and $(\eta^{n})_{n\in\{0,\dots,N\}}$ of random variables be as stated in Theorem 5.3. Then there exists $K_{\delta\eta}\in(0,\infty)$ independent of the step size $k$ such that

[TABLE]

Proof.

The nonnegativity follows immediately from the monotonicity of $f$. To prove the second inequality, we insert the scheme (5.3) and obtain

[TABLE]

For (6.7) we obtain

[TABLE]

because of the telescopic structure. Furthermore, it follows from Assumptions 4.1 and 4.4 that

[TABLE]

For (6.8) we apply Hölder’s inequality with $\rho=\max(2,p)$ and $\frac{1}{\rho}+\frac{1}{\rho^{\prime}}=1$ to obtain

[TABLE]

Then, from applications of the triangle inequality and Lemma 5.4, we get

[TABLE]

We apply the polynomial growth bound satisfied by $f$ and see that, for $p\in[2,\infty)$,

[TABLE]

is fulfilled, while for $p\in[1,2)$ we have

[TABLE]

In both cases the appearing terms are finite because of Assumption 4.4. Moreover, a further application of the triangle inequality yields

[TABLE]

Due to the linear growth bound (4.1) on $b$ and the a priori bound (5.4), it then follows that

[TABLE]

By application of Lemma 2.2 with constant $C_{\rho}$, we obtain

[TABLE]

Together with the linear growth bound (4.1) on $g$, this shows that

[TABLE]

Putting the estimates together proves the desired bound. ∎

We are now prepared to state and prove the main result of this section. While the main ingredients of the proof still consist of techniques introduced in [37, Sect. 4] for deterministic problems, the proof is somewhat more technical than the proof of Theorem 3.7. In particular, due to the presence of Lipschitz perturbations in the general problem (1.5) it is no longer possible to avoid an application of a Gronwall lemma. Moreover, as in [37, Sect. 4] we impose the following additional assumption on the multi-valued mapping $f$.

Assumption 6.3**.**

There exists $\gamma\in(0,\infty)$ such that, for every $v,w,z\in D(f)$, $f_{v}\in f(v)$, $f_{w}\in f(w)$, and $f_{z}\in f(z)$,

[TABLE]

In Lemma 3.2, we already proved that, if $f$ is the subdifferential of a convex potential, then Assumption 6.3 is satisfied with $\gamma=1$. For a further example, we refer to Section 7.

Theorem 6.4**.**

Let Assumptions 4.1 – 4.4 and Assumption 6.3 be satisfied. Let the step size $k=\frac{T}{N}$, $N\in{\mathbf{N}}$, be such that $8L_{b}k\in[0,1)$. Then there exists a constant $C\in(0,\infty)$ independent of $k$ such that

[TABLE]

Proof.

Let us first introduce some additional notation. We will denote the error between the exact solution $X$ to (1.5) and the numerical approximation ${\mathcal{X}}$ defined in (6.3) by $E(t):=X(t)-{\mathcal{X}}(t)$, $t\in[0,T]$. Furthermore, it will be convenient to split the error into two parts

[TABLE]

where we define

[TABLE]

${\mathbf{P}}$-almost surely for every $t\in(0,T]$. We expand the square of the norm of $E$ as

[TABLE]

In order to estimate the terms on the right-hand side of (6.11), we first observe in (6.9) that $E_{1}$ has absolutely continuous sample paths with $E_{1}(0)=0$. Hence we have $\frac{1}{2}\frac{\mathrm{d}}{\,\mathrm{d}t}|E_{1}(t)|^{2}=\langle\dot{E}_{1}(t),E_{1}(t)\rangle$ for almost every $t\in[0,T]$. Therefore, after integrating from [math] to $t\in(0,T]$ we get

[TABLE]

Furthermore, we also have

[TABLE]

Thus, after combining (6.12) and (6.13) we obtain

[TABLE]

For the first integral on the right-hand side of (6.14) we insert the derivative of $E_{1}$ and the definition of the error process $E$. This yields, for almost every $s\in(0,T]$,

[TABLE]

After recalling the definition of ${\mathcal{X}}$ we use Assumptions 4.1 and 6.3. Then, for almost every $s\in(t_{n-1},t_{n}]$ and all $n\in\{1,\ldots,N\}$, we get

[TABLE]

where the second term in the last step is non-positive due to the monotonicity of $f$ (cf. Definition 2.1). Moreover, because of the Lipschitz continuity of $b$, we have for almost every $s\in(0,T]$ that

[TABLE]

where we also made use of Young’s inequality. In addition, for every $n\in\{1,\ldots,N\}$ and $s\in(t_{n-1},t_{n}]$, we have that

[TABLE]

Therefore,

[TABLE]

Altogether, for every $t\in(t_{n-1},t_{n}]$ and $n\in\{1,\ldots,N\}$, we have shown that

[TABLE]

where we also inserted that $\int_{t_{n-1}}^{t}(t_{n}-s)\,\mathrm{d}s\leq\int_{t_{n-1}}^{t_{n}}(t_{n}-s)\,\mathrm{d}s=\frac{1}{2}k^{2}$ as well as $\int_{t_{n-1}}^{t}(t_{n}-s)^{2}\,\mathrm{d}s\leq\frac{1}{3}k^{3}$. It follows that, for every $n\in\{1,\ldots,N\}$ and $t\in(t_{n-1},t_{n}]$,

[TABLE]

Hence, together with Lemma 5.4 and Lemma 6.2 this shows that

[TABLE]

Next, we give an estimate for the second integral on the right-hand side of (6.14). For every $n\in\{1,\ldots,N\}$ and $t\in(t_{n-1},t_{n}]$ we decompose the integral as follows

[TABLE]

For every $i\in\{1,\ldots,n-1\}$ we then add and subtract $E_{2}(t_{i})$ in the second slot of the inner product in the first term on the right-hand side of (6.16). This gives

[TABLE]

After inserting the definition of $E_{2}$ from (6.10) the first integral is then equal to

[TABLE]

for all $i,n\in\{1,\ldots,N\}$, $i<n$, and $t\in(t_{n-1},t_{n}]$. Since $E_{1}(t_{i})-E_{1}(t_{i-1})=E(t_{i})-E(t_{i-1})-(E_{2}(t_{i})-E_{2}(t_{i-1}))$ is square-integrable and ${\mathcal{F}}_{t_{i}}$-measurable it therefore follows that

[TABLE]

for all $n\in\{1,\ldots,N\}$, $t\in(t_{n-1},t_{n}]$ and $t_{i}<t$. Hence, after taking expectations in (6.16) we arrive at

[TABLE]

Inserting the definitions (6.9) and (6.10) of $E_{1}$ and $E_{2}$ and applying Hölder’s inequality with $\rho=\max(2,p)$ and $\frac{1}{\rho}+\frac{1}{\rho^{\prime}}=1$, we get

[TABLE]

In the following, we will estimate $\Gamma_{1}$, $\Gamma_{2}$, $\Gamma_{3}$, and $\Gamma_{4}$ separately. For $\Gamma_{1}$ we obtain after an application of Hölder’s inequality for sums that

[TABLE]

If $p\in[2,\infty)$ then $\rho=p$ and $\rho^{\prime}=q$. In this case all integrals appearing are finite due to the bounds in Theorem 4.7 and Lemma 5.4. Moreover, if $p\in(1,2)$ then $\rho=\rho^{\prime}=2<q$. Then it follows from further applications of Hölder’s inequality and Jensen’s inequality that

[TABLE]

as well as

[TABLE]

Hence, we arrive at the same conclusion. If $p=1$ then the processes $(\eta(t))_{t\in[0,T]}$ and $(\eta^{n})_{n\in\{1,\ldots,N\}}$ are globally bounded due to the bound on $f$ in Assumption 4.1. Using Lemma 6.1 we see that

[TABLE]

Altogether, this yields

[TABLE]

for a suitable constant $C_{\Gamma}\in(0,\infty)$, which is independent of $k$. To estimate $\Gamma_{2}$ we argue analogously as in the case for $\Gamma_{1}$ to obtain that

[TABLE]

The first factor is bounded as we saw in the case for $\Gamma_{1}$. Furthermore, using Lemma 6.1, we have that

[TABLE]

Due to the a priori bound (5.4), it follows that there exists a constant $C_{\Gamma_{2}}\in(0,\infty)$, which does not depend on $k$ such that

[TABLE]

The estimates $\Gamma_{3}$ and $\Gamma_{4}$ follow analogously with the only new term that appears is of the form

[TABLE]

which is bounded due to Theorem 4.7 and the a priori bound (5.4). Therefore, there exist constants $C_{\Gamma_{3}},C_{\Gamma_{4}}\in(0,\infty)$ such that

[TABLE]

Hence, we obtain

[TABLE]

After taking expectations in (6.11) and inserting (6.14), (6.15), (6.17) as well as (6.4) from Lemma 6.1, we obtain for every $t\in(0,T]$ that

[TABLE]

The assertion then follows from an application of Gronwall’s lemma, see for example, [11, Appendix B]. ∎

Remark 6.5**.**

Up to this point, we only proved convergence for $X$ but not for $\eta$. However, from the existence of $X^{n}$ we also obtain that

[TABLE]

Analogously, we can write for the exact solution $\eta$ that

[TABLE]

Therefore, from the convergence of ${\mathcal{X}}$ to $X$ and the Lipschitz continuity of $b$ and $g$ we also obtain the estimate

[TABLE]

for every $n\in\{1,\ldots,N\}$.

7. Examples

7.1. Discontinuous drift coefficient

In this example, we show that Assumption 4.1 includes overdamped Langevin-type equations with a possibly discontinuous drift $f$. We consider the convex, nonnegative, yet not continuously differentiable function $\Phi(x):=|x|$, $x\in{\mathbf{R}}$, which has a multi-valued subdifferential $f\colon{\mathbf{R}}\to 2^{{\mathbf{R}}}$ defined by

[TABLE]

This mapping fulfills Assumption 4.1 for $p=1$. To be more precise, $f$ is a monotone function and there exists no proper monotone extension of its graph. In fact, the subdifferential of any proper, lower semi-continuous and convex function is a maximal monotone mapping by a well-known theorem of Rockafellar, cf. [45, Cor. 31.5.2] or [47, Satz 3.23].

Furthermore, we notice that $f_{x}x=\text{sgn}(x)x=|x|$ as well as $|f_{x}|\leq 1$ for every $x\in{\mathbf{R}}$ and $f_{x}\in f(x)$. This shows that $f$ fulfills all the conditions of Assumption 4.1. It remains to verify Assumption 6.3. Since $f$ is the subdifferential of $\Phi$ the variational inequality (3.2) is still satisfied in the sense that

[TABLE]

for all $x,y\in{\mathbf{R}}$ and $f_{x}\in f(x)$. Following the same steps as in the proof of Lemma 3.2 but replacing $f(v)$, $f(w)$, and $f(z)$ by arbitrary elements $f_{v}\in f(v)$, $f_{w}\in f(w)$, and $f_{z}\in f(z)$, respectively, shows that Assumption 6.3 is fulfilled. Therefore, the backward Euler–Maruyama method (1.6) is well-defined and yields an approximation of the exact solution $X$ of

[TABLE]

where $b\colon{\mathbf{R}}\to{\mathbf{R}}$ and $g\colon{\mathbf{R}}\to{\mathbf{R}}^{1,m}$ are Lipschitz continuous and $X_{0}\in L^{2}(\Omega)$. To be more precise, the piecewise linear interpolant ${\mathcal{X}}$ of the values $(X^{n})_{n\in\{0,\dots,N\}}$ defined in (6.3) fulfills

[TABLE]

for $C\in(0,\infty)$ that does not depend on the step size $k=\frac{T}{N}$. However, let us mention that the strong order of convergence of $1/4$ is not necessarily optimal in this particular example. We refer the reader to [9] for a corresponding result on the forward Euler–Maruyama method.

7.2. Stochastic $p$-Laplace equation

As a second example, we consider the discretization of the stochastic $p$-Laplace equation. A similar setting is studied in [5]. For a more detailed introduction to this class of problems, we refer the reader to this work and the references therein.

For $p\in[2,\infty)$ and $T\in(0,\infty)$ the stochastic $p$-Laplace equation is given by

[TABLE]

where ${\mathcal{D}}\subset{\mathbf{R}}^{n}$, $n\in{\mathbf{N}}$, is a bounded Lipschitz domain. By $W\colon[0,T]\times\Omega\to{\mathbf{R}}^{m}$, $m\in{\mathbf{N}}$, we denote a standard $({\mathcal{F}}_{t})_{t\geq 0}$-adapted Wiener process. We also assume that the initial value $u_{0}\colon{\mathcal{D}}\times\Omega\to{\mathbf{R}}$ fulfills

[TABLE]

Furthermore, let $\Psi\colon{\mathbf{R}}\to{\mathcal{L}}_{2}({\mathbf{R}}^{m};{\mathbf{R}})$ be a Lipschitz continuous mapping, where ${\mathcal{L}}_{2}({\mathbf{R}}^{m};{\mathbf{R}})$ denotes the space of Hilbert–Schmidt operators from ${\mathbf{R}}^{m}$ to ${\mathbf{R}}$. Note that the Nemytskii operator $\tilde{\Psi}\colon L^{2}({\mathcal{D}})\to{\mathcal{L}}_{2}({\mathbf{R}}^{m};L^{2}({\mathcal{D}}))$, given by $[\tilde{\Psi}(u)](x)=\Psi(u(x))$ for $u\in L^{2}({\mathcal{D}})$, is also Lipschitz continuous and will be of importance in the weak formulation below.

Further, let $W^{1,p}_{0}({\mathcal{D}})$ be the Sobolev space of weakly differentiable and $p$-fold integrable functions on ${\mathcal{D}}$ with vanishing trace on the boundary $\partial{\mathcal{D}}$, see [46, Section 1.2.3] or [39, Section 4.5] for a precise definition. The dual space of $W_{0}^{1,p}({\mathcal{D}})$ is denoted by $W^{-1,p}({\mathcal{D}})$ in the following. Then, the stochastic $p$-Laplace equation (7.1) has a solution $(u(t))_{t\in[0,T]}$ which is progressively measurable and an element of $L^{2}(\Omega;C([0,T];L^{2}({\mathcal{D}})))\cap L^{p}(\Omega;L^{p}(0,T;W_{0}^{1,p}({\mathcal{D}})))$. For further details we refer to [27, Example 4.1.9, Theorem 4.2.4].

For a spatial discretization of (7.1), we use a family of finite element spaces $(V_{h})_{h>0}$ such that $V_{h}\subset W_{0}^{1,p}({\mathcal{D}})$ for every $h>0$. Hereby, we interpret $h$ as a spatial refinement parameter. In the following, we consider a fixed parameter value $h>0$. By $d\in{\mathbf{N}}$ we then denote the dimension of the space $V_{h}$.

The spatially semi-discrete problem is to find a progressively measurable stochastic process $(u_{h}(t))_{t\in[0,T]}$ in the space $L^{2}(\Omega;C([0,T];L^{2}({\mathcal{D}})))\cap L^{p}(\Omega;L^{p}(0,T;V_{h}))$ such that

[TABLE]

for every $v_{h}\in V_{h}$ and $t\in[0,T]$. Hereby, $P_{h}\colon L^{2}({\mathcal{D}})\to V_{h}$ is the $L^{2}({\mathcal{D}})$-orthogonal projection onto $V_{h}$.

In order to apply our results from the previous sections, we rewrite (7.3) as a problem in ${\mathbf{R}}^{d}$. To this end, we consider a one-to-one relation between $V_{h}$ and ${\mathbf{R}}^{d}$ given by

[TABLE]

for a basis $\{\varphi_{1},\dots,\varphi_{d}\}$ of $V_{h}$. Through (7.4) we induce additional norms on ${\mathbf{R}}^{d}$ which are given by

[TABLE]

for every $x\in{\mathbf{R}}^{d}$. Observe that the norm $\|\cdot\|_{0}$ is also induced by the inner product

[TABLE]

where the mass matrix $M_{h}$ is symmetric and positive definite. Since all norms on ${\mathbf{R}}^{d}$ are equivalent, for each $i\in\{-1,0,1\}$ there exists $c_{i},C_{i}\in(0,\infty)$ such that

[TABLE]

for all $x\in{\mathbf{R}}^{d}$.

The $p$-Laplace operator in the spatially semi-discrete problem (7.3) can be written as $A_{h}\colon V_{h}\to V_{h}$ which is implicitly defined by

[TABLE]

for all $v_{h},w_{h}\in V_{h}$. By the same arguments as in [27, Example 4.1.9] one can easily verify that $A_{h}$ fulfills

[TABLE]

for all $v_{h},w_{h}\in V_{h}$. Then, for $x,y\in{\mathbf{R}}^{d}$ and associated $v_{x},v_{y}\in V_{h}$, we introduce mappings $\tilde{f}\colon{\mathbf{R}}^{d}\to{\mathbf{R}}^{d}$ and $\tilde{g}\colon{\mathbf{R}}^{d}\to{\mathbf{R}}^{d,m}$ implicitly by

[TABLE]

for $z\in{\mathbf{R}}^{m}$ and use these functions to define $f(x):=M_{h}\tilde{f}(x)$ as well as $g(x):=M_{h}^{\frac{1}{2}}\tilde{g}(x)$ for every $x\in{\mathbf{R}}^{d}$. As we assumed that $v_{x}\mapsto\tilde{\Psi}(v_{x})$ is Lipschitz continuous, there exists $L_{g}\in(0,\infty)$ such that

[TABLE]

for $x,y\in{\mathbf{R}}^{d}$ and $v_{x},v_{y}\in V_{h}$ fulfilling (7.4) and an orthonormal basis $\{e_{j}\}_{j\in\{1,\dots,m\}}$ of ${\mathbf{R}}^{m}$. Thus, $g$ fulfills Assumption 4.3. Due the integrability condition to (7.2) for $u_{0}$, it follows that $X_{0}$ fulfills Assumption 4.4.

Moreover, we see that $f$ is monotone, coercive, and bounded as we can write

[TABLE]

as well as

[TABLE]

and

[TABLE]

for all $x,y\in{\mathbf{R}}^{d}$ and $v_{x},v_{y}\in V_{h}$ fulfilling (7.4). Here, $\|\cdot\|_{{\mathcal{L}}({\mathbf{R}}^{m})}$ denotes the matrix norm in ${\mathbf{R}}^{m}$ which is induced by $|\cdot|$. Therefore, Assumption 4.1 is satisfied. To prove that $f$ fulfills Assumption 6.3 we note that the mapping $\Phi\colon V_{h}\to[0,\infty)$ given by

[TABLE]

is a potential of $A_{h}$, compare [46, Example 4.23]. Since $\Phi$ is convex it follows that

[TABLE]

where we use [13, Kapitel III, Lemma 4.10]. In the same way as in Lemma 3.2 we obtain that

[TABLE]

for all $v_{z},v_{x},v_{y}\in V_{h}$. Applying the definition of $f$, we then get

[TABLE]

for $x,y,z\in{\mathbf{R}}^{d}$ and $v_{x},v_{y},v_{z}\in V_{h}$ fulfilling (7.4). This shows that $f$ also fulfills Assumption 6.3.

Consequently, the results of the previous sections are applicable. More precisely, the backward Euler scheme (1.6) has a unique solution $(X^{n})_{n\in\{0,\dots,N\}}$ (cf. Theorem 5.3). Theorem 6.4 then states that the piecewise linear interpolant ${\mathcal{X}}$ of the values $(X^{n})_{n\in\{1,\dots,N\}}$ defined in (6.3) fulfills

[TABLE]

for $C\in(0,\infty)$ that does not depend on the step size $k$ where $X$ is the solution to the single-valued stochastic differential equation

[TABLE]

Observe that our proof does not yet rule out that the constant $C$ above depends on the dimension $d$ of the finite element space $V_{h}$. Hence, this is not a complete analysis of a full discretization of the stochastic partial differential equation (7.1) and a more detailed analysis is subject to future work. We refer to [5] for a related result in this direction.

Let us emphasize that, unlike the results in [5], we do not have to impose any temporal regularity assumption on the exact solution of (7.1) or on the solution of the semi-discrete problem (7.3). Since such regularity conditions are often not easily verified for quasi-linear stochastic partial differential equations we are confident that our approach could lead to interesting new insights in the numerical analysis of such infinite dimensional problems.

Acknowledgment

ME would like to thank the Berlin Mathematical School for the financial support. RK also gratefully acknowledges financial support by the German Research Foundation (DFG) through the research unit FOR 2402 – Rough paths, stochastic partial differential equations and related topics – at TU Berlin.