Numerical contractivity preserving implicit balanced Milstein-type schemes for SDEs with non-global Lipschitz coefficients

1.   Introduction

Assume $(\Omega, \mathcal{F}, \{\mathcal{F}_{t}\}_{t\geq0}, \mathbb{P})$ is a complete probability space with an increasing filtration $\{\mathcal{F}_{t}\}_{t\geq0}$ satisfying the usual conditions (that is, it is right continuous and increasing while $\mathcal{F}_{0}$ contains all $\mathbb{P}$-null sets). Let $\langle\cdot, \cdot\rangle$ denote the Euclidean inner product and $\Vert \cdot \Vert$ be the corresponding Euclidean vector norm in $\mathbb{R}^{n}$. The trace norm of a matrix $A\in\mathbb{R}^{n\times n}$, $n\in \mathbb{N}$ is denoted by $\Vert A \Vert : = \sqrt{trace(A^{T}A)}$. $\mathbb{E}$ denotes mathematical expectation. In this paper, we consider the following nonlinear systems of stochastic differential equations (SDEs) of $n$-dimensional Itô type given by

where $\mathbb{E}\Vert X_{0}\Vert^{2} < \infty$, $X(t)\in \mathbb{R}^{n}$, $W(t)$ is a scalar Brownian motion and the drift coefficient $f$ and diffusion coefficient g are Borel measurable real-valued vector functions in $\mathbb{R}^{n}$. Generally, analytical solutions to nonlinear SDEs (1.1) are seldom available, and resorting to numerical schemes for approximating SDEs are of significant interest in practice. Various SDEs arising from the field of applied science ^[1,2,3] rarely satisfy the restrictive global Lipschitz condition, such as, the stochastic Ginzburg Landau equation with a cubic nonlinear drift coefficient $f(x) = -4x-3x^{3}$, $x \in \mathbb{R}$. Unfortunately, the well-known Euler-Maruyama scheme generates divergent numerical approximations for SDEs with super-linearly growing coefficients ^[4]. Therefore, in order to avoid the numeric divergent phenomenon, numerous implicit schemes ^{[5,6,7,8,9,10,11,12,13,14,15]} and modifications of explicit schemes ^{[16,17,18,19,20,21,22,23,24,25,26,27]} attracted more and more attention for their numerical analysis of SDEs under non-globally Lipschitz conditions.

What we focus on in this paper is investigating whether two kinds of implicit balanced Milstein-type schemes can inherit numerically the relevant property of the mean square contractivity for nonlinear SDEs (1.1) with non-globally Lipschitz coefficients. To this end, let us first introduce the following definition of $p$-th moment stability for the SDEs (1.1) ^[28,29]. Suppose $Y(t)$ is the exact solution of the SDEs (1.1) with initial value $X(0) = Y_{0}$, where $\mathbb{E}\Vert Y_{0}\Vert^{2} < \infty$.

Definition 1.1. ^[28,29] The analytical solution of the SDEs (1.1) is called to be $p$-th moment stable if $\exists \bar{K}\in \mathbb{R}$

with $p$-th moment stability constant $\bar{K}$.

It should be noted that we call the analytical solution of the SDEs (1.1) to be strict $p$-th moment contractive if the $p$-th moment stability inequality (1.2) holds for $\bar{K} < 0$. More specifically, if $p = 2$ and $\bar{K} < 0$, we speak of strict mean square contractivity ^[29], (or exponential mean-square contractivity ^{[30,31,32,33]}). In general, strict $p$-th moment contractivity represents that initial perturbations have no significant impact on the long-term dynamic behavior of the SDEs (1.1). The $p$-th moment stability of nonlinear SDEs with $p$-th moment monotone coefficients was systematically investigated by Schurz in Lemma 2.8 ^[29]. The following Theorem gives a simplified overview of the nonlinear stability of the nonlinear SDEs (1.1) ^[29,34].

Theorem 1.2. ^[29,34] $X(t)$ and $Y(t)$ are analytical solutions of the SDEs (1.1) with different initial values $X_{0}$ and $Y_{0}$, respectively. Suppose that the drift and diffusion coefficients $f, g \in \rm{C}^{1}(\mathbb{R}^n)$ satisfy a respective global one-side Lipschitz condition and a global Lipschitz condition, i.e., there exists constants $\mu \in \mathbb{R}$ and $L > 0$, such that for $\forall X, Y \in \mathbb{R}^{n}$,

and

For $\forall t\in [0, +\infty)$,

where $a = 2\mu + L$.

The existence and uniqueness of the global solution to the SDEs (1.1) can be guaranteed ^[35,36]. Under the conditions of Theorem 1.2 and supposing $f (0) = 0$ and $g(0) = 0$, then for $\forall t\in [0, +\infty)$, $\mathbb{E}\Vert X(t)\Vert^{2}\leq e^{a t}\mathbb{E}\Vert X_{0}\Vert^{2}.$ Noting that when the diffusion coefficient $g = 0$, the SDEs (1.1) reduces to the corresponding deterministic ordinary differential equations (ODEs)

For any two solutions $X(t)$ and $Y(t)$ of ODEs (1.6) with initial data $X_{0}$ and $Y_{0}$, respectively, if the one-sided Lipschitz condition (1.3) holds with negative one-sided Lipschtiz constant $\mu$, then we have the contractive inequality

Nonlinear stability has been a central concept of the qualitative theory of ODEs ^[37]. For the numerical counterpart of the nonlinear stability of ODEs satisfying one-sided Lipschitz condition with one-sided Lipschitz constant $\mu < 0$, Dahlquist ^[38] presented the concept of $G$-stability for linear multistep methods (LMMs) and one-leg methods, while Butcher ^[39] introduced the concept of $B$-stability for implicit Runge-Kutta methods. We refer to the monograph ^[37] for more details about the contractivity of numerical methods for ODEs satisfying one-sided Lipschitz condition.

Similarly, in the case of SDEs, the strict mean square contractivity inequality (1.5) with the parameter $a < 0$ means an exponential decay of the mean square deviation between two solutions $X(t)$ and $Y(t)$ of the SDEs (1.1) with different initial data $X_{0}$ and $Y_{0}$, respectively. The numerical counterpart of the mean square contractivity for numerical schemes, which is omitted here, is defined in a similar manner as that of the exact solutions of the nonlinear stochastic systems in Definition 1.1. For numerical analysis of the mean square contractivity, Higham, Mao and Stuart ^[34] studied the stability of the backward Euler and split-step backward Euler methods. Yao and Gan ^[40] investigated the mean square contractivity of the drift-implicit Milstein and double-implicit Milstein schemes for nonlinear monotone SDEs. Exponential mean-square contractivity property of the stochastic Runge-Kutta methods ^[41], stochastic $\theta$- methods ^[42] and stochastic linear multistep methods (mainly mentioned two-step methods) ^[43] were discussed; however, it was noteworthy that the mean value theorem was utilized in the proof of the stability theorem of the last two numerical schemes ^[42,43]. Moment stability analysis of the two-point motion of drift-implicit $\theta$-methods (including the backward Euler method ^[34] when $\theta = 1$) for SDEs was analyzed systemically by Henri ^[29]. The aim of this paper is to focus on investigating whether the drift implicit balanced Milstein (DIBM) scheme and the semi-implicit balanced Milstein (SIBM) scheme (or double-implicit balanced Milstein scheme), which are considered as the modifications of the drift-implicit Milstein scheme and the double-implicit Milstein scheme ^[40], respectively, can also possess numeric property of mean square contractivity. Therefore, the theorems in this paper can be identified as the extension of ^[40]. The rest of the paper is organized as follows. In the next section, two types of implicit balanced Milstein schemes are introduced. In Section 3, sufficient conditions for the mean square contractivity for both of the mentioned implicit balanced Milstein-type schemes are derived. In Section 4, numerical experiments are given to verify the theoretical results. The last section presents some conclusions.

2.   Numerical schemes

This section mainly witnesses two types of implicit balanced Milstein schemes, which will be investigated in the following sections. As the first implicit balanced numerical method, a class of balanced implicit (BI) methods was proposed by Milstein, Platen and Schurz ^[44] for solving stiff SDEs, namely,

at some grid points $t_{k} = kh$, $k = 0, 1, \ldots$, on the time interval $[0, +\infty)$ with time step $h$. $x_{0} = X_{0}$ and $\Delta W_{k} = W(t_{k+1})-W(t_{k})$ denote the increment of Brownian motion, $C_{k} = c_{0}(x_{k})h+c_{1}(x_{k})|\Delta W_{k}|$.

Kahl and Schurz ^[45] presented a class of balanced Milstein (BM) schemes, namely,

where $x_{0} = X_{0}$, $L^{1}g(x) = \frac{\partial g(x)}{\partial x}g(x)$ with the $j$-th component

$C_{k} = c_{0}(x_{k})h+c_{2}(x_{k})\left(|\Delta W_{k}|^{2} -h\right)$, where control functions $c_{0}$ and $c_{2}$ satisfy the following condition ^[46,47].

Assumption 2.1. The control functions $c_{0}$ and $c_{2}$ are bounded $n\times n$-matrix-valued functions. For any real numbers $\alpha_{0}\in[0, \tilde{\alpha}_{1}]$, $\alpha_{2}\in[-\tilde{\alpha}_{2}, \tilde{\alpha}_{2}]$, where $\tilde{\alpha}_{1}\geq h$, $\tilde{\alpha}_{2}\geq ||\Delta W_{k}|^{2}-h|$ for any step-size $h$ under consideration and $x\in \mathbb{R}^{n}$, the $n\times n$ matrix is

where $I$ denotes the $n\times n$ identity matrix, is invertible and there exists a positive constant $K$ satisfying

For the choice of the control functions $c_{0}$ and $c_{2}$, Alcock and Burrage ^[48] investigated the choice of optimal parameter for the BI method (2.1). Wang and Liu ^[46] presented three typical criterions for one-dimension case. In practical computation, the control functions $c_{0}$ and $c_{2}$ are, in general matrix, often chosen as constants satisfying Assumption 2.1 ^[47,48]. For simplicity, assume that the control functions $c_{0}$ and $c_{2}$ will be chosen as constant matrices, satisfying $c_{0}$ and $c_{2}$ as positive definite or $c_{0}-c_{2}$ and $c_{2}$ as positive semi-definite ^[44,45], which satisfy Assumption 2.1.

In the following, let us introduce two kinds of implicit balanced Milstein-type schemes, e.g., the DIBM scheme and the SIBM scheme. The DIBM approximation ^[49] applied to SDEs (1.1) has the following form

where $x_{0} = X_{0}$, $C_{k} = c_{0}(x_{k})h+c_{2}(x_{k})\left(|\Delta W_{k}|^{2}-h\right)$. In addition, noticing the term in (2.4)

and bringing partial implicitness to the DIBM scheme (2.4) leads to the following numerical method

where $x_{0} = X_{0}$, $C_{k} = c_{0}(y_{k})h+c_{2}(y_{k})\left(|\Delta W_{k}|^{2}-h\right)$, $k = 0, 1, \ldots$. This method is named the SIBM scheme (or double-implicit balanced Milstein scheme) ^[49].

3.   Numerical contractivity analysis

In this section, we aim to investigate the numerical counterpart of the mean square contractivity for the above-mentioned two types of implicit balanced Milstein schemes, e.g., the DIBM scheme (2.4) and the SIBM scheme (2.5). It is proved that both schemes can well replicate the mean square contractivity of the refered nonlinear systems (1.1).

Assumption 3.1. Suppose there exists a constant $\omega$, such that for $\forall X, Y \in \mathbb{R}^{n}$,

Theorem 3.2. Under the assumptions of Theorem 1.2 and (3.1), assume $2h\mu K < 1$. Write $\bar{a} = 2\mu+KL$, $c_{1} = \frac{2\mu+KL+\frac{1}{2}hK\omega}{1-2h\mu K}K$. Let $\{x_{k}\}_{k\in \mathbb{N}}$ and $\{y_{k}\}_{k\in \mathbb{N}}$ be two parallel approximation sequences obtained by the DIBM scheme (2.4) starting from two distinct initial data $X_{0}$ and $Y_{0}$, respectively, then

where $\bar{a} > 0$, $c_{1} > 0$ or $\bar{a}\leq0$, $0 < h\leq\frac{-2\bar{a}}{K\omega}$, $c_{1}\leq0$.

Proof. By the DIBM scheme (2.4), we have

which leads to

Therefore, squaring both sides of the above equality yields

Taking expectation and using the one-side Lipschitz condition (1.3), the global Lipschitz condition (1.4) and inequality (3.1), we obtain

which yields

Consequently, taking account of the fact that $2h\mu K < 1$, we obtain

(ⅰ) If $\bar{a} = 2\mu+KL > 0$, we have $\frac{1+hLK^{2}+\frac{1}{2}h^{2}\omega K^{2}}{1-2h\mu K} > 1$ and

where $c_{1} = \frac{2\mu+KL+\frac{1}{2}hK\omega}{1-2h\mu K}K > 0$.

(ⅱ) If $\bar{a}\leq 0$, we have $0 < \frac{1+hLK^{2}+\frac{1}{2}h^{2}\omega K^{2}}{1-2h\mu K}\leq1$ for $0 < h\leq\frac{-2\bar{a}}{K\omega}$,

where $c_{1}\leq0$. □

Corollary 3.3. Under the assumptions of Theorem 3.2 and $f (0) = g(0) = 0$, then

where $\bar{a} > 0$, $c_{1} > 0$ or $\bar{a}\leq0$, $0 < h\leq\frac{-2\bar{a}}{K\omega}$, $c_{1}\leq0$.

Note that the inequality (3.3), which can be regarded as numerical analogue of the mean square stability inequality (1.5) for the analytic solutions of the SDEs (1.1), means that the DIBM scheme (2.4) is mean square stable. Specifically, when $\bar{a} = 2\mu+KL < 0$ and $0 < h < \frac{-2\bar{a}}{K\omega}$, inequality (3.3) represents the strict mean square contractivity of the DIBM scheme (2.4), which means that any two numerical trajectories of the stochastic dynamical system (1.1) converge to one other in mean square at an exponential rate and that perturbations in the initial data have no significant impact on numerical dynamic behavior along the entire time-scale $[0, +\infty)$. For strict mean square contractive approximation sequences $\{x_{k}\}_{k\in \mathbb{N}}$ and $\{y_{k}\}_{k\in \mathbb{N}}$, we have $\lim\limits_{k \to +\infty }{\mathbb{E}\Vert x_{k}-y_{k}\Vert^{2}} = 0$.

Let us give a minute for the description of the following theorem, which sheds light on the mean square contractivity of the SIBM scheme (2.5).

Theorem 3.4. Under the assumptions of Theorem 1.2 and Assumptions 3.1, suppose $2h\mu K < 1$. Let $\tilde{a} = \frac{1}{2}+K \left[ 2\mu + K\left(L +\frac{1}{2}\omega\right) \right]$, $c_{2} = \frac{\frac{1}{2}+K \left[ 2\mu + K\left(L +\frac{1}{2}\omega+\frac{3}{4}h\omega\right) \right]}{1-2h\mu K}$, then numerical solutions $\{x_{k}\}_{k\in \mathbb{N}}$ and $\{y_{k}\}_{k\in \mathbb{N}}$ with distinct initial values $X_{0}$ and $Y_{0}$, respectively, obtained by the SIBM scheme (2.5) satisfy

where $\tilde{a} > 0$, $c_{2} > 0$ or $\tilde{a}\leq0$, $0 < h\leq\frac{-4\tilde{a}}{3K^{2}\omega}$, $c_{2}\leq0$.

Proof. By the SIBM scheme (2.5), we have

Squaring both sides of the above equality yields

Utilizing the one-side Lipschitz condition (1.3), the global Lipschitz condition (1.4) and inequalities (3.1) and (3.2), we obtain

which leads to

Because $2h\mu K < 1$, consequently,

Let us consider two possible cases:

(ⅰ) If $\tilde{a} = \frac{1}{2}+K \left[ 2\mu + K\left(L +\frac{1}{2}\omega\right) \right] > 0$, we get $\frac {1+\frac{1}{2}h +hK^{2}\left(L+\frac{1}{2}\omega+\frac{3}{4}h\omega \right)} {1-2h\mu K} > 1$ and

where $c_{2} = \frac{\frac{1}{2}+K \left[ 2\mu + K\left(L +\frac{1}{2}\omega+\frac{3}{4}h\omega\right) \right]}{1-2h\mu K} > 0$.

(ⅱ) If $\tilde{a} \leq 0$, we have $0 < \frac {1+\frac{1}{2}h +hK^{2}\left(c+\frac{1}{2}\omega+\frac{3}{4}h\omega \right)} {1-2h\mu K}\leq1$ for $0 < h\leq\frac{-4\tilde{a}}{3K^{2}\omega}$ and

where $c_{2}\leq0$. □

Corollary 3.5. Under assumptions of Theorem 3.4 and $f (0) = 0$, $g(0) = 0$, then

where $\tilde{a} > 0$, $c_{2} > 0$ or $\tilde{a}\leq0$, $0 < h\leq\frac{-4\tilde{a}}{3K^{2}\omega}$, $c_{2}\leq0$.

The inequality (3.4) is indicative of the mean square stability of the SIBM scheme (2.5). Specifically, when $0 < h < \frac{-4\tilde{a}}{3K^{2}\omega}$ and $\tilde{a} = \frac{1}{2}+K \left[ 2\mu + K\left(L +\frac{1}{2}\omega\right) \right] < 0$, the inequality (3.4) manifests the strict mean square contractivity of the SIBM scheme (2.5). In this situation, we can easily have that $\lim\limits_{k \to +\infty }{\mathbb{E}\Vert x_{k}-y_{k}\Vert^{2}} = 0$.

4.   Numerical results

In this section, we illustrate intuitively the given theoretical analysis obtained in previous sections through numerical examples. Let us first consider the following one-dimension stochastic Ginzburg Landua equation with a cubic nonlinearity in the drift and linear diffusion ^{[32,41,43,50,51]}:

with different initial values $X_{0} = 0$ and $Y_{0} = 1$. The cubic drift coefficient is $f(x(t)) = Ax(t) + Bx^3(t)$ and the linear diffusion function is $g(x(t)) = Cx(t)$. We choose constants $A = -4$, $B = -3$, $C = 1$, $c_{0} = 4$ and $c_{2} = 1$. Clearly, the one-side Lipschitz condition (1.3) and global Lipschitz condition (1.4) hold with $\mu = -4$ and $L = 1$, and the problem (4.1) is strict mean square contractive with $a = -7$, according to Theorem 1.2. As shown in Figure 1, where the long-time development and evolution of the mean square deviation $\mathbb{E}\Vert x_{k}-y_{k}\Vert^{2}$ in logarithmic scale is depicted even for quite large step sizes $h = 1, 2, 5$ and $10$, both of the DIBM scheme (2.4) and the SIBM scheme (2.5) can well reproduce strict mean square contractivity. It is well consistent with the theoretical results established in Theorems 3.2 and 3.4.

As a second example, consider the Itô SDE with nonlinear diffusion ^[29,41,43]:

with distinct initial data $X_{0} = 0$ and $Y_{0} = 1$. The linear drift coefficient is $f(x(t)) = Ax(t)$ and the nonlinear diffusion function is $g(x(t)) = B\sin(x(t))$. We choose constants $A = -1$ and $B = 1$. According to Theorem 1.2, the problem (4.2) is strict mean square contractive with $\mu = -1$, $L = 1$ and $a = -1$. The numerical results, shown in Figure 2, confirm the validity of theoretical conclusions in the previous sections.

5.   Conclusions

Two types of implicit balanced Milstein schemes, e.g., the DIBM scheme and the SIBM scheme, were utilized to simulate the nonlinear SDEs (1.1) with non-global Lipschitz coefficients. We have systematically analyzed the numerical counterpart of mean square contractivity of the implicit balanced Milstein-type schemes for the underlying SDEs (1.1) under the assumptions of Theorem 1.2 and Assumptions 3.1. It was proved that both schemes considered can successfully inherit the property of mean square contractivity. Numerical experiments conformed to the theoretical results obtained in this paper.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research is supported by the Scientific Research Foundation of Hunan Provincial Education Department (19C0146, 21C1641) and the Scientific Research Project of Changsha Normal University (XJYB202338).

Conflict of interest

All authors declare that there are no competing interests.