# Analysis of a positivity-preserving splitting scheme for some nonlinear
stochastic heat equations

Charles-Edouard Bréhier  Université de Pau et des Pays de l’Adour, E2S UPPA,
CNRS, LMAP, Pau, France<EMAIL_ADDRESS>,  David Cohen
Department of Mathematical Sciences, Chalmers University of Technology and
University of Gothenburg, 41296 Gothenburg, Sweden<EMAIL_ADDRESS>and  Johan Ulander  Department of Mathematical Sciences, Chalmers University
of Technology and University of Gothenburg, 41296 Gothenburg, Sweden
<EMAIL_ADDRESS>
###### Abstract.

We construct a positivity-preserving Lie–Trotter splitting scheme with finite
difference discretization in space for approximating the solutions to a class
of nonlinear stochastic heat equations with multiplicative space-time white
noise. We prove that this explicit numerical scheme converges in the mean-
square sense, with rate $1/4$ in time and rate $1/2$ in space, under
appropriate CFL conditions. Numerical experiments illustrate the superiority
of the proposed numerical scheme compared with standard numerical methods
which do not preserve positivity.

AMS Classification. 60H35. 60M15. 65J08.

Keywords. Stochastic partial differential equations. Stochastic heat equation.
Splitting scheme. Positivity-preserving scheme. Mean-square convergence.

##  1\. Introduction

Starting with the seminal work [38] on an implicit scheme for stochastic
quasi-linear parabolic partial differential equations in $1995$, the field of
numerical analysis of stochastic partial differential equations (SPDEs) has
gained a huge interest during the last decades. We refer the interested
readers to [79, 25, 46, 26, 53] for references on the theory of SPDEs and to
[36, 73, 27, 69, 42, 84, 68, 80, 61, 85, 54, 66, 45, 44, 47, 28, 4, 5, 82, 22,
55, 53, 48, 43, 33, 83, 81, 65, 64, 56, 2, 78, 52, 37, 50, 29, 17, 6] for
references on the numerical analysis of SPDEs (with a particular focus on
works related to strong convergence for parabolic SPDEs).

In this work we propose and study a novel positivity-preserving numerical
scheme for a fully discrete approximation of the following nonlinear
Stochastic Heat Equation (SHE) with multiplicative space-time white noise

(1) $\left\\{\begin{aligned}
&\partial_{t}u(t,x)=\partial_{xx}^{2}u(t,x)+g(u(t,x))\dot{W}(t,x),\\\
&u(t,0)=u(t,1)=0,\\\ &u(0,x)=u_{0}(x),\end{aligned}\right.$

for $(t,x)\in[0,T]\times[0,1]$ and where $u_{0}\geq 0$ is continuous,
$g\colon\mathbb{R}\to\mathbb{R}$ is globally Lipschitz continuous, of class
$\mathcal{C}^{1}$ and satisfies $g(0)=0$, and $\dot{W}$ is a space-time white
noise, see Section 2 for precise definitions and assumptions. Taking $g(x)=x$
in equation (1) results in the celebrated parabolic Anderson model, see for
instance [18]. This equation is used to model (particle) branching processes,
hydrodynamics with random forcing, and serves as a model for turbulent
diffusions.

The positivity-preserving property of the exact solutions to the SPDE (1) is
the subject of extensive research: two of the first results in this direction
can be found in [63, 74], where this property is proven to be true for noise
of the form $u^{\gamma}\dot{W}$ (where $1\leq\gamma<3/2$) and for a
nonlinearity that is of at most linear growth. The case of a Lipschitz
nonlinearity $g$ is studied in, for example, [26, 70, 62]. For the sake of
completeness, we mention the paper [8] on positivity of SHE with random
initial conditions, the paper [76] on problems with spatially homogeneous
Wiener process, the paper [20] on the stochastic fractional heat equation, the
paper [19] on problems in $\mathbb{R}^{n}$, as well as the paper [23] on
systems of SHEs with a spatially correlated noise. Note that these references
are considering the space domain to be $\mathbb{R}$ or $\mathbb{R}^{n}$. To
the best of our current knowledge, there are no corresponding results for the
case of compact domains with homogeneous Dirichlet boundary conditions.

While standard time integrators for SPDEs, such as the Euler–Maruyama scheme
[27], the semi-implicit Euler–Maruyama scheme [36], and the stochastic
exponential Euler integrator [55] do converge when applied to the problem (1),
they do not preserve the positivity property of the exact solution. Note that
the semi-implicit Euler scheme and the exponential Euler integrator preserve
positivity in the deterministic case ($g\equiv 0$ in equation (1)).

In this work, we employ a splitting strategy for the time integration of the
SPDE (1). This results in an efficient and positivity-preserving explicit time
integrator. In essence, a splitting integrator decomposes the vector field of
the original evolution equation in several parts, such that the arising
subsystems are exactly integrated (or easily). Splitting schemes have been
extensively studied and successfully applied to deterministic differential
equations, see for instance [39, 10, 60] and references therein. Splitting
schemes are also very popular for an efficient time discretization of
stochastic (partial) differential equations. We refer the reader to the
following non-exhaustive list of articles: [58, 21, 34, 51, 7, 30, 3, 24, 67,
15, 16, 9, 59, 13, 12, 14].

The preservation of positivity by numerical methods have been investigated in
several references in both the deterministic and stochastic settings. Without
being exhaustive, we mention the following articles on positivity-preserving
schemes for stochastic differential equations: [72, 75, 40, 57, 1, 71, 49,
41]. Finally, let us mention the recent reference [86] on a positivity-
preserving numerical scheme for the linear stochastic heat equation with
finite dimensional noise. We are not aware of works on the numerical analysis
of positivity-preserving schemes for SPDEs driven by space-time white noise.

The fully-discrete Lie–Trotter splitting scheme, see equation (14), considered
in this article combines a finite difference approximation in space and the
explicit recursion

$u_{m+1}^{{\rm LT}}=\exp\left(\tau N^{2}D^{N}\right)\hat{u}_{m+1}^{{\rm LT}},$

where for $n=1,\ldots,N-1$ one has

$\hat{u}_{m+1,n}^{{\rm LT}}=\exp\left(\sqrt{N}f(u_{m,n}^{{\rm
LT}})\Delta_{m,n}W-\frac{Nf(u_{m,n}^{{\rm
LT}})^{2}\tau}{2}\right)u_{m,n}^{{\rm LT}},$

where $\tau=T/M>0$ denotes the time step size, $h=1/N$ is the mesh size,
$\Delta_{m,n}W$ denote space-time Wiener increments, $N^{2}D^{N}$ the
$(N-1)\times(N-1)$ matrix of the discrete Laplace operator, and $g(v)=vf(v)$.
Observe that the diffusion part of (1) is solved exactly, while the noise part
is solved exactly in the case of the parabolic Anderson model (where one has
$g(v)=v$ and $f(v)=1$ and thus the subsystem is a geometric Brownian motion).
This shares similarity with the works [31, 77] on stochastic differential
equations. For a general mapping $g$, we freeze the factor $f$ at the previous
time point and obtain a geometric Brownian motion in the spirit of the
exponential scheme proposed in [11] for finite dimensional problems.

The main results of the paper are the following:

  * •

We obtain a fully discrete explicit approximation of the stochastic heat
equation (1) that is positivity-preserving, see Proposition 4.

  * •

We show bounds for the second moment of the numerical approximation under a
CFL condition $\tau/h={\rm O}(1)$ in Proposition 5.

  * •

We prove strong convergence, with rate $1/4$, for the temporal discretization
under a CFL condition $\tau/h^{2}={\rm O}(1)$, see Theorem 6. The strong
convergence of the fully discrete scheme is provided in Corollary 7.

We leave the study of weak convergence of the proposed scheme to possible
future works. On top of that, we show positivity of the exact solution to the
SPDE (1) on compact domains. This follows naturally from the numerical
analysis of the proposed approximation, see Proposition 2. Let us mention that
the CFL conditions above are not due to the discretization of the Laplace
operator, since the linear part is solved exactly. They are due to the
discretization of the contribution of the space-time white noise in the
temporal evolution. Numerical experiments confirm that the CFL condition is
necessary when studying the mean-square convergence of the proposed scheme.

This paper is organized as follows. Section 2 presents the setting,
assumptions, and useful results on the considered SHE. We also recall results
on the finite difference discretization from [35]. Section 3 contains the
definition of the proposed Lie–Trotter splitting as well as the main results
of the paper. We postpone their proofs to Section 5. We dedicate Section 4 to
numerical experiments illustrating our qualitative and quantitative results on
the proposed splitting scheme. The last section 6 briefly presents an
extension to systems of nonlinear stochastic heat equations. Appendix A
contains a proof of an auxiliary inequality used in the proofs of the main
results.

##  2\. Setting

This section provides the necessary setting for the description of the
considered class of nonlinear stochastic heat equations as well as of its
solution. We recall the notion of a mild solution and a standard well-
posedness result for completeness. In addition, we recall the spatial
discretization by finite difference from [35].

For any real-valued continuous function $v\colon[0,1]\to{\mathbb{R}}$, let
$\|v\|_{\infty}=\underset{x\in[0,1]}{\max}~{}|v(x)|$.

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space, equipped with a
filtration $\bigl{(}\mathcal{F}_{t}\bigr{)}_{t\geq 0}$ which satisfies the
usual conditions. The expectation operator is denoted by $\mathbb{E}[\cdot]$.
In the sequel, $C$ denotes a generic constant that may vary from line to line.
We sometimes use subscripts on $C$ to indicate dependence on parameters.

###  2.1. Description of the SPDE

Let us first introduce the main assumptions needed for the numerical analysis
of the proposed time integrator for the stochastic heat equation.

######  Assumption 1.

The initial value $u_{0}\colon[0,1]\to{\mathbb{R}}$ is a function of class
$\mathcal{C}^{3}$, and satisfies the conditions $u_{0}(0)=u_{0}(1)=0$.

Note that the regularity assumption on the initial value above is for ease of
presentation. For weaker conditions, see [35] or [2].

When discussing positivity-preserving properties, a further condition is
needed.

######  Assumption 2.

The initial value $u_{0}\colon[0,1]\to{\mathbb{R}}$ satisfies $u_{0}(x)\geq 0$
for all $x\in[0,1]$.

For the nonlinearity in the considered SPDE, we make use of the following.

######  Assumption 3.

The mapping $g\colon{\mathbb{R}}\to{\mathbb{R}}$ is of class
$\mathcal{C}^{1}$, is globally Lipschitz continuous, and satisfies $g(0)=0$.

We denote by ${\rm L_{g}}$ the Lipschitz constant of $g$:

${\rm L_{g}}=\underset{v_{1},v_{2}\in{\mathbb{R}},v_{2}\neq
v_{1}}{\sup}~{}\frac{|g(v_{2})-g(v_{1})|}{|v_{2}-v_{1}|}.$

The moment bounds and the error estimates presented below depend on the value
of the Lipschitz constant ${\rm L_{g}}$. This is not indicated in order to
simplify the notation.

We then introduce the auxiliary mapping $f\colon{\mathbb{R}}\to{\mathbb{R}}$
defined for all $v\in{\mathbb{R}}\setminus\\{0\\}$ by

(2) $f(v)=\frac{g(v)}{v}=\int_{0}^{1}g^{\prime}(rv)\,\text{d}r$

and by $f(0)=g^{\prime}(0)$. Since $g^{\prime}$ is continuous by Assumption 3,
the mapping $f$ is continuous and bounded, and one has the upper bound
$\underset{v\in{\mathbb{R}}}{\sup}~{}|f(v)|\leq{\rm L_{g}}$.

For a fixed time horizon $T>0$, let
$W=\left\\{W(t,x)\,\colon\,t\in[0,T],x\in[0,1]\right\\}$ be an
$\mathcal{F}_{t}$-adapted Brownian sheet. We recall that a Brownian sheet is a
Gaussian random field with mean zero and covariance
${\mathbb{E}}[W(t,x)W(s,y)]=(t\wedge s)(x\wedge y)$ for all $s,t\in[0,T]$ and
$x,y\in[0,1]$, see for instance [79]. We consider the stochastic heat equation
in the Itô sense

(3) $\left\\{\begin{aligned}
&\text{d}u(t,x)=\partial_{xx}^{2}u(t,x)\,\text{d}t+g(u(t,x))\,\text{d}W(t,x),\\\
&u(t,0)=u(t,1)=0,\\\ &u(0,x)=u_{0}(x)\end{aligned}\right.$

for $t\in[0,T]$ and $x\in[0,1]$, where $u_{0}$ and $g$ satisfy Assumption 1
and Assumption 3, respectively.

In order to define a mild solution of the stochastic heat equation (3), we
introduce the heat kernel

$G(t,x,y)=\sum_{j=1}^{\infty}e^{-j^{2}\pi^{2}t}\sin(j\pi x)\sin(j\pi y),$

for $t\geq 0,x,y\in[0,1]$, which is the fundamental solution of the
(deterministic) heat equation with homogeneous Dirichlet boundary conditions:

$\left\\{\begin{aligned}
&\text{d}v(t,x)=\partial_{xx}^{2}v(t,x)\,\text{d}t,\\\ &v(t,0)=v(t,1)=0,\\\
&v(0,x)=\delta(x),\end{aligned}\right.$

where the initial value is the Dirac delta function.

A mild solution to the SPDE (3) is a random field
$\bigl{(}u(t,x)\bigr{)}_{t\in[0,T],x\in[0,1]}$ satisfying the following
integral equation almost surely: for all $t\in[0,T]$ and $x\in[0,1]$, one has

(4)
$u(t,x)=\int_{0}^{1}G(t,x,y)u_{0}(y)\,\text{d}y+\int_{0}^{t}\int_{0}^{1}G(t-s,x,y)g(u(s,y))\,\text{d}W(s,y).$

The stochastic integral in (4) is understood in the Itô–Walsh sense, see for
instance [26, 46, 79].

We collect some properties of the mild solution $u(t,x)$ to the stochastic
heat equation (3) in the following statement, see for instance [35,
Proposition 3.7].

######  Proposition 1.

Consider the stochastic heat equation (3) under Assumptions 1 and 3. There
exists a unique mild solution $\left(u(t,x)\right)_{t\in[0,T],x\in[0,1]}$ to
the SPDE (3). In addition, for all $T\in(0,\infty)$, there exists
$C_{T}\in(0,\infty)$ such that

$\underset{t\in[0,T]}{\sup}~{}\underset{x\in[0,1]}{\sup}~{}{\mathbb{E}}[|u(t,x)|^{2}]\leq
C_{T}(1+\|u_{0}\|_{\infty}^{2}).$

Finally, the solution satisfies the following mean-square regularity property:
for all $T\in(0,\infty)$, there exists $C_{T}\in(0,\infty)$ such that for all
$x_{1},x_{2}\in[0,1]$ and all $t_{1},t_{2}\in[0,T]$ one has

(5)
$\left({\mathbb{E}}[|u(t_{2},x_{2})-u(t_{1},x_{1})|^{2}]\right)^{\frac{1}{2}}\leq
C_{T}\bigl{(}|t_{2}-t_{1}|^{\frac{1}{4}}+|x_{2}-x_{1}|^{\frac{1}{2}}\bigr{)}.$

In this article, our objective is to propose and analyze consistent numerical
schemes which preserve the following property of the exact solution: if the
initial value $u_{0}$ is nonnegative, then the exact solution to the
stochastic heat equation, $u(t,\cdot)$, remains nonnegative for all $t>0$.

######  Proposition 2.

Consider the stochastic heat equation (3) together with Assumptions 1, 2 and
3. Then, for all $t\in(0,\infty)$ and all $x\in[0,1]$, almost surely, one has

$u(t,x)\geq 0.$

The proof of Proposition 2 above is postponed to Section 5.5. It is a
consequence of the analysis of the fully-discrete numerical scheme and
combines two arguments: on the one hand, the numerical scheme satisfies a
variant of Proposition 2, see Proposition 4 below, on the other hand, Theorem
6 gives a strong convergence result of the numerical approximation. Note that
similar results are known when considering the stochastic heat equation on the
real line, see for instance the works [63, 74] and the lecture notes [70]. We
are not aware of positivity-preserving results for SPDEs on bounded domains.

###  2.2. Spatial discretization

Let us recall the spatial discretization based on a finite difference
approximation on a uniform grid from [35]. For any integer $N\in{\mathbb{N}}$,
let $h=1/N$ be the space mesh size, and let $x_{n}=nh$ for $0\leq n\leq N$ be
the grid points. Let $\kappa^{N}\colon[0,1]\to\\{x_{0},\ldots,x_{N}\\}$, be
the mapping defined by $\kappa^{N}(x)=x_{n}$ for $x\in[x_{n},x_{n+1})$ if
$n\in\\{0,\ldots,N-1\\}$, and $\kappa^{N}(1)=\kappa^{N}(x_{N})=x_{N}=1$.

Throughout this article, we use the convention that for any vector
$v=\bigl{(}v_{n}\bigr{)}_{1\leq n\leq N-1}\in{\mathbb{R}}^{N-1}$, we append
discrete homogeneous Dirichlet boundary conditions $v_{0}=0$ and $v_{N}=0$
when needed.

We discretize the initial value $u_{0}$ of the stochastic heat equation (3) by
$u_{0,n}^{N}=u_{n}^{N}(0)=u(0,x_{n})$ for $0\leq n\leq N$. Note that discrete
homogeneous Dirichlet boundary conditions $u_{0,0}^{N}=u_{0,N}^{N}=0$ are
satisfied owing to Assumption 1. Let us then define a piecewise linear
extension $u^{N}(0,\cdot)\colon[0,1]\to{\mathbb{R}}$ satisfying
$u^{N}(0,x_{n})=u_{0,n}^{N}$ for all $n=0,\ldots,N$, meaning that for
$x\in(0,1)$ one has

$u^{N}(0,x)=N\bigl{(}\kappa^{N}(x)+h-x\bigr{)}u(0,\kappa^{N}(x))+N\bigl{(}x-\kappa^{N}(x)\bigr{)}u(0,\kappa^{N}(x)+h).$

Let $D^{N}=\bigl{(}D_{ij}^{N}\bigr{)}_{1\leq i,j\leq N-1}$ denote the matrix
coming from a standard finite difference discretization of the Laplace
operator at the grid points $x_{n}$ with homogeneous Dirichlet boundary
conditions. The matrix $D^{N}$ is thus given by

$D^{N}=\begin{pmatrix}-2&1&0&\ldots&0&0&0\\\ 1&-2&1&\ddots&0&0&0\\\
0&1&-2&\ddots&0&0&0\\\ \vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&0&0&\ddots&-2&1&0\\\ 0&0&0&\ddots&1&-2&1\\\
0&0&0&\ldots&0&1&-2\end{pmatrix}.$

We then introduce the discrete heat kernel
$G^{N}(t)=\bigl{(}G_{ij}^{N}(t)\bigr{)}_{1\leq i,j\leq N-1}=e^{tN^{2}D^{N}}$,
for $t\geq 0$. By convention, set $G_{00}^{N}(t)=G_{NN}^{N}(t)=1$,
$G_{0N}^{N}(t)=G_{N0}^{N}(t)=0$ and $G_{0j}^{N}(t)=G_{Nj}^{N}(t)=0$ for all
$j\in\\{1,\ldots,N-1\\}$, in order to satisfy homogeneous discrete Dirichlet
boundary conditions. Finally, we extend the definition of
$G^{N}(t)=\bigl{(}G_{ij}^{N}(t)\bigr{)}_{1\leq i,j\leq N-1}$ to
$\bigl{(}G^{N}(t,x,y)\bigr{)}_{t\geq 0,x,y\in[0,1]}$ by asking that
$G^{N}(t,x_{i},y_{j})=NG_{ij}^{N}(t)$ for $0\leq i,j\leq N$ and for $t\geq 0$
and $y\in[0,1]$

$G^{N}(t,x,y)=N\bigl{(}\kappa^{N}(x)+h-x\bigr{)}G^{N}(t,\kappa^{N}(x),\kappa^{N}(y))+N\bigl{(}x-\kappa^{N}(x)\bigr{)}G^{N}(t,\kappa^{N}(x)+h,\kappa^{N}(y))$

for $x\in(0,1)$ and $G^{N}(t,0,y)=G^{N}(t,1,y)=0$. As a result, the mapping
$(x,y)\mapsto G^{N}(t,x,y)$ is piecewise linear in $x$ and piecewise constant
in $y$ at all times $t$.

It is worth recalling the following well-known property of the discrete heat
kernel: one has $G_{ij}^{N}(t)\geq 0$ for all $i,j\in\\{1,\ldots,N-1\\}$ and
$t\geq 0$. As a consequence, one has $G^{N}(t,x,y)\geq 0$ for $x,y\in[0,1]$
and $t\geq 0$. This property follows from the fact that $tN^{2}D^{N}$ is a
Metzler matrix, see for instance [32] for a definition, and the exponential of
a Metzler matrix has only non-negative elements.

We are now in position to define the spatial discretization $u^{N}$, for
$N\in{\mathbb{N}}$, by the following integral equality

(6)
$u^{N}(t,x)=\int_{0}^{1}G^{N}(t,x,y)u^{N}(0,\kappa^{N}(y))\,\text{d}y+\int_{0}^{t}\int_{0}^{1}G^{N}(t-s,x,y)g(u^{N}(s,\kappa^{N}(y)))\,\text{d}W(s,y)$

for $t\geq 0$ and $x\in[0,1]$. Note that the mapping $x\in[0,1]\mapsto
u^{N}(t,x)$ is linear on $[x_{n},x_{n+1}]$, for every
$n\in\\{0,\ldots,N-1\\}$, for every $t\geq 0$. In addition, one has
$u^{N}(t,0)=u^{N}(t,1)=0$ for every $t\geq 0$. For a practical implementation
of the scheme, it is sufficient to compute $u_{n}^{N}(t)=u^{N}(t,x_{n})$ for
all $1\leq n\leq N-1$. This is performed as follows: for all $t\geq 0$ and
$1\leq n\leq N-1$ one has

$u_{n}^{N}(t)=\sum_{j=1}^{N-1}G_{nj}^{N}(t)u_{0,j}^{N}+\sqrt{N}\sum_{j=1}^{N-1}\int_{0}^{t}G_{nj}^{N}(t-s)g(u_{j}^{N}(s))\,\text{d}W_{j}^{N}(s),$

where

$W_{n}^{N}(t)=\sqrt{N}\bigl{(}W(t,x_{n+1})-W(t,x_{n})\bigr{)}.$

By definition of a Wiener sheet, observe that the processes
$\bigl{(}W_{1}^{N}(t)\bigr{)}_{t\geq
0},\ldots,\bigl{(}W_{N-1}^{N}(t)\bigr{)}_{t\geq 0}$ are independent standard
real-valued Wiener processes, for any $N\in{\mathbb{N}}$.

Introduce the ${\mathbb{R}}^{N-1}$-valued process $u^{N}$ defined by
$u^{N}(t)=\bigl{(}u_{n}^{N}(t)\bigr{)}_{1\leq n\leq N-1}$ for all $t\geq 0$.
This process is solution of the following stochastic differential equation

(7)
$\text{d}u^{N}(t)=N^{2}D^{N}u^{N}(t)\,\text{d}t+\sqrt{N}g(u^{N}(t))\,\text{d}W^{N}(t)$

with initial value $u^{N}(0)=\bigl{(}u_{0}^{N}\bigr{)}_{1\leq n\leq N-1}$,
where the notation
$\bigl{(}g(u^{N}(t))\,\text{d}W^{N}(t)\bigr{)}_{n}=g(u_{n}^{N}(t))\,\text{d}W_{n}^{N}(t)$
is used.

Let us recall the following convergence result for the spatial discretization,
see [35, Theorem 3.1].

######  Proposition 3.

Consider the stochastic heat equation (3) with a nonlinearity $g$ satisfying
Assumption 3. Denote by $\left(u(t,x)\right)_{t\in[0,T],x\in[0,1]}$ its exact
solution and by $\left(u^{N}(t,x)\right)_{t\in[0,T],x\in[0,1]}$ the numerical
approximation by finite differences with mesh size $h=1/N$. For all
$T\in(0,\infty)$ and any initial value $u_{0}$ satisfying Assumption 1, there
exists $C_{T}(u_{0})\in(0,\infty)$ such that for all $h=1/N$ with
$N\in{\mathbb{N}}$ one has

(8)
$\underset{t\in[0,T]}{\sup}~{}\underset{x\in[0,1]}{\sup}~{}\bigl{(}{\mathbb{E}}[|u^{N}(t,x)-u(t,x)|^{2}]\bigr{)}^{\frac{1}{2}}\leq
C_{T}(u_{0})h^{\frac{1}{2}}.$

In the error analysis below, the following auxiliary result from [2] (see
Proposition $2.4$) on the temporal regularity of $u^{N}$ is used: there exists
$C_{T}(u_{0})\in(0,\infty)$ such that for all $t,s\in[0,T]$, one has

(9)
$\underset{N\in{\mathbb{N}}}{\sup}~{}\underset{x\in[0,1]}{\sup}~{}{\mathbb{E}}[|u^{N}(t,x)-u^{N}(s,x)|^{2}]\leq
C_{T}(u_{0})|t-s|^{\frac{1}{2}}.$

##  3\. The positivity-preserving splitting scheme

In the core part of this paper, we present and study the strong convergence of
an efficient and positivity-preserving time integrator for the stochastic heat
equation (3).

Let $T\in(0,\infty)$ and divide the interval $[0,T]$ into $M\in{\mathbb{N}}$
subintervals $[t_{m},t_{m+1}]$ of length $\tau=T/M$, where $t_{m}=m\tau$ for
$m\in\\{0,\ldots,M\\}$. Introduce the mapping
$\ell^{M}\colon[0,T]\to\\{t_{0},\ldots,t_{M}\\}$, defined by
$\ell^{M}(t)=t_{m}$ for all $t\in[t_{m},t_{m+1})$, if
$m\in\\{0,\ldots,M-1\\}$, and $\ell^{M}(T)=\ell^{M}(t_{M})=t_{M}=T$.

We propose a fully-discrete explicit scheme based on a Lie–Trotter splitting
strategy producing approximations $u_{m}^{{\rm LT}}=\bigl{(}u_{m,n}^{{\rm
LT}}\bigr{)}_{1\leq n\leq N-1}$ of the finite difference approximation
$u^{N}(t_{m})=\bigl{(}u_{n}^{N}(t_{m})\bigr{)}_{1\leq n\leq N-1}$ at the grid
times $t_{m}$, $m=0,\ldots,M$. We set the initial value to be $u_{0,n}^{{\rm
LT}}=u_{n}^{N}(0)=u_{0,n}^{N}$ for all $1\leq n\leq N-1$. As above, one has
$u_{m,0}^{{\rm LT}}=0$ and $u_{m,N}^{{\rm LT}}=0$ for all
$m\in\\{0,\ldots,M\\}$. In this way, homogeneous Dirichlet boundary conditions
are satisfied by the numerical scheme at all times.

We explain the construction of the scheme in Section 3.1. We then describe the
main results of this article: the positivity-preserving property of the
splitting scheme (Proposition 4) and the mean-square convergence in time with
order $1/4$ (Theorem 6 and Corollary 7).

###  3.1. Description of the time integrator

Let us describe how the splitting scheme is constructed. Given the numerical
solution $u_{m}^{{\rm LT}}=\left(u_{m,n}^{{\rm LT}}\right)_{1\leq n\leq N-1}$
at grid time $t_{m}=m\tau$ for $0\leq m\leq M-1$, the solution $u_{m+1}^{{\rm
LT}}$ at the next grid time $t_{m+1}=t_{m}+\tau$ is constructed by
successively solving two subsystems in ${\mathbb{R}}^{N-1}$:

  * •

first, the linear Itô SDE system

(10) $\text{d}v_{m,n}^{M,N,1}(t)=\sqrt{N}v_{m,n}^{M,N,1}(t)f(u_{m,n}^{{\rm
LT}})\,\text{d}W_{n}^{N}(t),$

for $n\in\\{1,\ldots,N-1\\}$ and $t\in[t_{m},t_{m+1}]$, with initial value
$v_{m,n}^{M,N,1}(t_{m})=u_{m,n}^{{\rm LT}}$, where we recall (see equation (2)
in Section 2) that the auxiliary function $f$ is such that $g(v)=vf(v)$ for
all $v\in{\mathbb{R}}$;

  * •

second, the linear ODE system

(11) $\text{d}v_{m}^{M,N,2}(t)=N^{2}D^{N}v_{m}^{M,N,2}(t)\,\text{d}t,$

for $t\in[t_{m},t_{m+1}]$, with initial value
$v_{m,n}^{M,N,2}(t_{m})=v_{m,n}^{M,N,1}(t_{m+1})$.

Observe that the solutions of the two subsystems above are known: the solution
of the SDE (10) is given by

(12) $v_{m,n}^{M,N,1}(t)=\exp\left(\sqrt{N}f(u_{m,n}^{{\rm
LT}})\bigl{(}W_{n}^{N}(t)-W_{n}^{N}(t_{m})\bigr{)}-\frac{Nf(u_{m,n}^{{\rm
LT}})^{2}(t-t_{m})}{2}\right)u_{m,n}^{{\rm LT}},$

for all $t\in[t_{m},t_{m+1}]$, and the solution of the ODE (11) is given by

(13) $v_{m}^{M,N,2}(t)=e^{(t-t_{m})N^{2}D^{N}}v_{m}^{M,N,1}(t_{m+1}),$

for all $t\in[t_{m},t_{m+1}]$.

Gathering the expressions above gives the following expression for the
proposed Lie–Trotter splitting scheme

(14) $u_{m+1}^{{\rm LT}}=e^{\tau
N^{2}D^{N}}\left(\exp\Bigl{(}\sqrt{N}f(u_{m,n}^{{\rm
LT}})\Delta_{m,n}W-\frac{Nf(u_{m,n}^{{\rm
LT}})^{2}\tau}{2}\Bigr{)}u_{m,n}^{{\rm LT}}\right)_{1\leq n\leq N-1},$

where $\Delta_{m,n}W=W_{n}^{N}(t_{m+1})-W_{n}^{N}(t_{m})$. Observe that the
random variables $\bigl{(}\Delta W_{m,n}\bigr{)}_{0\leq m\leq M-1,1\leq n\leq
N-1}$ are independent standard real-valued Gaussian random variables.

The splitting scheme formula (14) can also be written as

$u_{m+1,n}^{{\rm
LT}}=\sum_{k=1}^{N-1}G_{nk}^{N}(\tau)\exp\Bigl{(}\sqrt{N}f(u_{m,k}^{{\rm
LT}})\Delta_{m,n}W-\frac{Nf(u_{m,k}^{{\rm
LT}})^{2}\tau}{2}\Bigr{)}u_{m,k}^{{\rm LT}},$

for all $m\in\\{0,\ldots,M-1\\}$ and $n\in\\{1,\ldots,N-1\\}$.

One of the key properties of the proposed splitting scheme is the following:
if the initial value $u_{0}^{{\rm LT}}=\bigl{(}u_{0,n}^{{\rm
LT}}\bigr{)}_{1\leq n\leq N-1}$ only has nonnegative elements, then for all
$m\in\\{1,\ldots,M\\}$ the numerical solution $u_{m}^{{\rm
LT}}=\bigl{(}u_{m,n}^{{\rm LT}}\bigr{)}_{1\leq n\leq N-1}$ at time
$t_{m}=m\tau$ also only has nonnegative elements almost surely. In other
words, the proposed scheme is positivity-preserving. This is stated in the
next proposition.

######  Proposition 4.

Let $M\in{\mathbb{N}}$ and $N\in{\mathbb{N}}$ be arbitrary integers and let
$T\in(0,\infty)$. Let Assumptions 1, 2 and 3 be satisfied. Let the sequence
$u_{0}^{{\rm LT}},\ldots,u_{M}^{{\rm LT}}$ be given by the splitting scheme
(14), with $h=1/N$ and $\tau=T/M$, with initial value $u_{0,n}^{{\rm
LT}}=u_{0}(x_{n})\geq 0$ for all $n\in\\{1,\ldots,N\\}$. Then, almost surely,
one has

$u_{m,n}^{{\rm LT}}\geq 0,$

for all $m\in\\{1,\ldots,M\\}$ and $n\in\\{1,\ldots,N-1\\}$.

###### Proof of Proposition 4.

The proof proceeds by recursion on the time index $m$.

  * •

Note that $u_{0,n}^{{\rm LT}}=u(0,x_{n})\geq 0$ for all
$n\in\\{0,\ldots,N\\}$.

  * •

Assume that the property $u_{m,n}^{{\rm LT}}\geq 0$, for all
$n\in\\{1,\ldots,N-1\\}$, holds at time $t_{m}=m\tau$. We prove that under
this assumption, it also holds at time $t_{m+1}=(m+1)\tau$.

The argument is straightforward: the solutions of the subsystems (10) and (11)
are nonnegative at all times when they have nonnegative initial values. More
precisely, first one has

$v_{m,n}^{M,N,1}(t_{m+1})=e^{\sqrt{N}f(u_{m,n}^{{\rm
LT}})\Delta_{m,n}W-N\frac{f(u_{m,n}^{{\rm LT}})^{2}\tau}{2}}u_{m,n}^{{\rm
LT}}\geq 0$

for all $n\in\\{1,\ldots,N-1\\}$. Second, using the inequality
$G_{nk}^{N}(\tau)\geq 0$ (see Section 2.2), one has

$u_{m+1,n}^{{\rm
LT}}=v_{m,n}^{M,N,2}(t_{m+1})=\sum_{k=1}^{N-1}G_{nk}^{N}(\tau)v_{m,k}^{M,N,1}(t_{m+1})\geq
0.$

Thus the positivity property of the numerical solution holds at time
$t_{m+1}=(m+1)\tau$.

As a consequence, the property $u_{m,n}^{{\rm LT}}\geq 0$, for all
$n\in\\{1,\ldots,N-1\\}$, holds for any $m\in\\{0,\ldots,M\\}$. The proof of
Proposition 4 is completed. ∎

###  3.2. Convergence results

Let us now prove that the proposed numerical scheme provides accurate
approximation of the exact solution. In this article, we show mean-square
error estimates and give orders of convergence with respect to $\tau=T/M$ and
$h=1/N$.

We impose a CFL stability condition in the sequel to ensure stability and
convergence of the Lie–Trotter splitting scheme (14) when applied to the
stochastic heat equation (3); more precisely, we introduce conditions of the
type $\tau\leq\gamma h$ or $\tau\leq\gamma h^{2}$ in the statements below, for
some (nonrandom) arbitrary parameter $\gamma\in(0,\infty)$. The conditions on
$\tau$ and $h$ above are equivalent to the conditions $\gamma M\geq TN$ and
$\gamma M\geq TN^{2}$ on $M$ and $N$ respectively.

Owing to Proposition 3, it is sufficient to focus on the error $u_{m,n}^{{\rm
LT}}-u_{n}^{N}(t_{m})$ to obtain estimates for the total error $u_{m,n}^{{\rm
LT}}-u(t_{m},x_{n})$. Proposition 5 shows moment bounds of the numerical
solution and is used to prove our main result in Theorem 6. As a corollary we
obtain convergence of $u^{{\rm LT}}_{m,n}$ to the exact solution
$u(t_{m},x_{n})$ at the grid points using results from [35].

Note that Assumption 2 on the positivity of the initial value is not needed in
the statements on the moment bounds and on the convergence of the scheme
below.

######  Proposition 5.

Assume that Assumptions 1 and 3 are satisfied. Let the sequence $u_{0}^{{\rm
LT}},\ldots,u_{M}^{{\rm LT}}$ be given by the Lie–Trotter splitting scheme
(14).

For all $\gamma\in(0,\infty)$ and all $T\in(0,\infty)$, there exists
$C_{\gamma,T}\in(0,\infty)$ such that for all $\tau=T/M$ and $h=1/N$
satisfying the condition $\tau\leq\gamma h$, one has

(15) $\underset{0\leq m\leq M}{\sup}~{}\underset{1\leq n\leq
N-1}{\sup}~{}{\mathbb{E}}\left[|u_{m,n}^{{\rm LT}}|^{2}\right]\leq
C_{\gamma,T}\bigl{(}1+\|u_{0}\|_{\infty}^{2}\bigr{)}.$

The proof of this proposition also provides moment bounds for a space-time
continuous version $u^{{\rm LT}}(t,x)$, defined by equation (21), of the
Lie–Trotter splitting scheme (14).

We are now in position to state the main convergence result of this article.
For ease of presentation, we only consider errors at space-time grid points.

######  Theorem 6.

Assume that Assumptions 1 and 3 are satisfied. Let the sequence $u_{0}^{{\rm
LT}},\ldots,u_{M}^{{\rm LT}}$ be given by the Lie–Trotter scheme (14), and let
$\bigl{(}u^{N}(t)\bigr{)}_{t\geq 0,0\leq n\leq N}$ be given by the spatial
semi-discretization scheme (7).

For all $\gamma\in(0,\infty)$ and $T\in(0,\infty)$, there exists
$C_{\gamma,T}(u_{0})\in(0,\infty)$ such that for all $\tau=T/M$ and $h=1/N$
satisfying the condition $\tau\leq\gamma h$, one has

(16) $\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{m,n}^{{\rm
LT}}-u_{n}^{N}(t_{m})|^{2}]\right)^{\frac{1}{2}}\leq
C_{\gamma,T}(u_{0})\left(\tau^{\frac{1}{4}}+\left(\frac{\tau}{h}\right)^{\frac{1}{2}}\right).$

In addition, for all $\tau=T/M$ and $h=1/N$ satisfying the condition
$\tau\leq\gamma h^{2}$, one has

(17) $\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{m,n}^{{\rm
LT}}-u_{n}^{N}(t_{m})|^{2}]\right)^{\frac{1}{2}}\leq
C_{\gamma,T}(u_{0})\tau^{\frac{1}{4}}.$

Proving (17) from the error estimate (16) under the stronger condition
$\tau\leq\gamma h^{2}$ is straightforward.

Combining Theorem 6 and Proposition 3, one directly obtains error estimates
for the fully-discrete scheme.

######  Corollary 7.

Consider the setting and assumptions of Theorem 6. For all
$\gamma\in(0,\infty)$ and $T\in(0,\infty)$, there exists
$C_{\gamma,T}(u_{0})\in(0,\infty)$ such that for all $\tau=T/M$ and $h=1/N$
satisfying the condition $\tau\leq\gamma h^{2}$, one has

(18) $\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{m,n}^{{\rm
LT}}-u(t_{m},x_{n})|^{2}]\right)^{\frac{1}{2}}\leq
C_{\gamma,T}(u_{0})h^{\frac{1}{2}}.$

We postpone the proofs of the above results to Section 5.

##  4\. Numerical experiments

In this section we provide numerical experiments to support and verify the
above theoretical results. Recall that $\tau=T/M>0$ is the time step size and
$h=1/N>0$ is the space mesh size. We compare the proposed Lie–Trotter
splitting scheme (14), denoted LT below, to the following classical time
integrators when applied to the spatially discretized system (7):

  * •

the Euler–Maruyama scheme (denoted EM below), see for instance [27]

$u^{\rm EM}_{m+1}=u^{\rm EM}_{m}+\tau N^{2}D^{N}u^{\rm
EM}_{m}+\sqrt{N}g(u^{\rm EM}_{m})\Delta_{m}W,$

  * •

the semi-implicit Euler–Maruyama scheme (denoted SEM below), see for instance
[36]

$u^{\rm SEM}_{m+1}=u^{\rm SEM}_{m}+\tau N^{2}D^{N}u^{\rm
SEM}_{m+1}+\sqrt{N}g(u^{\rm SEM}_{m})\Delta_{m}W,$

  * •

the stochastic exponential Euler integrator (denoted SEXP below), see for
instance [55]

$u^{\rm SEXP}_{m+1}=e^{\tau N^{2}D^{N}}\left(u^{\rm SEXP}_{m}+\sqrt{N}g(u^{\rm
SEXP}_{m})\Delta_{m}W\right).$

###  4.1. Preservation of the positivity

We start by illustrating the positivity-preserving property of the Lie–Trotter
scheme (LT) and show the lack of positivity-preserving behavior for the
Euler–Maruyama scheme (EM), the semi-implicit Euler–Maruyama scheme (SEM), and
the stochastic exponential scheme (SEXP). To do this, we use the same noise
samples for all time integrators when applied to the space-discretization of
the SPDE (3) as described in Section 2.2 with the initial condition
$u_{0}=\sin(\pi x)$ and final time $T=20$. We consider this problem with the
three choices of multiplicative term given by $g(v)=\lambda v$,
$g(v)=\lambda\ln\left(1+v\right)$, and $g(v)=\lambda\left(v+\sin(v)\right)$.
The real-valued parameter $\lambda>0$ is introduced to avoid the need to run
numerical experiments with very long time horizons $T$ in order to obtain
negative values for the numerical schemes SEXP and SEM. We remark that
$g(v)=\lambda\ln\left(1+v\right)$ is well-behaved for $v\geq 0$ but problems
may occur if $v\leq-1$. Since the proposed LT scheme is guaranteed to preserve
positivity, this is not problematic. However, this could happen for the time
integrators SEXP, SEM or EM. The numerical results are presented in Tables 1
and 2, where the notation $k/50$ indicates that $k$ out of $50$ samples remain
positive.

In Table 1, we let $g(v)=2.5v$ and we consider $50$ sample paths for each of
the time integrators for several choices of the discretization parameters
$\tau$ and $h$. Table 1 confirms that the LT scheme preserves positivity. This
is not the case for SEXP, SEM and EM. We observe that fewer samples of SEXP
and SEM contain negative values for small time steps $\tau$. This is expected
as each of the time integrators SEXP, SEM, and even EM, converges (for every
fixed $h$) to the exact, everywhere positive, solution of the space-
discretized system of SDEs in equation (7).

$(\tau,h)$ | LT | SEXP | SEM | EM
---|---|---|---|---
$(10^{-3},10^{-2})$ | $50/50$ | $0/50$ | $0/50$ | $0/50$
$(10^{-4},10^{-3})$ | $50/50$ | $50/50$ | $50/50$ | $0/50$
$(10^{-5},10^{-3})$ | $50/50$ | $50/50$ | $50/50$ | $0/50$
Table 1. Proportion of samples containing only positive values out of $50$
simulated sample paths for the Lie–Trotter splitting scheme (LT), the
stochastic exponential Euler integrator (SEXP), the semi-implicit
Euler–Maruyama scheme (SEM), and Euler–Maruyama scheme (EM) for the diffusion
coefficient $g(v)=2.5v$ and several choices of discretization parameters
$\tau$ and $h$.

In Table 2 we instead fix the discretization parameters $\tau=10^{-5}$ and
$h=10^{-3}$ and consider different types of multiplicative terms $g(v)$. We
again use $50$ samples in each of the entries of Table 2. From the results of
Table 2, one can observe the poor performance of the EM scheme in all cases.
This table also illustrates the fact that increasing the size of the
multiplicative term prevents SEM and SEXP to remain positive. It should be
clear that increasing the value of $\lambda$ even more, or the length of the
time interval, would hinder the numerical solutions to stay positive for all
time integrators except for the proposed Lie–Trotter splitting scheme.

$g(v)$ | LT | SEXP | SEM | EM
---|---|---|---|---
$2.5\ln(1+v)$ | $50/50$ | $50/50$ | $50/50$ | $0/50$
$3.5\ln(1+v)$ | $50/50$ | $50/50$ | $50/50$ | $0/50$
$5\ln(1+v)$ | $50/50$ | $47/50$ | $26/50$ | $0/50$
$2.5v$ | $50/50$ | $50/50$ | $50/50$ | $0/50$
$3.5v$ | $50/50$ | $50/50$ | $50/50$ | $0/50$
$5v$ | $50/50$ | $4/50$ | $50/50$ | $0/50$
$2.5\left(v+\sin(v)\right)$ | $50/50$ | $44/50$ | $50/50$ | $0/50$
$3.5\left(v+\sin(v)\right)$ | $50/50$ | $0/50$ | $0/50$ | $0/50$
$5\left(v+\sin(v)\right)$ | $50/50$ | $0/50$ | $0/50$ | $0/50$
Table 2. Proportion of samples containing only positive values out of $50$
simulated sample paths for the Lie–Trotter splitting scheme (LT), the
stochastic exponential Euler integrator (SEXP), the semi-implicit
Euler–Maruyama scheme (SEM), and the Euler–Maruyama scheme (EM) for several
choices of diffusion terms $g(v)$. The discretization parameters are
$\tau=10^{-5}$ and $h=10^{-3}$.

###  4.2. Mean-square errors

For the next numerical experiment, we discretize the stochastic heat equation
(3) with initial value $u_{0}(x)=\sin(\pi x)$ by a finite-difference scheme in
space with mesh size $h=2^{-8}$. The resulting system of stochastic
differential equations (7) is then discretized by the time integrators LT,
SEXP, and SEM. The classical EM scheme is not appropriate in this setting and
numerical results are thus not presented. The following choices for the
function $g$ are considered: $g(v)=v$ and $g(v)=\frac{v}{1+v^{2}}$ and
$g(v)=\ln(1+v)$, for $v\geq 0$, and $g(v)=v\exp(-v^{2})$. Figure 1 displays,
in a loglog plot, the mean-square errors

$\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\bigl{(}{\mathbb{E}}[|u_{m,n}^{\text{num}}-u^{\text{ref}}(t_{m},x_{n})|^{2}]\bigr{)}^{\frac{1}{2}}$

measured at the space-time grid points $(t_{m},x_{n})$ for the time interval
$[0,0.5]$. The reference solution $u^{\text{ref}}$ is computed using the LT
splitting scheme with time step size $\tau_{\text{ref}}=2^{-16}$. Here, $200$
samples have been used to approximate the expectations. We have checked that
the Monte Carlo error is negligible to observe mean-square convergence. In
this figure, one can observe a rate of convergence $1/2$ instead of $1/4$ in
the mean-square error estimates (17) for the splitting scheme in Theorem 6.
This is related to the mean-square error estimates (16) and the role of the
CFL condition $\tau\leq\gamma h^{2}$ to obtain (17).

(a) $g(v)=v$

(b) $g(v)=\frac{v}{(1+v^{2})}$

(c) $g(v)=\ln(1+v)$

(d) $g(v)=v\exp(-v^{2})$

Figure 1. Mean-square errors on the time interval $[0,0.5]$ of the splitting
scheme (LT), the stochastic exponential Euler integrator (SEXP), and the semi-
implicit Euler–Maruyama scheme (SEM). Mesh size $h=2^{-8}$ and average over
$200$ samples.

To illustrate this, we compute the mean-square errors of the Lie–Trotter
splitting scheme when applied to the finite difference discretization of the
stochastic heat equation with different values of the mesh size, namely
$h=2^{-4},2^{-6},2^{-8},2^{-10}$. This is presented only for the two
nonlinearities $g(v)=1.5v$ and $g(v)=1.5\frac{v}{1+v^{2}}$. We have used $200$
samples to approximate the expectations. The other parameters are the same as
in the previous numerical experiments. The results are presented in Figure 2.
In these experiments we observe upper bounds which are not uniform with
respect to $h$, in fact we observe the contribution of the error term
$\tau^{\frac{1}{2}}h^{-\frac{1}{2}}$ in the mean-square error estimates (16).

(a) $g(v)=v$

(b) $g(v)=\frac{v}{(1+v^{2})}$

Figure 2. Mean-square errors on the time interval $[0,0.5]$ of the splitting
scheme for several values of the spatial mesh $h$. Average over $200$ samples.

In the final numerical experiment, we consider the same parameters as above
and the function $g(v)=v^{1.25}$. Observe that this nonlinearity is not
globally Lipschitz continuous and is thus not covered by the the results from
Section 3.2. A convergence plot for the splitting scheme (14) is provided in
Figure 3. As above, we observe a mean-square order of convergence $1/2$, but
which should not be uniform with respect to $h$, similarly to what is observed
in Figure 2. To prove such rate of convergence is beyond the scope of this
paper and will be the subject of a future work.

Figure 3. Mean-square errors on the time interval $[0,0.5]$ of the splitting
scheme (LT) when applied to the stochastic heat equation (3) with
$g(v)=v^{1.25}$. Mesh size $h=2^{-8}$ and average over $200$ samples.

##  5\. Proofs of the main results

The objective of this section is to provide the proof of the results stated in
Section 3.2, namely the moment bounds in Proposition 5 and the mean-square
error estimates in Theorem 6 and in Corollary 7. We also prove Proposition 2,
which ensures positivity of the exact solution. Preliminary auxiliary tools
are given in Sections 5.1 and 5.2, before proceeding with the detailed proofs.

###  5.1. Auxiliary process

In this section, for any $M\in{\mathbb{N}}$ and $N\in{\mathbb{N}}$, we define
an auxiliary stochastic process $\bigl{(}u^{{\rm
LT}}(t,x)\bigr{)}_{t\in[0,T],x\in[0,1]}$ satisfying $u^{{\rm
LT}}(t_{m},x_{n})=u_{m,n}^{{\rm LT}}$ for all $m\in\\{0,\ldots,M\\}$ and
$n\in\\{1,\ldots,N-1\\}$. The auxiliary process $u^{{\rm LT}}$ is piecewise
continuous with respect to the spatial variable $x$, while its temporal
evolution on each interval $(t_{m},t_{m+1})$ follows a stochastic differential
equation similar to (10).

Recall that the auxiliary mappings
$\kappa^{N}\colon[0,1]\to\\{x_{0},\ldots,x_{N}\\}$ and
$\ell^{M}\colon[0,T]\to\\{t_{0},\ldots,t_{M}\\}$ are defined in Sections 2.2
and 4 respectively.

Let $n\in\\{1,\ldots,N-1\\}$ and $m\in\\{0,\ldots,M-1\\}$, then for all
$t\in[t_{m},t_{m+1}]$ set

(19) $u_{m,n}^{{\rm
LT}}(t)=\sum_{k=1}^{N-1}G_{nk}^{N}(\ell^{M}(t)-t_{m})v_{m,k}^{M,N,1}(t),$

where $v_{m,n}^{M,N,1}(t)=\exp\left(\sqrt{N}f(u_{m,n}^{{\rm
LT}})\bigl{(}W_{n}^{N}(t)-W_{n}^{N}(t_{m})\bigr{)}-\frac{Nf(u_{m,n}^{{\rm
LT}})^{2}(t-t_{m})}{2}\right)u_{m,n}^{{\rm LT}}$ is the explicit expression
(12) of the solution at time $t\in[t_{m},t_{m+1}]$ of the auxiliary stochastic
subsystem (10) used in the construction of the splitting integrator. Observe
that $u_{m,n}^{{\rm LT}}(t)=v_{m,n}^{M,N,1}(t)$ for all $t\in[t_{m},t_{m+1})$,
and, in particular, that $u_{m,n}^{{\rm LT}}(t_{m})=u_{m,n}^{{\rm LT}}$.
Moreover, by the construction of the splitting scheme, see (14), it holds that
$u_{m,n}^{{\rm LT}}(t_{m+1})=u_{m+1,n}^{{\rm LT}}$.

As a result, for any $n\in\\{1,\ldots,N-1\\}$, the mapping $u_{n}^{{\rm
LT}}\colon t\in[0,T]\mapsto u_{n}^{{\rm LT}}(t)$ defined such that
$u_{n}^{{\rm LT}}(t)=u_{m,n}^{{\rm LT}}(t)$ for $t\in[t_{m},t_{m+1}]$ is well-
defined. It is continuous on each interval $[t_{m},t_{m+1})$, and one has
$u_{n}^{{\rm LT}}(t_{m})=u_{m,n}^{{\rm LT}}$ for all $m\in\\{0,\ldots,M\\}$.

We claim that the following identity holds: for all $M\in{\mathbb{N}}$ and
$N\in{\mathbb{N}}$, for all $n\in\\{1,\ldots,N-1\\}$ and $t\in[0,T]$, one has

(20) $u_{n}^{{\rm LT}}(t)=\sum_{k=1}^{N-1}G_{nk}^{N}(\ell^{M}(t))u_{k,0}^{{\rm
LT}}+\sqrt{N}\int_{0}^{t}\sum_{k=1}^{N-1}G_{nk}^{N}(\ell^{M}(t)-\ell^{M}(s))u_{k}^{{\rm
LT}}(s)f(u_{k}^{{\rm LT}}(\ell^{M}(s))\,\text{d}W_{k}^{N}(s).$

The proof is based on a straightforward recursion argument.

Recall from Section 2.2 that one has the identities
$NG_{nk}^{N}(t)=G^{N}(t,x_{n},x_{k})$ and
$\sqrt{N}\text{d}W_{n}^{N}(t)=N\bigl{(}W(t,x_{n+1})-W(t,x_{n})\bigr{)}$. We
are now in position to provide the definition of the auxiliary process
$u^{{\rm LT}}$: for $t\in[0,T]$ and $x\in[0,1]$, define

(21) $\displaystyle u^{{\rm LT}}(t,x)$
$\displaystyle=\int_{0}^{1}G^{N}(t,x,y)u_{0}(\kappa^{N}(y))\,\text{d}y$
$\displaystyle\quad+\int_{0}^{t}\int_{0}^{1}G^{N}(\ell^{M}(t)-\ell^{M}(s),x,y)u^{{\rm
LT}}(s,\kappa^{N}(y))f(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))\,\text{d}W(s,y).$

In the identity (21) above, it is worth recalling that $x\in[0,1]\mapsto
G^{N}(t,x,y)$ is a piecewise linear mapping, whereas $y\in[0,1]\mapsto
G^{N}(t,x,y)$ is a piecewise constant mapping, with
$G^{N}(t,x_{n},x_{k})=NG_{nk}^{N}(t)$ for all $1\leq n,k\leq N-1$ and
$t\in[0,T]$.

Combining (20) and (21), one obtains the identity $u^{{\rm
LT}}(t,x_{n})=u_{n}^{{\rm LT}}(t)$ for all $t\in[0,T]$ and
$n\in\\{1,\ldots,N-1\\}$, and therefore one obtains the required property
$u^{{\rm LT}}(t_{m},x_{n})=u_{n}^{{\rm LT}}(t_{m})=u_{m,n}^{{\rm LT}}$. Note
that, for any $t\in[0,T]$, the mapping $x\in[0,1]\mapsto u^{{\rm LT}}(t,x)$ is
piecewise linear, more precisely it is linear on each subinterval
$[x_{n},x_{n+1}]$.

###  5.2. Auxiliary inequalities

In this subsection we state several inequalities used in the convergence
analysis of the splitting scheme.

  * •

For any continuous function $v\colon[0,1]\to{\mathbb{R}}$, one has (see for
instance [35, Eq. (3.5)])

(22) $\underset{N\in{\mathbb{N}}}{\sup}~{}\underset{t\geq
0}{\sup}~{}\underset{x\in[0,1]}{\sup}~{}\Big{|}\int_{0}^{1}G^{N}(t,x,y)v(\kappa^{N}(y))\,\text{d}y\Big{|}\leq\underset{x\in[0,1]}{\sup}~{}|v(x)|.$

  * •

For all $T\in(0,\infty)$, there exists $C_{T}\in(0,\infty)$ such that for all
$t\in(0,\infty)$ one has (see for instance [2, Lemma 2.3])

(23)
$\underset{N\in{\mathbb{N}}}{\sup}~{}\underset{x\in[0,1]}{\sup}~{}\int_{0}^{1}|G^{N}(t,x,y)|^{2}\,\text{d}y\leq\frac{C_{T}}{\sqrt{t}}.$

  * •

For all $T\in(0,\infty)$, there exists $C_{T}\in(0,\infty)$ such that for all
$t\in(0,T]$ and all $M\in{\mathbb{N}}$ one has

(24)
$\underset{N\in{\mathbb{N}}}{\sup}~{}\underset{x\in[0,1]}{\sup}~{}\int_{0}^{t}\int_{0}^{1}\big{|}G^{N}(t-s,x,y)-G^{N}(t-\ell^{M}(s),x,y)\big{|}^{2}\,\text{d}y\,\text{d}s\leq
C_{T}\sqrt{\tau}.$

Since we are not aware of a detailed proof of the inequality (24) in the
literature, we provide a proof in Appendix A. Note that the proof is similar
to the proof of [2, Lemma 2.3].

Let us also recall the following discrete Grönwall inequality, see for
instance [48, Lemma A.4]: assume that a sequence $\bigl{(}a_{m}\bigr{)}_{0\leq
m\leq M}$ of nonnegative numbers satisfies the inequality

$a_{m}\leq A+C\tau\sum_{k=0}^{m-1}\frac{a_{k}}{\sqrt{t_{m}-t_{k}}},$

where we recall that $t_{k}=k\tau=\frac{kT}{M}$, for some $A,C\in(0,\infty)$.
Then, there exists $C_{T}\in(0,\infty)$, depending only on $C$ and on $T$,
such that one has

(25) $\underset{0\leq m\leq M}{\sup}~{}a_{m}\leq C_{T}A.$

###  5.3. Moment bounds

The objective of this section is to prove Proposition 5. Recall that this
requires to impose the condition $\tau\leq\gamma h$ where we recall that
$\tau=T/M$, $h=1/N$ and where $\gamma\in(0,\infty)$ is an arbitrary parameter.

###### Proof of Proposition 5.

Using the definition (21) of the auxiliary process $u^{{\rm LT}}$, for all
$m\in\\{1,\ldots,M\\}$ and $n\in\\{1,\ldots,N-1\\}$, one has

$\displaystyle u_{m,n}^{{\rm LT}}$ $\displaystyle=u^{{\rm LT}}(t_{m},x_{n})$
$\displaystyle=\int_{0}^{1}G^{N}(t_{m},x_{n},y)u_{0}(\kappa^{N}(y))\,\text{d}y$
$\displaystyle\quad+\int_{0}^{t}\int_{0}^{1}G^{N}(t_{m}-\ell^{M}(s),x_{n},y)u^{{\rm
LT}}(s,\kappa^{N}(y))f(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))\,\text{d}W(s,y).$

Using Itô’s isometry formula, one obtains

$\displaystyle{\mathbb{E}}[|u_{m,n}^{{\rm LT}}|^{2}]$
$\displaystyle={\mathbb{E}}[\big{|}\int_{0}^{1}G^{N}(t_{m},x_{n},y)u_{0}(\kappa^{N}(y))\,\text{d}y\big{|}^{2}]$
$\displaystyle\quad+\int_{0}^{t}\int_{0}^{1}|G^{N}(t_{m}-\ell^{M}(s),x_{n},y)|^{2}{\mathbb{E}}[|u^{{\rm
LT}}(s,\kappa^{N}(y))|^{2}|f(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))|^{2}]\,\text{d}y\,\text{d}s.$

On the one hand, using the auxiliary inequality (22) and Assumption 1, one
obtains

${\mathbb{E}}[\big{|}\int_{0}^{1}G^{N}(t_{m},x_{n},y)u_{0}(\kappa^{N}(y))\,\text{d}y\big{|}^{2}]\leq\|u_{0}\|_{\infty}^{2}.$

On the other hand, recall that Assumption 3 implies that $f$ is bounded by
${\rm L_{g}}$. In addition, for all $k\in\\{0,\ldots,m-1\\}$ and all
$s\in[t_{k},t_{k+1})$, one has

${\mathbb{E}}[|u^{{\rm
LT}}(s,\kappa^{N}(y))|^{2}]={\mathbb{E}}[|v_{k,n}^{M,N,1}(s)|^{2}]$

where $n\in\\{1,\ldots,N-1\\}$ is such that $\kappa^{N}(y)=x_{n}$ and
$\bigl{(}v_{k,n}^{M,N,1}(s)\bigr{)}_{s\in[t_{k},t_{k+1}]}$ is the solution of
the auxiliary stochastic subsystem (10). Using the expression (12) for the
solution of (10) and the tower property of conditional expectation, one
obtains the upper bound

${\mathbb{E}}[|v_{k,n}^{M,N,1}(s)|^{2}]={\mathbb{E}}[e^{\frac{Nf(u_{k,n}^{{\rm
LT}})^{2}(s-t_{k})}{2}}|u_{k,n}^{{\rm LT}}|^{2}]\leq e^{\frac{N\tau{\rm
L_{g}}^{2}}{2}}{\mathbb{E}}[|u_{k,n}^{{\rm LT}}|^{2}]\leq e^{\frac{{\rm
L_{g}}^{2}\gamma}{2}}{\mathbb{E}}[|u_{k,n}^{{\rm LT}}|^{2}],$

using the boundedness of $f$ and the condition $N\tau\leq\gamma$.

Using the auxiliary inequality (23), gathering the upper bounds above yields
the following inequality: for all $m\in\\{1,\ldots,M\\}$ one has

$\underset{1\leq n\leq N-1}{\sup}~{}{\mathbb{E}}[|u_{m,n}^{{\rm
LT}}|^{2}]\leq\|u_{0}\|_{\infty}^{2}+C_{\gamma,T}\tau\sum_{k=0}^{m-1}\frac{1}{\sqrt{t_{m}-t_{k}}}\underset{1\leq
n\leq N-1}{\sup}~{}{\mathbb{E}}[|u_{k,n}^{{\rm LT}}|^{2}].$

Using the discrete Grönwall inequality (25) then gives

(26) $\underset{0\leq m\leq M}{\sup}~{}\underset{1\leq n\leq
N-1}{\sup}~{}{\mathbb{E}}[|u_{m,n}^{{\rm LT}}|^{2}]\leq
C_{\gamma,T}\|u_{0}\|_{\infty}^{2},$

where $C_{\gamma,T}\in(0,\infty)$ is independent of $M$, $N$ and
$\|u_{0}\|_{\infty}^{2}$. This shows moment bounds of the numerical solution
at the grid. It remains to extend this moment bound for $u^{{\rm
LT}}(t,x_{n})$ when $t$ is no longer assumed to be a grid point $t_{m}$.

For all $t\in[0,T)$ and $n\in\\{0,\ldots,N-1\\}$, let $m\in\\{0,\ldots,M-1\\}$
be such that $t_{m}=\ell^{M}(t)$, using the same arguments as above one has

${\mathbb{E}}[|u^{{\rm
LT}}(t,x_{n})|^{2}]={\mathbb{E}}[|v_{m,n}^{M,N,1}(t)|^{2}]\leq e^{\frac{{\rm
L_{g}}^{2}\gamma}{2}}{\mathbb{E}}[|u_{k,n}^{{\rm LT}}|^{2}]\leq
C_{\gamma,T}\|u_{0}\|_{\infty}^{2},$

where the inequality (26) is used in the last step. As a consequence, one has

(27) $\underset{t\in[0,T]}{\sup}~{}\underset{1\leq n\leq
N-1}{\sup}~{}{\mathbb{E}}[|u^{{\rm LT}}(t,x_{n})|^{2}]\leq
C_{\gamma,T}\|u_{0}\|_{\infty}^{2}.$

Finally, since $x\mapsto u^{{\rm LT}}(t,x)$ is linear on each subinterval
$[x_{n},x_{n+1}]$, one obtains

(28)
$\underset{t\in[0,T]}{\sup}~{}\underset{x\in[0,1]}{\sup}~{}{\mathbb{E}}[|u^{{\rm
LT}}(t,x)|^{2}]\leq\underset{t\in[0,T]}{\sup}~{}\underset{1\leq n\leq
N-1}{\sup}~{}{\mathbb{E}}[|u^{{\rm LT}}(t,x_{n})|^{2}]\leq
C_{\gamma,T}\|u_{0}\|_{\infty}^{2}.$

The proof of Proposition 5 is thus completed. ∎

A straightforward consequence of Proposition 5 is the following result.

######  Lemma 8.

Let Assumption 1 and Assumption 3 be satisfied. Let $\left(u^{{\rm
LT}}(t,x)\right)_{t\in[0,T],x\in[0,1]}$ be given by the mild formula (21).

For all $\gamma\in(0,\infty)$ and all $T\in(0,\infty)$, there exists
$C_{\gamma,T}\in(0,\infty)$ such that for all $\tau=T/M$ and $h=1/N$
satisfying the condition $\tau\leq\gamma h$, for all $m\in\\{0,\ldots,M-1\\}$
and all $t\in[t_{m},t_{m+1})$, one has

(29) $\underset{1\leq n\leq N-1}{\sup}~{}\bigl{(}{\mathbb{E}}[|u^{{\rm
LT}}(t,x_{n})-u^{{\rm LT}}(t_{m},x_{n})|^{2}]\bigr{)}^{\frac{1}{2}}\leq
C_{\gamma,T}(1+\|u_{0}\|_{\infty})\left(\frac{\tau}{h}\right)^{\frac{1}{2}}.$

###### Proof of Lemma 8.

Let $n\in\\{1,\ldots,N-1\\}$ and $m\in\\{0,\ldots,M-1\\}$, then for all
$t\in[t_{m},t_{m+1})$ one has

$\displaystyle u^{{\rm LT}}(t,x_{n})-u^{{\rm LT}}(t_{m},x_{n})$
$\displaystyle=v_{m,n}^{M,N,1}(t)-v_{m,n}^{M,N,1}(t_{m})$
$\displaystyle=\sqrt{N}\int_{t_{m}}^{t}v_{m,n}^{M,N,1}(s)f(u_{m,n}^{{\rm
LT}})\,\text{d}W_{n}^{N}(s)$ $\displaystyle=\sqrt{N}\int_{t_{m}}^{t}u^{{\rm
LT}}(s,x_{n})f(u_{m,n}^{{\rm LT}})\,\text{d}W_{n}^{N}(s),$

where we recall that the auxiliary process
$\bigl{(}v_{m,n}^{M,N,1}(t)\bigr{)}_{t_{m}\leq t\leq t_{m+1}}$ is defined by
the auxiliary subsystem (10) which gives the first step of the splitting
procedure, see Section 3.1.

Since the mapping $f$ is bounded, using Itô’s isometry formula, the condition
$\tau N\leq\gamma$ and the moment bounds (15) from Proposition 5, one obtains

${\mathbb{E}}[|u^{{\rm LT}}(t,x_{n})-u^{{\rm LT}}(t_{m},x_{n})|^{2}]\leq{\rm
L_{g}}^{2}N\tau{\mathbb{E}}[|u_{m,n}^{{\rm LT}}|^{2}]\leq{\rm
L_{g}}^{2}C_{\gamma,T}\bigl{(}1+\|u_{0}\|_{\infty}^{2}\bigr{)}\tau h^{-1}.$

The proof of Lemma 8 is thus completed. ∎

###  5.4. Convergence analysis

This section is devoted to the proof of the mean-square convergence of the
splitting scheme given in Theorem 6.

###### Proof of Theorem 6.

Recall that $u_{m,n}^{{\rm LT}}=u^{{\rm LT}}(t_{m},x_{n})$ for all
$n\in\\{1,\ldots,N-1\\}$ and $m\in\\{0,\ldots,M\\}$, where $\bigl{(}u^{{\rm
LT}}(t,x)\bigr{)}_{t\in[0,1],x\in[0,1]}$ is the process defined by (19).

For all $n\in\\{1,\ldots,N-1\\}$ and $m\in\\{1,\ldots,M\\}$, let us define

$E_{m,n}=u^{N}(t_{m},x_{n})-u_{m,n}^{{\rm LT}}\quad\text{and}\qquad
E_{m}=\underset{1\leq n\leq N-1}{\sup}~{}{\mathbb{E}}[|E_{m,n}|^{2}].$

Using the expression (6) for $u^{N}(t,x)$ and the expression (21) for $u^{{\rm
LT}}(t,x)$, one obtains the following decomposition of the error: for all
$n\in\\{1,\ldots,N-1\\}$ and $m\in\\{1,\ldots,M\\}$, one has

$\displaystyle E_{m,n}$ $\displaystyle=u^{N}(t_{m},x_{n})-u^{{\rm
LT}}(t_{m},x_{n})$
$\displaystyle=\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)g(u^{N}(s,\kappa^{N}(y)))\,\text{d}W(s,y)$
$\displaystyle-\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-\ell^{M}(s),x,y)u^{{\rm
LT}}(s,\kappa^{N}(y))f(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))\,\text{d}W(s,y)$
$\displaystyle=E_{m,n}^{(1)}+E_{m,n}^{(2)},$

where we set

$\displaystyle E_{m,n}^{(1)}$
$\displaystyle=\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)\bigl{[}g(u^{N}(s,\kappa^{N}(y)))-u^{{\rm
LT}}(s,\kappa^{N}(y))f(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))\bigr{]}\,\text{d}W(s,y),$ $\displaystyle
E_{m,n}^{(2)}$
$\displaystyle=\int_{0}^{t_{m}}\int_{0}^{1}\bigl{[}G^{N}(t_{m}-s,x,y)-G^{N}(t_{m}-\ell^{M}(s),x,y)\bigr{]}u^{{\rm
LT}}(s,\kappa^{N}(y))f(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))\,\text{d}W(s,y).$

Let us first deal with the error term $E_{m,n}^{(1)}$. Recall that
$g(u)=uf(u)$, therefore one has the decomposition
$E_{m,n}^{(1)}=E_{m,n}^{(1,1)}+E_{m,n}^{(1,2)}+E_{m,n}^{(1,3)}$, where

$\displaystyle E_{m,n}^{(1,1)}$
$\displaystyle=\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)\bigl{[}g(u^{N}(s,\kappa^{N}(y)))-g(u^{N}(\ell^{M}(s),\kappa^{N}(y)))\bigr{]}\,\text{d}W(s,y)$
$\displaystyle E_{m,n}^{(1,2)}$
$\displaystyle=\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)\bigl{[}g(u^{N}(\ell^{M}(s),\kappa^{N}(y)))-g(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y)))\bigr{]}\,\text{d}W(s,y)$ $\displaystyle
E_{m,n}^{(1,3)}$
$\displaystyle=\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)\bigl{[}u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))-u^{{\rm LT}}(s,\kappa^{N}(y))\bigr{]}f(u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y)))\,\text{d}W(s,y).$

Using Itô’s isometry formula, the global Lipschitz continuity assumption on
$g$, one obtains

$\displaystyle{\mathbb{E}}[|E_{m,n}^{(1,1)}|^{2}]$ $\displaystyle\leq{\rm
L_{g}}^{2}\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)^{2}{\mathbb{E}}[|u^{N}(s,\kappa^{N}(y))-u^{N}(\ell^{M}(s),\kappa^{N}(y))|^{2}]\,\text{d}y\,\text{d}s$
$\displaystyle\leq C_{T}(u_{0}){\rm
L_{g}}^{2}\sqrt{\tau}\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)^{2}\,\text{d}y\,\text{d}s$
$\displaystyle\leq C_{T}(u_{0})\sqrt{\tau},$

where we have used the temporal regularity estimate (9) for $u^{N}$ and the
auxiliary inequality (23).

Similarly, using Itô’s isometry formula, the global Lipschitz continuity
assumption on $g$, one obtains

$\displaystyle{\mathbb{E}}[|E_{m,n}^{(1,2)}|^{2}]$ $\displaystyle\leq{\rm
L_{g}}^{2}\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)^{2}{\mathbb{E}}[|u^{N}(\ell^{M}(s),\kappa^{N}(y))-u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))|^{2}]\,\text{d}y\,\text{d}s$
$\displaystyle\leq
C\sum_{k=0}^{m-1}E_{k}\int_{t_{k}}^{t_{k+1}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)^{2}\,\text{d}y\,\text{d}s.$

Using the inequality (23), for all $k\in\\{0,\ldots,m-1\\}$, one has

$\displaystyle\int_{t_{k}}^{t_{k+1}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)^{2}\,\text{d}y\,\text{d}s$
$\displaystyle\leq\int_{t_{k}}^{t_{k+1}}\frac{C_{T}}{\sqrt{t_{m}-s}}\,\text{d}s$
$\displaystyle=2C_{T}\bigl{(}\sqrt{t_{m}-t_{k}}-\sqrt{t_{m}-t_{k+1}}\bigr{)}$
$\displaystyle=2C_{T}\sqrt{t_{m}-t_{k}}\Bigl{(}1-\sqrt{1-\frac{\tau}{t_{m}-t_{k}}}\Bigr{)}$
$\displaystyle\leq\frac{2C_{T}\tau}{\sqrt{t_{m}-t_{k}}},$

where we have used the inequality $1-\sqrt{1-z}\leq z$ for all $z\in[0,1]$ in
the last step. Therefore one has

${\mathbb{E}}[|E_{m,n}^{(1,2)}|^{2}]\leq
C_{T}\tau\sum_{k=0}^{m-1}\frac{E_{k}}{\sqrt{t_{m}-t_{k}}}.$

Finally, for the third term, using Itô’s isometry formula and the boundedness
of $f$, one obtains

$\displaystyle{\mathbb{E}}[|E_{m,n}^{(1,3)}|^{2}]$ $\displaystyle\leq
L_{g}^{2}\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)^{2}{\mathbb{E}}[|u^{{\rm
LT}}(\ell^{M}(s),\kappa^{N}(y))-u^{{\rm
LT}}(s,\kappa^{N}(y))|^{2}]\,\text{d}y\,\text{d}s$ $\displaystyle\leq
C_{\gamma,T}(u_{0})\tau
h^{-1}\int_{0}^{t_{m}}\int_{0}^{1}G^{N}(t_{m}-s,x,y)^{2}\,\text{d}y\,\text{d}s$
$\displaystyle\leq C_{\gamma,T}(u_{0})\tau h^{-1}$

using the temporal regularity estimate (29) from Lemma 8 for $u^{{\rm LT}}$
and the auxiliary inequality (23).

Let us now deal with the error term $E_{m,n}^{(2)}$. Using Itô’s formula, the
boundedness of $f$ and the moment bounds (15) from Proposition 5, one obtains

$\displaystyle{\mathbb{E}}[|E_{m,n}^{(2)}|^{2}]$ $\displaystyle\leq{\rm
L_{g}}^{2}\int_{0}^{t_{m}}\int_{0}^{1}\big{|}G^{N}(t-\ell^{M}(s),x,y)-G^{N}(t_{m}-\ell^{M}(s),x,y)\big{|}^{2}{\mathbb{E}}[|u^{{\rm
LT}}(s,\kappa^{N}(y))|^{2}]\,\text{d}y\,\text{d}s$ $\displaystyle\leq
C_{\gamma,T}(u_{0})\int_{0}^{t_{m}}\int_{0}^{1}\big{|}G^{N}(t-\ell^{M}(s),x,y)-G^{N}(t_{m}-\ell^{M}(s),x,y)\big{|}^{2}\,\text{d}y\,\text{d}s$
$\displaystyle\leq C_{\gamma,T}(u_{0})\sqrt{\tau},$

owing to the auxiliary inequality (24) in the last step.

Gathering the estimates, for all $m\in\\{1,\ldots,M\\}$, one has

$E_{m}\leq C_{\gamma,T}(u_{0})\left(\sqrt{\tau}+\tau
h^{-1}\right)+C_{T}\tau\sum_{k=0}^{m-1}\frac{E_{k}}{\sqrt{t_{m}-t_{k}}}.$

Applying the discrete Grönwall inequality (25) (see Section 5.2) then yields

$\underset{0\leq m\leq M}{\sup}~{}E_{m}\leq
C_{\gamma,T}(u_{0})\left(\sqrt{\tau}+\frac{\tau}{h}\right).$

This gives the error estimate (16). When the condition $\tau\leq\gamma h^{2}$
is satisfied, one has $\tau h^{-1}\leq\sqrt{\gamma}\tau^{\frac{1}{2}}$ and one
has the error estimate (17). This concludes the proof of Theorem 6. ∎

Let us also provide the proof of Corollary 7.

###### Proof of Corollary 7.

It suffices to combine the error estimate (8) from Proposition 3 for the
spatial discretization error, and the error estimate (16) from Theorem 6 for
the temporal discretization error. One then obtains the error estimate for the
splitting scheme

$\displaystyle\bigl{(}{\mathbb{E}}[|u_{m,n}^{{\rm
LT}}-u(t_{m},x_{n})|^{2}]\bigr{)}^{\frac{1}{2}}$
$\displaystyle\leq\bigl{(}{\mathbb{E}}[|u_{m,n}^{{\rm
LT}}-u^{N}(t_{m},x_{n})|^{2}]\bigr{)}^{\frac{1}{2}}+\bigl{(}{\mathbb{E}}[|u^{N}(t_{m},x_{n})-u(t_{m},x_{n})|^{2}]\bigr{)}^{\frac{1}{2}}$
$\displaystyle\leq
C_{\gamma,T}(u_{0})\tau^{\frac{1}{4}}+C_{T}(u_{0})h^{\frac{1}{2}}$
$\displaystyle\leq
C_{\gamma,T}(u_{0})\gamma^{\frac{1}{4}}h^{\frac{1}{2}}+C_{T}(u_{0})h^{\frac{1}{2}},$

under the condition $\tau\leq\gamma h^{2}$. This gives the error estimate (18)
and concludes the proof of Corollary 7. ∎

###  5.5. Proof of Proposition 2

We conclude this section with the proof of the positivity property of the
exact solution to the stochastic heat equation (3) on a bounded domain.

###### Proof of Proposition 2.

Owing to Corollary 7 and to the temporal regularity estimate (5) satisfied by
the solution $u$ of the SPDE in equation (3), one obtains the following result
(recall that $\tau=T/M$ and $h=1/N$): there exists
$C_{\gamma,T}(u_{0})\in(0,\infty)$ such that for all $N\in{\mathbb{N}}$ and
$M\in{\mathbb{N}}$, such that $M\geq\frac{TN^{2}}{\gamma}$, for all
$t\in[0,T]$ and $x\in[0,1]$, one has

(30) $\bigl{(}{\mathbb{E}}[|u(t,x)-u^{{\rm
LT}}(\ell^{M}(t),\kappa^{N}(x))|^{2}]\bigr{)}^{\frac{1}{2}}\leq
C_{\gamma,T}(u_{0})N^{-\frac{1}{2}}.$

Let $t\in[0,T]$ and $x\in[0,1]$ be fixed, then there exists a sequence
$\bigl{(}N_{k})_{k\in{\mathbb{N}}}$ such that $N_{k}\to\infty$ and $u^{{\rm
LT}}(\ell^{M_{k}}(t),\kappa^{N_{k}}(x))$ converges to $u(t,x)$ almost surely.
Since $u^{{\rm LT}}(\ell^{M_{k}}(t),\kappa^{N_{k}}(x))\geq 0$ almost surely
owing to Proposition 4, one obtains $u(t,x)\geq 0$ almost surely. ∎

##  6\. Generalization to systems

In this section, we briefly describe how to generalize the construction of the
splitting scheme (14) and the analysis above to stochastic systems of the type

(31) $\left\\{\begin{aligned}
&\text{d}u_{1}(t,x)=\partial_{xx}^{2}u_{1}(t,x)\,\text{d}t+g_{1}(u_{1}(t,x),u_{2}(t,x))\,\text{d}W_{1}(t,x),\\\
&\text{d}u_{2}(t,x)=\partial_{xx}^{2}u_{2}(t,x)\,\text{d}t+g_{2}(u_{1}(t,x),u_{2}(t,x))\,\text{d}W_{2}(t,x),\\\
&u_{1}(t,0)=u_{1}(t,1)=0,\quad u_{2}(t,0)=u_{2}(t,1)=0,\\\
&u_{1}(0,x)=u_{1,0}(x),\quad u_{2}(0,x)=u_{2,0}(x),\end{aligned}\right.$

for $(t,x)\in[0,T]\times[0,1]$, where
$g_{1},g_{2}\colon{\mathbb{R}}^{2}\to{\mathbb{R}}$ are globally Lipschitz
continuous mappings, with initial values $u_{1,0},u_{2,0}$ satisfying
Assumptions 1 and 2. The two evolution equations are driven by space-time
white noise. The Wiener sheets $W_{1}$ and $W_{2}$ can either be equal or
independent. For ease of presentation we only deal with systems of two
equations, while considering systems of arbitrary size would also be possible.

In this setting, to obtain solutions which only have nonnegative values, it is
necessary to replace Assumption 3 by the following.

######  Assumption 4.

The mappings $g_{1},g_{2}\colon{\mathbb{R}}^{2}\to{\mathbb{R}}$ are of class
$\mathcal{C}^{1}$ and globally Lipschitz continuous. In addition, they satisfy
$g_{1}(0,v_{2})=0$ and $g_{2}(v_{1},0)=0$ for all
$(v_{1},v_{2})\in{\mathbb{R}}^{2}$.

One then has the following generalization of Proposition 2.

######  Proposition 9.

Consider the SPDE system (31). Let Assumption 4 be satisfied and assume that
the initial values $u_{1,0},u_{2,0}$ satisfy Assumptions 1 and 2. Then, for
all $t\in(0,\infty)$ and all $x\in[0,1]$, almost surely, one has

$u_{1}(t,x)\geq 0~{},\quad u_{2}(t,x)\geq 0.$

As in Sections 2 and 3, the mesh size and the time-step sizes are denoted by
$h=1/N$ and $\tau=T/M$ respectively, and the space and time grid points are
denoted by $x_{n}=nh$ and $t_{m}=m\tau$, with $0\leq n\leq N$ and $0\leq m\leq
M$. In addition, introduce the mappings
$f_{1},f_{2}\colon{\mathbb{R}}^{2}\to{\mathbb{R}}$ defined by

$\displaystyle f_{1}(v_{1},v_{2})$
$\displaystyle=\frac{g_{1}(v_{1},v_{2})}{v_{1}}=\int_{0}^{1}\partial_{v_{1}}g_{1}(rv_{1},v_{2})\,\text{d}r~{},\quad
f_{2}(v_{1},v_{2})$
$\displaystyle=\frac{g_{2}(v_{1},v_{2})}{v_{2}}=\int_{0}^{1}\partial_{v_{2}}g_{2}(v_{1},rv_{2})\,\text{d}r.$

Owing to Assumption 4, the mappings $f_{1}$ and $f_{2}$ are bounded and
continuous mappings. Finally, for all $t\geq 0$ and $n\in\\{1,\ldots,N-1\\}$
define

$W_{1,n}^{N}(t)=\sqrt{N}\bigl{(}W_{1}(t,x_{n+1})-W_{1}(t,x_{n})\bigr{)}~{},\quad
W_{2,n}^{N}(t)=\sqrt{N}\bigl{(}W_{2}(t,x_{n+1})-W_{2}(t,x_{n})\bigr{)}$

and define the noise increments

$\Delta_{m,n}W_{1}=W_{1,n}^{N}(t_{m+1})-W_{1,n}^{N}(t_{m})~{},\quad\Delta_{m,n}W_{2}=W_{2,n}^{N}(t_{m+1})-W_{2,n}^{N}(t_{m})$

for all $n\in\\{1,\ldots,N-1\\}$ and $m\in\\{0,\ldots,M-1\\}$.

Using the finite difference method and the same notation as in Section 2.2,
one obtains the spatial semi-discretization scheme for the SPDE system (31)
with mesh size $h$ as follows:

(32) $\left\\{\begin{aligned}
\text{d}u_{1}^{N}(t)&=N^{2}D^{N}u_{1}^{N}(t)\,\text{d}t+\sqrt{N}g_{1}(u_{1}^{N}(t),u_{2}^{N}(t))\,\text{d}W_{1}^{N}(t)\\\
\text{d}u_{2}^{N}(t)&=N^{2}D^{N}u_{2}^{N}(t)\,\text{d}t+\sqrt{N}g_{2}(u_{1}^{N}(t),u_{2}^{N}(t))\,\text{d}W_{2}^{N}(t).\end{aligned}\right.$

We are now in position to state the definition of the fully-discrete scheme
based on a Lie–Trotter splitting strategy and inspired by (14) for the
approximation of solutions of (31): for all $m\in\\{0,\ldots,M-1\\}$, set

(33) $\left\\{\begin{aligned} u_{1,m+1}^{{\rm LT}}&=e^{\tau
N^{2}D^{N}}\left(\exp\Bigl{(}\sqrt{N}f_{1}(u_{1,m,n}^{{\rm
LT}},u_{2,m,n}^{{\rm LT}})\Delta_{m,n}W_{1}-\frac{Nf_{1}(u_{1,m,n}^{{\rm
LT}},u_{2,m,n}^{{\rm LT}})^{2}\tau}{2}\Bigr{)}u_{1,m,n}^{{\rm
LT}}\right)_{1\leq n\leq N-1}\\\ u_{2,m+1}^{{\rm LT}}&=e^{\tau
N^{2}D^{N}}\left(\exp\Bigl{(}\sqrt{N}f_{2}(u_{1,m,n}^{{\rm
LT}},u_{2,m,n}^{{\rm LT}})\Delta_{m,n}W_{2}-\frac{Nf_{2}(u_{1,m,n}^{{\rm
LT}},u_{2,m,n}^{{\rm LT}})^{2}\tau}{2}\Bigr{)}u_{2,m,n}^{{\rm
LT}}\right)_{1\leq n\leq N-1},\end{aligned}\right.$

with initial values $u_{1,0}^{{\rm LT}}=\bigl{(}u_{1,0}(x_{n})\bigr{)}_{1\leq
n\leq N-1}$ and $u_{1,0}^{{\rm LT}}=\bigl{(}u_{2,0}(x_{n})\bigr{)}_{1\leq
n\leq N-1}$.

The scheme (33) is positivity-preserving in the following sense.

######  Proposition 10.

Let $M\in{\mathbb{N}}$ and $N\in{\mathbb{N}}$ be arbitrary integers and let
$T\in(0,\infty)$. Let Assumption 4 be satisfied, and assume that the initial
values $u_{1,0},u_{2,0}$ satisfy Assumptions 1 and 2. Let the sequence
$u_{1,0}^{{\rm LT}},\ldots,u_{1,M}^{{\rm LT}}$ and $u_{2,0}^{{\rm
LT}},\ldots,u_{2,M}^{{\rm LT}}$ be given by the splitting scheme (33), with
$h=1/N$ and $\tau=T/M$, with initial values $u_{1,0,n}^{{\rm
LT}}=u_{1,0}(x_{n})\geq 0$ and $u_{1,0,n}^{{\rm LT}}=u_{2,0}(x_{n})\geq 0$ for
all $n\in\\{1,\ldots,N\\}$. Then, almost surely, one has

$u_{1,m,n}^{{\rm LT}}\geq 0~{},\quad u_{2,m,n}^{{\rm LT}}\geq 0,$

for all $m\in\\{1,\ldots,M\\}$ and $n\in\\{1,\ldots,N-1\\}$.

The proof of Proposition 10 is a straightforward modification of the proof of
Proposition 4. Moreover, one has the following variant of Proposition 5.

######  Proposition 11.

Let Assumption 4 be satisfied and assume that the initial values
$u_{1,0},u_{2,0}$ satisfy Assumptions 1 and 2. Let the sequences
$u_{1,0}^{{\rm LT}},\ldots,u_{1,M}^{{\rm LT}}$ and $u_{2,0}^{{\rm
LT}},\ldots,u_{2,M}^{{\rm LT}}$ be given by the Lie–Trotter splitting scheme
(33).

For all $\gamma\in(0,\infty)$ and all $T\in(0,\infty)$, there exists
$C_{\gamma,T}\in(0,\infty)$ such that for all $\tau=T/M$ and $h=1/N$
satisfying the condition $\tau\leq\gamma h$, one has

(34) $\underset{0\leq m\leq M}{\sup}~{}\underset{1\leq n\leq
N-1}{\sup}~{}{\mathbb{E}}\left[|u_{1,m,n}^{{\rm
LT}}|^{2}\right]+\underset{0\leq m\leq M}{\sup}~{}\underset{1\leq n\leq
N-1}{\sup}~{}{\mathbb{E}}\left[|u_{2,m,n}^{{\rm LT}}|^{2}\right]\leq
C_{\gamma,T}\bigl{(}1+\|u_{1,0}\|_{\infty}^{2}+\|u_{2,0}\|_{\infty}^{2}\bigr{)}.$

Finally, one has the following generalization of Theorem 6.

######  Theorem 12.

Let Assumption 4 be satisfied and assume that the initial values
$u_{1,0},u_{2,0}$ satisfy Assumptions 1 and 2. Let the sequences
$u_{1,0}^{{\rm LT}},\ldots,u_{1,M}^{{\rm LT}}$ and $u_{2,0}^{{\rm
LT}},\ldots,u_{2,M}^{{\rm LT}}$ be given by the Lie–Trotter splitting scheme
(33), and let $\bigl{(}u_{1}^{N}(t)\bigr{)}_{t\geq 0,0\leq n\leq N}$ and
$\bigl{(}u_{2}^{N}(t)\bigr{)}_{t\geq 0,0\leq n\leq N}$ be given by the spatial
semi-discretization scheme (32).

For all $\gamma\in(0,\infty)$ and $T\in(0,\infty)$, there exists
$C_{\gamma,T}(u_{1,0},u_{2,0})\in(0,\infty)$ such that for all $\tau=T/M$ and
$h=1/N$ satisfying the condition $\tau\leq\gamma h$, one has

(35) $\displaystyle\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{1,m,n}^{{\rm
LT}}-u_{1,n}^{N}(t_{m})|^{2}]\right)^{\frac{1}{2}}$ $\displaystyle\leq
C_{\gamma,T}(u_{1,0},u_{2,0})\left(\tau^{\frac{1}{4}}+\left(\frac{\tau}{h}\right)^{\frac{1}{2}}\right)$
$\displaystyle\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{2,m,n}^{{\rm
LT}}-u_{2,n}^{N}(t_{m})|^{2}]\right)^{\frac{1}{2}}$ $\displaystyle\leq
C_{\gamma,T}(u_{1,0},u_{2,0})\left(\tau^{\frac{1}{4}}+\left(\frac{\tau}{h}\right)^{\frac{1}{2}}\right).$

In addition, for all $\tau=T/M$ and $h=1/N$ satisfying the condition
$\tau\leq\gamma h^{2}$, one has

(36) $\displaystyle\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{1,m,n}^{{\rm
LT}}-u_{1,n}^{N}(t_{m})|^{2}]\right)^{\frac{1}{2}}$ $\displaystyle\leq
C_{\gamma,T}(u_{1,0},u_{2,0})\tau^{\frac{1}{4}}$ $\displaystyle\underset{0\leq
m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{2,m,n}^{{\rm
LT}}-u_{2,n}^{N}(t_{m})|^{2}]\right)^{\frac{1}{2}}$ $\displaystyle\leq
C_{\gamma,T}(u_{1,0},u_{2,0})\tau^{\frac{1}{4}}.$

The proofs of Proposition 11 and of Theorem 12 are omitted since they follow
from the same arguments as those of Proposition 5 and of Theorem 6. Finally,
one obtains the following variant of Corollary 7

######  Corollary 13.

Consider the setting and assumptions of Theorem 12. For all
$\gamma\in(0,\infty)$ and $T\in(0,\infty)$, there exists
$C_{\gamma,T}(u_{1,0},u_{2,0})\in(0,\infty)$ such that for all $\tau=T/M$ and
$h=1/N$ satisfying the condition $\tau\leq\gamma h^{2}$, one has

(37) $\displaystyle\underset{0\leq m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{1,m,n}^{{\rm
LT}}-u_{1}(t_{m},x_{n})|^{2}]\right)^{\frac{1}{2}}$ $\displaystyle\leq
C_{\gamma,T}(u_{1,0},u_{2,0})h^{\frac{1}{2}}$ $\displaystyle\underset{0\leq
m\leq M}{\sup}~{}\underset{0\leq n\leq
N}{\sup}~{}\left({\mathbb{E}}[|u_{2,m,n}^{{\rm
LT}}-u_{2}(t_{m},x_{n})|^{2}]\right)^{\frac{1}{2}}$ $\displaystyle\leq
C_{\gamma,T}(u_{1,0},u_{2,0})h^{\frac{1}{2}}.$

To conclude this presentation of the positivity-preserving Lie–Trotter
splitting scheme (33) for the approximation of solutions of the SPDE system
(31), we report some numerical experiments.

The first numerical experiment illustrates the positivity-preserving property
of the Lie–Trotter splitting scheme (LT) when applied to the system of SPDEs
(31) driven by two independent noise. The initial values are taken to be
$u_{1,0}=u_{2,0}=\sin(\pi x)$, the final time is $T=5$ and the multiplicative
terms are $g_{1}(v_{1},v_{2})=7\sin(v_{1})\cos(v_{2})$ and
$g_{2}(v_{1},v_{2})=7\cos(v_{1})\sin(v_{2})$. The discretization parameters
are $\tau=2^{-2}$ and $h=2^{-8}$. The proportion of samples containing only
positive values out of $500$ simulated samples for all considered time
integrators are presented in Table 3.

LT (first,second) | SEXP (first,second) | SEM (first,second) | EM (first,second)
---|---|---|---
$500/500,500/500$ | $500/500,499/500$ | $498/500,496/500$ | $0/500,0/500$
Table 3. Proportion of samples containing only positive values out of $500$
simulated sample paths for the Lie–Trotter splitting scheme (LT), the
stochastic exponential Euler integrator (SEXP), the semi-implicit
Euler–Maruyama scheme (SEM), and the Euler–Maruyama scheme (EM). First and
second component. The multiplicative terms are
$g_{1}(v_{1},v_{2})=7\sin(v_{1})\cos(v_{2})$ and
$g_{2}(v_{1},v_{2})=7\cos(v_{1})\sin(v_{2})$. The discretization parameters
are $\tau=2^{-2}$ and $h=2^{-8}$.

The second numerical experiment illustrates the mean-square convergence of the
Lie–Trotter splitting scheme when applied to systems of nonlinear SHEs. Figure
4 presents, in a loglog plot, the mean-square errors measured at the space-
time grid for the time interval $[0,0.5]$. The discretization parameters are
$h=2^{-8}$ and $\tau=2^{-4},2^{-5},\ldots,2^{-16}$ (the last one being used
for the reference solution). We have used $200$ samples to approximate the
expected values. The expected mean-square orders of convergence is observed in
this figure.

Figure 4. Mean-square errors of the Lie–Trotter splitting scheme (first
component denoted by LT1, second by LT2) when applied to the system of
stochastic heat equations with multiplicative terms
$g_{1}(v_{1},v_{2})=\sin(v_{1})\cos(v_{2})$ and
$g_{2}(v_{1},v_{2})=\cos(v_{1})\sin(v_{2})$. Mesh size $h=2^{-8}$ and average
over $200$ samples.

##  Appendix A Proof of auxiliary inequalities

###### Proof of the auxiliary inequality (24).

Let us recall some notation. For all $N\in{\mathbb{N}}$, all $t\geq 0$ and
$x,y\in[0,1]$, one has

$G^{N}(t,x,y)=\sum_{j=1}^{N-1}e^{-\lambda_{j}^{N}t}\varphi_{j}^{N}(x)\varphi_{j}(\kappa^{N}(y)),$

where $\lambda_{j}^{N}=4N^{2}\sin\bigl{(}\frac{j\pi}{2N}\bigr{)}^{2}$,
$\varphi_{j}(\cdot)=\sqrt{2}\sin(j\pi\cdot)$ and $\varphi_{j}^{N}$ is the
linear interpolation of $\varphi_{j}$ at the space grid points $x_{n}=nh$ for
$n=1,\ldots,N-1$.

Using the orthogonality property

$\int_{0}^{1}\varphi_{j}(\kappa^{N}(y))\varphi_{k}(\kappa^{N}(y))\,\text{d}y=\delta_{jk},$

one obtains

$\displaystyle\int_{0}^{t}\int_{0}^{1}\big{|}G^{N}(t-s,x,y)$
$\displaystyle-G^{N}(t-\ell^{M}(s),x,y)\big{|}^{2}\,\text{d}y\,\text{d}s$
$\displaystyle=\int_{0}^{t}\int_{0}^{1}\big{|}\sum_{j=1}^{N-1}\bigl{(}e^{-\lambda_{j}^{N}(t-s)}-e^{-\lambda_{j}^{N}(t-\ell^{M}(s))}\bigr{)}\varphi_{j}^{N}(x)\varphi_{j}(\kappa^{N}(y))\big{|}^{2}\,\text{d}y\,\text{d}s$
$\displaystyle=\int_{0}^{t}\sum_{j=1}^{N-1}\bigl{(}e^{-\lambda_{j}^{N}(t-s)}-e^{-\lambda_{j}^{N}(t-\ell^{M}(s))}\bigr{)}^{2}\varphi_{j}^{N}(x)^{2}\,\text{d}s$
$\displaystyle\leq
2\int_{0}^{t}\sum_{j=1}^{N-1}\bigl{(}e^{-\lambda_{j}^{N}(t-s)}-e^{-\lambda_{j}^{N}(t-\ell^{M}(s))}\bigr{)}^{2}\,\text{d}s$
$\displaystyle\leq
2\sum_{j=1}^{N-1}\int_{0}^{t}e^{-2\lambda_{j}^{N}(t-s)}\bigl{(}1-e^{-\lambda_{j}^{N}(s-\ell^{M}(s))}\bigr{)}^{2}\,\text{d}s$
$\displaystyle\leq
C\sum_{j=1}^{N-1}\frac{\max(1,\lambda_{j}^{N}\tau)^{2}}{\lambda_{j}^{N}}.$

One checks that there exists $c\in(1,\infty)$ such that for all $N\geq 1$ and
$j\in\\{1,\ldots,N-1\\}$ one has

$c^{-1}\leq\frac{\lambda_{j}^{N}}{j^{2}}\leq c.$

Let $L\in{\mathbb{N}}$ be an arbitrary positive integer. Owing to the
inequalities above, one obtains

$\displaystyle\sum_{j=1}^{N-1}\frac{\max(1,\lambda_{j}^{N}\tau)^{2}}{\lambda_{j}^{N}}$
$\displaystyle\leq C\sum_{j=1}^{\infty}\frac{\max(1,j^{2}\tau)^{2}}{j^{2}}$
$\displaystyle\leq C\sum_{j=1}^{L}j^{2}\tau^{2}+C\sum_{j=L+1}^{\infty}j^{-2}$
$\displaystyle\leq C\tau^{2}L^{3}+CL^{-1},$

using standard comparison of series and integrals arguments. Choosing
$L=\lfloor\tau^{-\frac{1}{2}}\rfloor\geq 1$ (where $\lfloor\cdot\rfloor$
denotes the integer part), and recalling that $\tau\in(0,1)$, one obtains

$\int_{0}^{t}\int_{0}^{1}\big{|}G^{N}(t-s,x,y)-G^{N}(t-\ell^{M}(s),x,y)\big{|}^{2}\,\text{d}y\,\text{d}s\leq
C\tau^{\frac{1}{2}}.$

The value of $C$ is independent of $N\in{\mathbb{N}}$, $t\in(0,T]$ and
$x\in[0,1]$. The proof of the auxiliary inequality (24) is thus completed. ∎

## Acknowledgements

The work of CEB is partially supported by the project SIMALIN
(ANR-19-CE40-0016) operated by the French National Research Agency. The work
of DC and JU is partially supported by the Swedish Research Council (VR)
(projects nr. $2018-04443$). The computations were performed on resources
provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N,
Umeå University and at UPPMAX, Uppsala University.

## References

  * [1] K. Abiko and T. Ishiwata.  Positivity-preserving numerical schemes for stochastic differential equations.  Jpn. J. Ind. Appl. Math., 39(3):1095–1108, 2022.
  * [2] R. Anton, D. Cohen, and L. Quer-Sardanyons.  A fully discrete approximation of the one-dimensional stochastic heat equation.  IMA J. Numer. Anal., 40(1):247–284, 2020.
  * [3] V. Barbu and M. Röckner.  A splitting algorithm for stochastic partial differential equations driven by linear multiplicative noise.  Stoch. Partial Differ. Equ. Anal. Comput., 5(4):457–471, 2017.
  * [4] A. Barth and A. Lang.  Simulation of stochastic partial differential equations using finite element methods.  Stochastics, 84(2-3):217–231, 2012.
  * [5] A. Barth and A. Lang.  $L^{p}$ and almost sure convergence of a Milstein scheme for stochastic partial differential equations.  Stochastic Process. Appl., 123(5):1563–1587, 2013.
  * [6] C. Bauzet, F. Nabet, K. Schmitz, and A. Zimmermann.  Convergence of a finite-volume scheme for a heat equation with a multiplicative lipschitz noise, 2022.
  * [7] C. Bayer and H. Oberhauser.  Splitting methods for SPDEs: from robustness to financial engineering, optimal control, and nonlinear filtering.  In Splitting methods in communication, imaging, science, and engineering, Sci. Comput., pages 499–539. Springer, Cham, 2016.
  * [8] F. E. Benth.  On the positivity of the stochastic heat equation.  Potential Anal., 6(2):127–148, 1997.
  * [9] A. Berg, D. Cohen, and G. Dujardin.  Lie-Trotter splitting for the nonlinear stochastic Manakov system.  J. Sci. Comput., 88(1):Paper No. 6, 31, 2021.
  * [10] S. Blanes and F. Casas.  A concise introduction to geometric numerical integration.  Monographs and Research Notes in Mathematics. CRC Press, Boca Raton, FL, 2016.
  * [11] M. Bossy, J.-F. Jabir, and K. Martínez.  On the weak convergence rate of an exponential Euler scheme for SDEs governed by coefficients with superlinear growth.  Bernoulli, 27(1):312–347, 2021.
  * [12] C.-E. Bréhier and D. Cohen.  Strong rates of convergence of a splitting scheme for Schrödinger equations with nonlocal interaction cubic nonlinearity and white noise dispersion.  SIAM/ASA J. Uncertain. Quantif., 10(1):453–480, 2022.
  * [13] C.-E. Bréhier and D. Cohen.  Analysis of a splitting scheme for a class of nonlinear stochastic Schrödinger equations, 2023.
  * [14] C.-E. Bréhier, D. Cohen, and G. Giordano.  Splitting schemes for FitzHugh-Nagumo stochastic partial differential equations, 2022.
  * [15] C.-E. Bréhier, J. Cui, and J. Hong.  Strong convergence rates of semidiscrete splitting approximations for the stochastic Allen-Cahn equation.  IMA J. Numer. Anal., 39(4):2096–2134, 2019.
  * [16] C.-E. Bréhier and L. Goudenège.  Weak convergence rates of splitting schemes for the stochastic Allen-Cahn equation.  BIT, 60(3):543–582, 2020.
  * [17] O. Butkovsky, K. Dareiotis, and M. Gerencsér.  Optimal rate of convergence for approximations of spdes with non-regular drift, 2021.
  * [18] R. A. Carmona and S. A. Molchanov.  Parabolic Anderson problem and intermittency.  Mem. Amer. Math. Soc., 108(518):viii+125, 1994.
  * [19] L. Chen and J. Huang.  Comparison principle for stochastic heat equation on $\mathbb{R}^{d}$.  Ann. Probab., 47(2):989–1035, 2019.
  * [20] L. Chen and K. Kim.  On comparison principle and strict positivity of solutions to the nonlinear stochastic fractional heat equations.  Ann. Inst. Henri Poincaré Probab. Stat., 53(1):358–388, 2017\.
  * [21] S. Cox and J. van Neerven.  Convergence rates of the splitting scheme for parabolic linear stochastic Cauchy problems.  SIAM J. Numer. Anal., 48(2):428–451, 2010.
  * [22] S. Cox and J. van Neerven.  Pathwise Hölder convergence of the implicit-linear Euler scheme for semi-linear SPDEs with multiplicative noise.  Numer. Math., 125(2):259–345, 2013.
  * [23] J. Cresson, M. Efendiev, and S. Sonner.  On the positivity of solutions of systems of stochastic PDEs.  ZAMM Z. Angew. Math. Mech., 93(6-7):414–422, 2013.
  * [24] J. Cui, J. Hong, Z. Liu, and W. Zhou.  Strong convergence rate of splitting schemes for stochastic nonlinear Schrödinger equations.  J. Differential Equations, 266(9):5625–5663, 2019.
  * [25] G. Da Prato and J. Zabczyk.  Stochastic equations in infinite dimensions, volume 152 of Encyclopedia of Mathematics and its Applications.  Cambridge University Press, Cambridge, second edition, 2014.
  * [26] R. Dalang, D. Khoshnevisan, C. Mueller, D. Nualart, and Y. Xiao.  A minicourse on stochastic partial differential equations, volume 1962 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin, 2009.  Held at the University of Utah, Salt Lake City, UT, May 8–19, 2006, Edited by Khoshnevisan and Firas Rassoul-Agha.
  * [27] A. M. Davie and J. G. Gaines.  Convergence of numerical schemes for the solution of parabolic stochastic partial differential equations.  Math. Comp., 70(233):121–134, 2001.
  * [28] A. Deya.  Numerical schemes for rough parabolic equations.  Appl. Math. Optim., 65(2):253–292, 2012.
  * [29] A. Deya and R. Marty.  A full discretization of the rough fractional linear heat equation.  Electron. J. Probab., 27:Paper No. 122, 41, 2022.
  * [30] R. Duboscq and R. Marty.  Analysis of a splitting scheme for a class of random nonlinear partial differential equations.  ESAIM Probab. Stat., 20:572–589, 2016.
  * [31] U. Erdoğan and G. J. Lord.  A new class of exponential integrators for sdes with multiplicative noise.  IMA J. Numer. Anal., 39(2):820–846, 2019.
  * [32] L. Farina and S. Rinaldi.  Positive Linear Systems: Theory and Applications.  John Wiley & Sons, Incorporated, New York, 2000.
  * [33] M. Gerencsér and I. Gyöngy.  Finite difference schemes for stochastic partial differential equations in Sobolev spaces.  Appl. Math. Optim., 72(1):77–100, 2015.
  * [34] W. Grecksch and H. Lisei.  Approximation of stochastic nonlinear equations of Schrödinger type by the splitting method.  Stoch. Anal. Appl., 31(2):314–335, 2013.
  * [35] I. Gyöngy.  Lattice approximations for stochastic quasi-linear parabolic partial differential equations driven by space-time white noise. I.  Potential Anal., 9(1):1–25, 1998.
  * [36] I. Gyöngy.  Lattice approximations for stochastic quasi-linear parabolic partial differential equations driven by space-time white noise. II.  Potential Anal., 11(1):1–37, 1999.
  * [37] I. Gyöngy and A. Millet.  Accelerated finite elements schemes for parabolic stochastic partial differential equations.  Stoch. Partial Differ. Equ. Anal. Comput., 8(3):580–624, 2020.
  * [38] I. Gyöngy and D. Nualart.  Implicit scheme for quasi-linear parabolic partial differential equations perturbed by space-time white noise.  Stochastic Process. Appl., 58(1):57–72, 1995.
  * [39] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration, volume 31 of Springer Series in Computational Mathematics.  Springer, Heidelberg, 2010.  Structure-preserving algorithms for ordinary differential equations, Reprint of the second (2006) edition.
  * [40] N. Halidias.  Construction of positivity preserving numerical schemes for some multidimensional stochastic differential equations.  Discrete Contin. Dyn. Syst. Ser. B, 20(1):153–160, 2015.
  * [41] N. Halidias and I. S. Stamatiou.  Boundary preserving explicit scheme for the Aït-Sahalia model.  Discrete Contin. Dyn. Syst. Ser. B, 28(1):648–664, 2023.
  * [42] E. Hausenblas.  Approximation for semilinear stochastic evolution equations.  Potential Anal., 18(2):141–186, 2003.
  * [43] M. Hutzenthaler and A. Jentzen.  Numerical approximations of stochastic differential equations with non-globally Lipschitz continuous coefficients.  Mem. Amer. Math. Soc., 236(1112):v+99, 2015.
  * [44] A. Jentzen and P. E. Kloeden.  The numerical approximation of stochastic partial differential equations.  Milan J. Math., 77:205–244, 2009.
  * [45] A. Jentzen and P. E. Kloeden.  Overcoming the order barrier in the numerical approximation of stochastic partial differential equations with additive space-time noise.  Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 465(2102):649–667, 2009.
  * [46] D. Khoshnevisan.  Analysis of stochastic partial differential equations, volume 119 of CBMS Regional Conference Series in Mathematics.  Published for the Conference Board of the Mathematical Sciences, Washington, DC; by the American Mathematical Society, Providence, RI, 2014.
  * [47] M. Kovács, S. Larsson, and F. Lindgren.  Strong convergence of the finite element method with truncated noise for semilinear parabolic stochastic equations with additive noise.  Numer. Algorithms, 53(2-3):309–320, 2010.
  * [48] R. Kruse.  Strong and weak approximation of semilinear stochastic evolution equations, volume 2093 of Lecture Notes in Mathematics.  Springer, Cham, 2014.
  * [49] Z. Lei, S. Gan, and Z. Chen.  Strong and weak convergence rates of logarithmic transformed truncated EM methods for SDEs with positive solutions.  J. Comput. Appl. Math., 419:Paper No. 114758, 21, 2023.
  * [50] Y. Li, C.-W. Shu, and S. Tang.  A local discontinuous Galerkin method for nonlinear parabolic SPDEs.  ESAIM Math. Model. Numer. Anal., 55(suppl.):S187–S223, 2021.
  * [51] J. Liu.  A mass-preserving splitting scheme for the stochastic Schrödinger equation with multiplicative noise.  IMA J. Numer. Anal., 33(4):1469–1479, 2013.
  * [52] Z. Liu and Z. Qiao.  Strong approximation of monotone stochastic partial differential equations driven by white noise.  IMA J. Numer. Anal., 40(2):1074–1093, 2020.
  * [53] G. J. Lord, C. E. Powell, and T. Shardlow.  An introduction to computational stochastic PDEs.  Cambridge Texts in Applied Mathematics. Cambridge University Press, New York, 2014.
  * [54] G. J. Lord and T. Shardlow.  Postprocessing for stochastic parabolic partial differential equations.  SIAM J. Numer. Anal., 45(2):870–889, 2007.
  * [55] G. J. Lord and A. Tambue.  Stochastic exponential integrators for the finite element discretization of SPDEs for multiplicative and additive noise.  IMA J. Numer. Anal., 33(2):515–543, 2013.
  * [56] G. J. Lord and A. Tambue.  Stochastic exponential integrators for a finite element discretisation of SPDEs with additive noise.  Appl. Numer. Math., 136:163–182, 2019.
  * [57] X. Mao, F. Wei, and T. Wiriyakraikul.  Positivity preserving truncated Euler-Maruyama method for stochastic Lotka-Volterra competition model.  J. Comput. Appl. Math., 394:Paper No. 113566, 17, 2021.
  * [58] R. Marty.  On a splitting scheme for the nonlinear Schrödinger equation in a random medium.  Commun. Math. Sci., 4(4):679–705, 2006.
  * [59] R. Marty.  Local error of a splitting scheme for a nonlinear Schrödinger-type equation with random dispersion.  Commun. Math. Sci., 19(4):1051–1069, 2021.
  * [60] R. I. McLachlan and G. R. W. Quispel.  Splitting methods.  Acta Numer., 11:341–434, 2002.
  * [61] A. Millet and P.-L. Morien.  On implicit and explicit discretization schemes for parabolic SPDEs in any dimension.  Stochastic Process. Appl., 115(7):1073–1106, 2005.
  * [62] G. R. Moreno Flores.  On the (strict) positivity of solutions of the stochastic heat equation.  Ann. Probab., 42(4):1635–1643, 2014.
  * [63] C. Mueller.  On the support of solutions to the heat equation with noise.  Stochastics Stochastics Rep., 37(4):225–245, 1991.
  * [64] J. D. Mukam and A. Tambue.  A note on exponential Rosenbrock-Euler method for the finite element discretization of a semilinear parabolic partial differential equation.  Comput. Math. Appl., 76(7):1719–1738, 2018.
  * [65] J. D. Mukam and A. Tambue.  Strong convergence analysis of the stochastic exponential Rosenbrock scheme for the finite element discretization of semilinear SPDEs driven by multiplicative and additive noise.  J. Sci. Comput., 74(2):937–978, 2018.
  * [66] T. Müller-Gronbach and K. Ritter.  An implicit Euler scheme with non-uniform time discretization for heat equations with multiplicative noise.  BIT, 47(2):393–418, 2007.
  * [67] J. L. Padgett and Q. Sheng.  Convergence of an operator splitting scheme for abstract stochastic evolution equations.  In Advances in mathematical methods and high performance computing, volume 41 of Adv. Mech. Math., pages 163–179. Springer, Cham, 2019.
  * [68] R. Pettersson and M. Signahl.  Numerical approximation for a white noise driven SPDE with locally bounded drift.  Potential Anal., 22(4):375–393, 2005.
  * [69] J. Printems.  On the discretization in time of parabolic stochastic partial differential equations.  M2AN Math. Model. Numer. Anal., 35(6):1055–1078, 2001.
  * [70] L. Ryzhik.  Lecture notes for Introduction to Spde, spring 2016, May 2016.
  * [71] C. Scalone.  Positivity preserving stochastic $\theta$-methods for selected SDEs.  Appl. Numer. Math., 172:351–358, 2022.
  * [72] H. Schurz.  Basic concepts of numerical analysis of stochastic differential equations explained by balanced implicit theta methods.  In Stochastic differential equations and processes, volume 7 of Springer Proc. Math., pages 1–139. Springer, Heidelberg, 2012.
  * [73] T. Shardlow.  Numerical methods for stochastic parabolic PDEs.  Numer. Funct. Anal. Optim., 20(1-2):121–145, 1999.
  * [74] T. Shiga.  Two contrasting properties of solutions for one-dimensional stochastic partial differential equations.  Canad. J. Math., 46(2):415–437, 1994.
  * [75] L. Szpruch, X. Mao, D. J. Higham, and J. Pan.  Numerical simulation of a strongly nonlinear Ait-Sahalia-type interest rate model.  BIT, 51(2):405–425, 2011.
  * [76] G. Tessitore and J. Zabczyk.  Strict positivity for stochastic heat equations.  Stochastic Process. Appl., 77(1):83–98, 1998.
  * [77] I. Tubikanec, M. Tamborrino, P. Lansky, and E. Buckwar.  Qualitative properties of different numerical methods for the inhomogeneous geometric Brownian motion.  J. Comput. Appl. Math., 406:Paper No. 113951, 29, 2022.
  * [78] C. von Hallern and A. Rößler.  An analysis of the Milstein scheme for SPDEs without a commutative noise condition.  In Monte Carlo and quasi-Monte Carlo methods, volume 324 of Springer Proc. Math. Stat., pages 503–521. Springer, Cham, 2020.
  * [79] J. B. Walsh.  An introduction to stochastic partial differential equations.  In École d’été de probabilités de Saint-Flour, XIV—1984, volume 1180 of Lecture Notes in Math., pages 265–439. Springer, Berlin, 1986.
  * [80] J. B. Walsh.  Finite element methods for parabolic stochastic PDE’s.  Potential Anal., 23(1):1–43, 2005.
  * [81] X. Wang.  Strong convergence rates of the linear implicit Euler method for the finite element discretization of SPDEs with additive noise.  IMA J. Numer. Anal., 37(2):965–984, 2017.
  * [82] X. Wang and S. Gan.  A Runge-Kutta type scheme for nonlinear stochastic partial differential equations with multiplicative trace class noise.  Numer. Algorithms, 62(2):193–223, 2013.
  * [83] X. Wang and R. Qi.  A note on an accelerated exponential Euler method for parabolic SPDEs with additive noise.  Appl. Math. Lett., 46:31–37, 2015.
  * [84] Y. Yan.  Semidiscrete Galerkin approximation for a linear stochastic parabolic partial differential equation driven by an additive noise.  BIT, 44(4):829–847, 2004.
  * [85] Y. Yan.  Galerkin finite element methods for stochastic parabolic partial differential equations.  SIAM J. Numer. Anal., 43(4):1363–1384, 2005.
  * [86] X. Yang, Z. Yang, and C. Zhang.  Stochastic heat equation: numerical positivity and almost surely exponential stability.  Comput. Math. Appl., 119:312–318, 2022.