Backward Euler–Maruyama Method for the Random Periodic Solution of a Stochastic Differential Equation with a Monotone Drift

In this paper, we study the existence and uniqueness of the random periodic solution for a stochastic differential equation with a one-sided Lipschitz condition (also known as monotonicity condition) and the convergence of its numerical approximation via the backward Euler–Maruyama method. The existence of the random periodic solution is shown as the limit of the pull-back flows of the SDE and the discretized SDE, respectively. We establish a convergence rate of the strong error for the backward Euler–Maruyama method with order of convergence 1/2.


Introduction
Periodicity is widely exhibited in a large number of natural phenomena like oscillations, waves, or even lying behind many complicated ensembles such as biological and economic systems. However, periodic behaviors are often found to be subject to random perturbation or under the influence of noise. Physicists have attempted to study random perturbations to periodic solutions for some time by considering a first linear approximation or asymptotic expansions in small noise regime, but this approach restricted its applicability to the small fluctuation (c.f. Van Kampen [13], Weiss and Knoblock [16]). It was only until recently that the random periodic solution was endowed with a proper definition (c.f. Zhao and Zheng [19], Feng, Zhao and Zhou [8]), which is compatible with definitions of both the stationary solution (also termed as random fixed points) and the deterministic periodic solution. It gives a rigorous and clearer understanding to physically interesting problems of certain random phenomena with a periodic nature and also represents a long time limit of the underlying random dynamical system.
Let us recall the definition of the random periodic solution for stochastic semiflows given in [8]. Let H be a separable Banach space. Denote by ( , F, P, (θ s ) s∈R ) a metric dynamical system and θ s : → is assumed to be measurably invertible for all s ∈ R. Denote := {(t, s) ∈ R 2 , s ≤ t}. Consider a stochastic semi-flow u : × × H → H , which satisfies the following standard condition for all r ≤ s ≤ t, r , s, t ∈ R, for a.e. ω ∈ . We do not assume the map u(t, s, ω) : H → H to be invertible for (t, s) ∈ , ω ∈ .
Building on this new concept, there have been more recent progresses toward understanding the random periodicity of various stochastic systems. The existence of random periodic solutions to stochastic differential equations (SDEs) and stochastic partial differential equations (SPDEs) are initially studied in [8] and [4], with additive noise. Instead of following the traditional geometric method of establishing the Poincaré mapping, a new analytical method for coupled infinite horizon forward-backward integral equations is introduced. It was then followed by the study on the anticipating random periodic solutions (c.f. Feng, Wu and Zhao: [6] and [7]). Regarding applications, Chekroun, Simonnet and Ghil [3] employed random periodic results to climate dynamics, and Wang [14] observed random periodicity behavior in the study of bifurcations of stochastic reaction diffusion equations.
In general, random periodic solutions cannot be solved explicitly. One may treat the numerical approximation that stay sufficient close to the true solution as a good substitute to study stochastic dynamics. It is worth mentioning here that this is a numerical approximation of an infinite time horizon problem. The classical numerical approaches including the Euler-Marymaya method and a modified Milstein method to simulate random period solutions of a dissipative system with global Lipschitz condition have been investigated in [5], which is the first paper that numerical schemes were used to approximate the random period trajectory.
In this paper, we study the random periodic solutions of stochastic differential equations with weakened conditions on the drift term compared to [5] and simulate them via the backward Euler-Maruyama method. Let W : R × → R d be a standard two-sided Wiener process on the probability space ( , F, P), with the filtration defined by F t Throughout this paper, we shall use | · | for the Euclidean norm, u := E[|u| 2 ] and u p : We are interested in the R d -valued random periodic solution to a SDE of the form where ξ is a F t 0 -measurable random initial condition. In addition, A, f , and g, and ξ satisfy the following assumptions: The linear operator A : R d → R d is self-adjoint and positive definite.
Assumption 1 implies the existence of a positive, increasing sequence (λ i ) i∈ [d] ⊂ R such that 0 < λ 1 ≤ λ 2 ≤ · · · λ d , and of an orthonormal basis (e i ) i∈ [d] of R d such that Assumption 2 The mapping f : R × R d → R d is continuous and periodic in time with period τ . Moreover, there exists a C f ∈ (0, ∞) such that for all u, u 1 , u 2 ∈ R d and t ∈ [0, τ ).
It is well known that under these assumptions the solution X t 0 · : [t 0 , T ] × → R d to (3) is uniquely determined by the variation-of-constants formula

The Pull-Back
We know there exists a standard P-preserving ergodic Wiener shift θ such that θ t (ω)(s) = W t+s − W t for s, t ∈ R. We will show that when k → ∞, the pullback X −kτ t (ξ ) has a limit X * t in L 2 ( ) and X * t is the random periodic solution of SDE (3), satisfying To achieve it, we need additional assumptions on ξ and f .
Assumption 5 There exists a constant C ξ such that ξ < C ξ .

Assumption 6
There exists a constantĈ f such that Assumption 6 together with Assumptions 1 to 3 ensures the existence of a global semiflow generated from SDE (3) with additive noise [11]. Section 3 is devoted to the first main result, which claims the existence and uniqueness of random periodic solutions to the SDE (3) under the one-sided Lipschitz condition on the drift.

Theorem 7
Under Assumptions 1 to 6, there exists a unique random periodic solution In Sect. 4, we derive additional properties of the solution such as the uniform boundedness for a higher moment of X −kτ t and solution regularity under an additional Assumption 13, which imposes superlinearity of f and assumes a larger lowerbound for λ 1 compared to Assumption 4. Those properties will play an important role in proving the order of convergence of the backward Euler-Maruyama in Theorem 19.

The Backward Euler-Maruyama
For stiff ordinary differential equations, the implicit method is preferred due to its good performance even on a time grid with a large step size [15]. For its stochastic counterpart such as (3), we shall approximate the solution using the backward Euler-Maruyama method, the simplest version of implicit methods for SDEs.
Let us fix an equidistant partition T h := { jh, j ∈ Z} with stepsize h ∈ (0, 1). Note that T h stretches along the real line because eventually we are dealing with an infinite time horizon problem in the form of (5). Then to simulate the solution to (3) starting at −kτ , the backward Euler-Maruyama method on T h is given by the recursion for all j ∈ N, where the initial valueX −kτ −kτ = ξ , and W −kτ + jh : , and similar arguments for the g term. The implementation of (7) requires solving a nonlinear equation at each iteration. Theorem 16 ensures the well-posedness of difference equation (7) under Assumptions 1 to 4. We explore the random periodicity of its solution in Sect. 5 and prove the second main result in our paper: Theorem 8 Under Assumptions 1 to rm 5, for any h ∈ (0, 1) with τ = nh, n ∈ N, the backward Euler-Maruyama method (7) admits a random period solution on T h .
We also determine a strong order 1/2 for the backward Euler-Maruyama method in Theorem 19 and Corollary 20. Compared to Theorem 3.4 and Theorem 4.2 in [5] which imposed condition on the size of h (to be sufficient small) because of the implementation of explicit numerical methods, we benefit a flexible choice of stepsize h from using the backward Euler-Maruyama method even in the infinite horizon case.
Finally we assess the performance of the backward Euler-Maruyama method via a numerical experiment and compare it with the one of the classical Euler-Maruyama method under various steps. The result shows that the backward Euler-Maruyama method is able to converge to the random periodic solution when the stepsize is fairly large while Euler-Maruyama method diverges.

Preliminaries
In this section, we present a few useful mathematical tools for later use. then If in addition, the function a is non-decreasing, then Lemma 10 (The Grönwall inequality: a discrete version [17,18]) Consider two nonnegative sequences (u n ) n∈N , (a n ) n∈N ⊂ R which for some given w ∈ [0, ∞) satisfy u n ≤ a n + w Then, for all n ∈ N, it also holds true that Also the crucial but simple equality for the analysis of the backward Euler-Maruyama is

Existence and Uniqueness of the Random Periodic Solution
In this section, we focus on the existence and uniqueness of the random periodic solution to SDE (3). To achieve it, we first show there is a uniform bound for the second moment of its solution under necessary assumptions.

Lemma 11
For SDE (3) with given initial condition ξ and satisfying Assumptions 1 to 5, we have where K 2 :=

Proof of Lemma 11
Applying Itô formula to e 2λ 1 t X −kτ t (ξ ) 2 and taking the expectation yield Note that 2(λ 1 I − A) is non-positive definite. Then making use of Assumptions 2 and 3 gives and K 3 := 2C f . Note that K 3 ≤ 2λ 1 because of Assumption 4. By the Grönwall inequality, we have that Note that K 1 e 2λ 1 kτ + K 2 = ξ 2 . By Assumption 5, it leads to Then we explore the solution dependence on initial conditions.
In addition, if Assumption 4 holds, then for every > 0, there exists a t ≥ −kτ such that it holds whenevert ≥ t.

Proof of Lemma 12
. From (4), we have that Similar as the proof of Lemma 11, we apply Itô formula to e 2λ 1 t |E −kτ t | 2 , take the expectation, make use of Assumption 2 and get Applying Eqn. (10) gives the desired inequality. The claim in (14) follows if Assumption 4 holds.
With Lemmas 11, 12 and Assumption 6, the main result Theorem 7 can be shown by following the same argument in the proof of Theorem 2.4 in [5].

More Results on the Solution
In this section, we explore some properties of the solution to 3 for analysis later.

Assumption 13
There exists a constant q ∈ (1, ∞) and a positive L such that for t 1 , t 2 ∈ [0, τ ) and u 1 , u 2 ∈ R d . In addition, there exists a positive number p ∈ [4q − 2, ∞) such that The first property we will show is the uniform boundedness for the p-th moment of the SDE solution.

Proof of Proposition 14
From the proof of Lemma 11, we know that Then applying Itô formula to e pλ 1 t |X −kτ t | p = e 2λ 1 t |X −kτ t | 2 p/2 and taking into consideration 2(λ 1 I − A) being non-positive definite give

Now by the Young inequality
and the inequality from fundamental calculus, Following a similar argument as in Proposition 5.4 and 5.5 [2], we can easily get the following bounds for analysis later.

Proposition 15 Let Assumptions 1 to 5 and 13 hold. Then there exists a positive constant C q,A, f which depends on q, d, A,C f only, such that
for all t 1 , t 2 ≥ −kτ . Moreover,

The Random Periodic Solution of the Backward Euler-Maruyama Scheme
In this section, we will prove that the backward Euler-Maruyama method (7) admits a unique discretized random period solution. To achieve this, let us first show the existence and uniqueness of solution to the targeted scheme.
The next Lemma claims there is a uniform bound for the second moment of the numerical solution under necessary assumptions.

Lemma 17
Under Assumptions 1 to 5, for any h ∈ (0, 1), it holds for the backward Euler-Maruyama method (7) on T h that Proof of Lemma 17 First note that from (11), we have that for any N ∈ N From (7), we have that Note that E g (N − 1)h W −kτ +(N −1)h ,X −kτ −kτ +(N −1)h = 0. Taking the expectation of both sides of (22) and making use of Assumption 2 give Then cancelling the same term on both side gives Let α := 2C f +σ 2 2(λ 1 −C f ) . Rearranging the terms above gives By iteration, this leads to Because of Assumptions 4 and 5, the term on the right-hand side above can be bounded by ξ 2 + α, which is independent of k, N and h.
The next result shows two numerical solutions starting from different initial conditions can be arbitrarily close after sufficiently many iterations.  (11) again, which allows us to examine the following term: Following a similar argument as in the proof of Lemma 17, this leads to By iteration, we have Because of λ 1 > C f , the assertion follows.

Error Analysis
Theorem 19 Under Assumptions 1 to 5 and 13, for any h ∈ (0, 1) with τ = nh, n ∈ N, there exists a constant C that depends on q, A, f , g and d such that the backward Euler-Maruyama method (7) approximates the true solution of (3) on T h with

Proof of Theorem 19
First note that Define By the Young's inequality and Assumption 2, we are able to choose 2 By Proposition 15, we know there exists a constant C depending on q, A, f and g such that Note that β is bounded because of Proposition 15. Then from (11) and the estimate above, we have that The inequality above can be rearranged to By iteration and assumingX −kτ −kτ = X −kτ −kτ = ξ , we have Finally due to Assumption 4 (alternatively, Assumption 13), we have e N 2 ≤ βh (λ 1 −C f ) 2 . Then the assertion follows. The pull-back patĥ X −30 (t, θ −t ω) generated by backward Euler-Maruyama method evaluation ofX 0 (r , θ −r ω). To allow convergence, we look at the path pattern from t = 2 to t = 5 in Fig. 3. Apparently we have obtained a periodic pull-back path as expected, which in turn shows the random periodicity of the original path. Finally, we test the order of convergence of the backward Euler-Maruyama method and compare the performance with (forward) Euler-Maruyama method. For its approximation, we first generated a reference solution with a small step size of h ref = 2 −15 . This reference solution was then compared to numerical solutions with larger step  Fig. 4. We plot the Monte Carlo estimates of the root-mean-squared errors versus the underlying temporal step size, i.e., the number i on the x-axis indicates the corresponding simulation is based on the temporal step size h = 2 −i . Both methods give the order of convergence above 1, which is beyond the theoretical order of convergence. When the stepsize is large, say, h = 2 −4 , the Euler-Maruyama method has the error 0.048 which is almost five times of the error 0.011 from the backward Euler-Maruyama method. Indeed if we relax the stepsize to h = 2 −3 , the Euler-Maruyama diverges while the backward Euler-Maruyama method still converges as expected. This further supports Theorem 8 and the advantage of backward Euler-Maruyama method: the backward Euler-Maruyama method converges regardless the size of stepsize (h < 1).
Funding This work is supported by the Alan Turing Institute for funding this work under EPSRC Grant EP/N510129/1 and EPSRC for funding though the Project EP/S026347/1, titled "Unparameterised multimodal data, high order signatures, and the mathematics of data science."

Data Availability
The datasets generated during and/or analyzed during the current study are available in the Github repository, https://github.com/yuewu57/RPS_BackwardEuler.

Declarations
Competing interests Authors are required to disclose financial or non-financial interests that are directly or indirectly related to the work submitted for publication.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.