A note on Asymptotic mean-square stability of stochastic linear two-step methods for SDEs

TL;DR

This paper analyzes the asymptotic mean-square stability of various two-step stochastic numerical schemes for SDEs, providing conditions for stability and comparing their performance in scalar and multi-dimensional cases.

Contribution

It derives necessary and sufficient stability conditions for two-step schemes applied to SDEs, including scalar and systems, and compares stability properties of different methods.

Findings

BDF2 scheme is unconditionally stable for any step-size.

AB2 and AM2 methods are unconditionally stable.

Numerical experiments confirm theoretical stability conditions.

Abstract

In this note we study the asymptotic mean-square stability for two-step schemes applied to a scalar stochastic differential equation (sde) and applied to systems of sdes. We derive necessary and sufficient conditions for the asymptotic MS-stability of the methods in terms of the parameters of the schemes. The stochastic Backward Differentiation Formula (BDF2) scheme is asymptotically mean-square stable for any step-size whereas the two-step Adams-Bashforth (AB2) and Adams-Moulton (AM2) methods are unconditionally stable. The improved versions of the schemes do not perform better w.r.t their stability behavior in the scalar case, as expected, but the situation is different in more dimensions. Numerical experiments confirm theoretical results.

Tables

Table 1. Table 1 . Coefficients of two-step Maruyama schemes and improved two-step schemes as in ( 2.2 ) and ( 2.3 ) with α 0 = γ 1 = 1 . subscript 𝛼 0 subscript 𝛾 1 1 \alpha_{0}=\gamma_{1}=1.

Method	${\alpha }_1$	${\alpha }_2$	${\beta }_0$	${\beta }_1$	${\beta }_2$	${\gamma }_2$	${\eta }_1$	${\eta }_2$
AB2	$-1$	$0$	$0$	$3/2$	-$1/2$	$0$	$-$	$-$
AB2I	$-1$	$0$	$0$	$3/2$	-$1/2$	$0$	$0$	$-1/2$
AM2	$-1$	$0$	$5/12$	$8/12$	-$1/12$	$0$	$-$	$-$
AM2I	$-1$	$0$	$5/12$	$8/12$	-$1/12$	$0$	$-5/12$	$-1/12$
BDF2	$-4/3$	$1/3$	$2/3$	$0$	$0$	$-1/3$	$-$	$-$
BDF2I	$-4/3$	$1/3$	$2/3$	$0$	$0$	$-1/3$	$-2/3$	$1/3$

Table 2. Table 2 . Parameters of two-step schemes as in ( 2.4 ).

Method	$a$	$b$	$c$	$d$
AB2	$1+\left(3/2\right)x$	$y$	$-x/2$	$0$
AB2I	$1+\left(3/2\right)x$	$y\left(1+x\right)$	$-x/2$	$-xy/2$
AM2	$\frac{1+\left(8/12\right)x}{1-\left(5/12\right)x}$	$\frac{y}{1-\left(5/12\right)x}$	$\frac{-x/12}{1-\left(5/12\right)x}$	0
AM2I	$\frac{1+\left(8/12\right)x}{1-\left(5/12\right)x}$	$\frac{y+\left(7/12\right)xy}{1-\left(5/12\right)x}$	$\frac{-x/12}{1-\left(5/12\right)x}$	$\frac{-xy/12}{1-\left(5/12\right)x}$
BDF2	$\frac{4/3}{1-\left(2/3\right)x}$	$\frac{y}{1-\left(2/3\right)x}$	$\frac{-1/3}{1-\left(2/3\right)x}$	$\frac{-y/3}{1-\left(2/3\right)x}$
BDF2I	$\frac{4/3}{1-\left(2/3\right)x}$	$\frac{y+xy/3}{1-\left(2/3\right)x}$	$\frac{-1/3}{1-\left(2/3\right)x}$	$\frac{-y/3}{1-\left(2/3\right)x}$

Equations

$$dX(t)=F(t,X(t))dt+G(t,X(t))dW(t),X(t_0)=X_0,$$

$$\sum _{j=0}^2{\alpha }_j{X}_{i-j}=h\sum _{j=0}^2{\beta }_j{F}_{i-j}+\sum _{j=1}^2\sum _{r=1}^m{\gamma }_j{G}_{r,i-j}\sqrt{h}{\xi }_{r,i-j},i=2,3,\dots ,$$

$$\sum _{j=0}^2{\alpha }_j{X}_{i-j}=h\sum _{j=0}^2{\beta }_j{F}_{i-j}+\sum _{j=1}^2\sum _{r=1}^m\left({\gamma }_j{G}_{r,i-j}\sqrt{h}{\xi }_{r,i-j}+({\gamma }_j+{\eta }_j)(F^{\prime }{G}_r{)}_{i-j}{h}^{3/2}{\xi }_{r,i-j}\right),$$

$$dX(t)=\lambda X(t)dt+\mu X(t)dW(t),X(t_0)=X_0,$$

$$\sum _{j=0}^2{\alpha }_j{X}_{i-j}=h\sum _{j=0}^2{\beta }_j\lambda X_{i-j}+\sum _{j=1}^2{\gamma }_j\mu X_{i-j}\sqrt{h}{\xi }_{i-j},i=2,3,\dots ,$$

$$\sum _{j=0}^2{\alpha }_j{X}_{i-j}=h\sum _{j=0}^2{\beta }_j\lambda X_{i-j}+\sum _{j=1}^2\left({\gamma }_j\mu X_{i-j}\sqrt{h}{\xi }_{i-j}+({\gamma }_j+{\eta }_j)\lambda \mu X_{i-j}{h}^{3/2}{\xi }_{i-j}\right),i=2,3,\dots ,$$

$$X_i=a{X}_{i-1}+c{X}_{i-2}+b{X}_{i-1}{\xi }_{i-1}+d{X}_{i-2}{\xi }_{i-2},i=2,3,\dots ,$$

$$a=\frac{-{\alpha }_1+h{\beta }_1\lambda }{{\alpha }_0-h{\beta }_0\lambda },b=\frac{\sqrt{h}{\gamma }_1\mu }{{\alpha }_0-h{\beta }_0\lambda },c=\frac{-{\alpha }_2+h{\beta }_2\lambda }{{\alpha }_0-h{\beta }_0\lambda },d=\frac{\sqrt{h}{\gamma }_2\mu }{{\alpha }_0-h{\beta }_0\lambda }$$

$$b^{\ast }=b+\frac{\lambda \mu ({\gamma }_1+{\eta }_1)h^{3/2}}{{\alpha }_0-h{\beta }_0\lambda },d^{\ast }=d+\frac{\lambda \mu ({\gamma }_2+{\eta }_2)h^{3/2}}{{\alpha }_0-h{\beta }_0\lambda }.$$

$$\mathcal{S}=\left[\begin{array}{cccc}|a{|}^2+|b{|}^2& \overline{a}c& a\overline{c}& |c{|}^2+|d{|}^2+ab\overline{d}+\overline{a}\overline{b}d\\ \overline{a}& 0& \overline{c}& b\overline{d}\\ a& c& 0& \overline{b}d\\ 1& 0& 0& 0\end{array}\right],$$

$$\rho (\mathcal{S})<1.$$

$$P(z)=z^4+p_1{z}^3+p_2{z}^2+p_{3}z+p_4$$

$$p_1=-|a{|}^2-|b{|}^2,p_2=-2|c{|}^2-|d{|}^2-2\Re (ab\overline{d})-2\Re (a^2\overline{c}),$$

$$p_3=-2\Re (\overline{a}bc\overline{d})-|a{|}^2|c{|}^2+|b{|}^2|c{|}^2,p_4=|c{|}^4+|c{|}^2|d{|}^2.$$

$$P^{\#}(z)=z^4\overline{P}\left(\frac{1}{\overline{z}}\right)=p_4{z}^4+p_3{z}^3+p_2{z}^2+p_{1}z+1;$$

$${\Delta }_4(P,P^{\#})=\left[\begin{array}{cccc}1& 0& 0& 0\\ p_1& 1& 0& 0\\ p_2& p_1& 1& 0\\ p_3& p_2& p_1& 1\end{array}\right]{\left[\begin{array}{cccc}1& 0& 0& 0\\ p_1& 1& 0& 0\\ p_2& p_1& 1& 0\\ p_3& p_2& p_1& 1\end{array}\right]}^T-\left[\begin{array}{cccc}{p}_4& 0& 0& 0\\ p_3& p_4& 0& 0\\ p_2& p_3& p_4& 0\\ p_1& p_2& p_3& p_4\end{array}\right]{\left[\begin{array}{cccc}{p}_4& 0& 0& 0\\ p_3& p_4& 0& 0\\ p_2& p_3& p_4& 0\\ p_1& p_2& p_3& p_4\end{array}\right]}^T$$

$$\text{det}{\Delta }_k(P,P^{\#})>0,k=1,\dots ,4.$$

$$P_k(z)=\frac{1}{z}\left(P_{k-1}(z)-{\nu }_{k-1}{Q}_{k-1}(z)\right),k\ge 1,$$

$$Q_k(z)=Q_{k-1}(z)-{\overline{\nu }}_{k-1}{P}_{k-1}(z),k\ge 1,$$

$${\nu }_k=\frac{P_k(0)}{Q_k(0)},k\ge 0.$$

$$\text{det}{\Delta }_k(P,Q)=(1-|{\nu }_{k-1}^2|)|Q_{k-1}(0){|}^2,k\ge 1,$$

$$|{\nu }_k|<1,k=0,\dots ,3.$$

$${\nu }_0=p_4,{\nu }_1=\frac{p_3-p_4{p}_1}{1-|p_4{|}^2},{\nu }_2=\frac{(1-|p_4{|}^2)(p_2-p_4{p}_2)-(p_3-p_4{p}_1)(p_1-p_4{p}_3)}{(1-|p_4{|}^2{)}^2-|p_3-p_4{p}_1{|}^2},$$

$${\nu }_3=\frac{p_1-p_4{p}_3-\frac{(p_3-p_4{p}_1)(p_2-p_2{p}_4)}{1-|p_4{|}^2}-\frac{(1-|p_4{|}^2)(p_2-p_4{p}_2)-(p_3-p_4{p}_1)(p_1-p_4{p}_3)}{(1-|p_4{|}^2{)}^2-|p_3-p_4{p}_1{|}^2}(p_1-p_4{p}_3-\frac{(p_3-p_4{p}_1)(p_2-p_4{p}_2)}{1-|p_4{|}^2})}{1-|p_4{|}^2-\frac{|p_3-p_4{p}_1{|}^2}{1-|p_4{|}^2}-\frac{(1-|p_4{|}^2)(p_2-p_4{p}_2)-(p_3-p_4{p}_1)(p_1-p_4{p}_3)}{(1-|p_4{|}^2{)}^2-|p_3-p_4{p}_1{|}^2}(p_2-p_4{p}_2-\frac{(p_3-p_4{p}_1)(p_1-p_4{p}_3)}{1-|p_4{|}^2})}$$

$$\left\{\begin{array}[]{l}|p_{4}|<1,\\ |p_{3}-p_{4}p_{1}|<1-|p_{4}|^{2},\\ |(1-|p_{4}|^{2})(p_{2}-p_{4}p_{2})-(p_{3}-p_{4}p_{1})(p_{1}-p_{4}p_{3})|<(1-|p_{4}|^{2})^{2}-|p_{3}-p_{4}p_{1}|^{2},\\ \left|p_{1}-p_{4}p_{3}-\frac{(p_{3}-p_{4}p_{1})(p_{2}-p_{2}p_{4})}{1-|p_{4}|^{2}}-\frac{(1-|p_{4}|^{2})(p_{2}-p_{4}p_{2})-(p_{3}-p_{4}p_{1})(p_{1}-p_{4}p_{3})}{(1-|p_{4}|^{2})^{2}-|p_{3}-p_{4}p_{1}|^{2}}(p_{1}-p_{4}p_{3}-\frac{(p_{3}-p_{4}p_{1})(p_{2}-p_{4}p_{2})}{1-|p_{4}|^{2}})\right|\\ <\left|1-|p_{4}|^{2}-\frac{|p_{3}-p_{4}p_{1}|^{2}}{1-|p_{4}|^{2}}-\frac{(1-|p_{4}|^{2})(p_{2}-p_{4}p_{2})-(p_{3}-p_{4}p_{1})(p_{1}-p_{4}p_{3})}{(1-|p_{4}|^{2})^{2}-|p_{3}-p_{4}p_{1}|^{2}}(p_{2}-p_{4}p_{2}-\frac{(p_{3}-p_{4}p_{1})(p_{1}-p_{4}p_{3})}{1-|p_{4}|^{2}})\right|\end{array}\right.$$

$$(SC)\left\{\begin{array}[]{l}p_{4}<1,\\ |p_{3}-p_{4}p_{1}|<1-(p_{4})^{2},\\ |(1-(p_{4})^{2})(p_{2}-p_{4}p_{2})-(p_{3}-p_{4}p_{1})(p_{1}-p_{4}p_{3})|<(1-(p_{4})^{2})^{2}-|p_{3}-p_{4}p_{1}|^{2},\\ |(1+p_{4})(p_{1}-p_{4}p_{3})-p_{2}(p_{3}-p_{4}p_{1})|<(1-(p_{4})^{2})(1+p_{2}+p_{4})-(p_{1}+p_{3})(p_{3}-p_{4}p_{1})\end{array}\right.$$

$$(SCJ)\left\{\begin{array}[]{l}p_{4}<1,\\ |p_{1}+p_{3}|<1+p_{2}+p_{4},\\ \left|p_{2}(1-p_{4})(1-(p_{4})^{2})-(p_{3}-p_{4}p_{1})(p_{1}-p_{4}p_{3})\right|<(1-(p_{4})^{2})^{2}-(p_{3}-p_{1}p_{4})^{2}.\end{array}\right.$$

$$\left|p_2(1-p_4)+p_4(1-(p_4{)}^2)+p_1(p_4{p}_1-p_3)\right|<p_2{p}_4(1-p_4)+1-(p_4{)}^2+p_3(p_1{p}_4-p_3).$$

$$a=\frac{-{\alpha }_1+{\beta }_{1}x}{1-{\beta }_{0}x},b=\frac{y}{1-{\beta }_{0}x},c=\frac{-{\alpha }_2+{\beta }_{2}x}{1-{\beta }_{0}x},d=\frac{{\gamma }_{2}y}{1-{\beta }_{0}x}$$

$$b^{\ast }=b+\frac{(1+{\eta }_1)xy}{1-{\beta }_{0}x},d^{\ast }=d+\frac{({\gamma }_2+{\eta }_2)xy}{1-{\beta }_{0}x},$$

Full text

AMS subject classification 2010: 60H10, 65C20, 65L20

1. 1 Introduction.

2. 2 Linear Test equation.

3. 3 Stability conditions for two-step Maruyama methods to the scalar test equation.

4. 3.1 Two-step Adams-Bashforth and Adams-Moulton Maruyama methods.

5. 3.2 Schemes for hereditary systems.

6. 4 Linear MS-stability.

7. 4.1 MS-stability of Adams-Bashforth Maruyama scheme.

8. 4.2 MS-stability of Adams-Moulton Maruyama scheme.

9. 4.3 Two-step BDF Maruyama scheme.

10. 5 Experiments.

11. 6 Linear system of SDEs and multi-dimensional noise.

12. 6.1 Stability of two-step methods for linear system of SDEs driven by multi-dimensional noise.

13. 6.2 A linear system of SDEs driven by a single noise term.

14. 6.3 A linear system of SDEs driven by two noise terms.

1. Introduction.

Consider the general type $d$-dimensional Itô stochastic differential equation (sde)

[TABLE]

driven by the $m$-dimensional Wiener process, where the coefficients $F:[0,T]\times{\mathbb{R}}^{d}\mapsto{\mathbb{R}}^{d},G:[0,T]\times{\mathbb{R}}^{d}\mapsto{\mathbb{R}}^{d\times m}$ are such that there exists a unique path-wise strong solution of (1.1), cf. [Mao07, Ch. 2.3]. We will also study complex-valued functions $F,G,X.$ The two-step Maryuama method with an equidistant step-size $h$ for the approximations $X_{i}\approx X(t_{i})$ of the solution of (1.1) read

[TABLE]

and the improved two-step Maryuama method is given by

[TABLE]

for $i=2,3,\ldots,$ in the case of systems with commutative noise; here $\{\xi_{r,i}\}_{i\in{\mathbb{N}}_{0}},r=1,\ldots,m,$ are sequences of i.i.d. standard normal r.v.s, $(\alpha_{j},\beta_{j}),j=0,\ldots,2$ and $(\gamma_{j},\eta_{j}),j=1,2$ are appropriate parameters and $f_{i-j}$ denotes $f(t_{i-j},X_{i-j})$ for appropriate functions $f$ as above. For convergence properties of (1.2) and (1.3) see [BW06], [BW07]. Here, we are interested in mean-square asymptotic properties of the above numerical approximations. We perform a linear stability analysis using linear time-invariant test equations; in [BHBW06] sufficient conditions are given for asymptotic mean-square stability of (1.2) applying appropriate Lyapunov-type functionals. We provide in our main result, Theorem 3, necessary and sufficient conditions for the asymptotic mean-square stability of (1.2) and (1.3) following a different approach.

In Section 2 we use the scalar linear test-equation to study the stability properties of the two-step Maruyama methods. The stability matrix ${\mathcal{S}}$ of the two-step methods is analyzed. Section 3 provides our main result regarding stability conditions for two-step Maruyama methods and applications of it. The linear mean-square stability of the methods is studied in Section 4 and experiments are made in Section 5. Section 6 is devoted to systems of linear test equations with multi-dimensional noise.

2. Linear Test equation.

Consider the scalar linear test-equation with multiplicative noise

[TABLE]

where the coefficients $\lambda,\mu\in{\mathbb{C}}$ and assume w.l.o.g. that $X_{0}$ is non-random. The two-step Maryuama method with an equidistant step-size $h$ for the approximations $X_{i}\approx X(t_{i})$ of the solution of (2.1) read, (apply (1.2) with $F_{i}=\lambda X_{i},G_{i}=\mu X_{i},m=1$)

[TABLE]

and the improved two-step Maryuama method is given by

[TABLE]

where $\{\xi_{i}\}_{i\in{\mathbb{N}}_{0}}$ is a sequence of i.i.d. standard normal r.v.s and $(\alpha_{j},\beta_{j}),j=0,\ldots,2$ and $(\gamma_{j},\eta_{j}),j=1,2$ are appropriate parameters.

The recurrences (2.2) and (2.3) can be rewritten in the form

[TABLE]

where for (2.2) the complex coefficients $a,b,c$ and $d$ read

[TABLE]

and for (2.3) the complex coefficients $a,c$ are the same and $b,d$ read

[TABLE]

In Table 1 we list the coefficients $\alpha_{i},\beta_{i},\gamma_{i},\eta_{i}$ for different two-step Maruyama methods.

The stability or transition matrix ${\mathcal{S}}$ of the two-step method (2.4) reads

[TABLE]

where $\overline{z}$ stands for the conjugate of $z\in{\mathbb{C}}.$ The zero solution of the difference equations (2.4) is asymptotically mean-square stable iff the spectral radius of the mean-square stability matrix ${\mathcal{S}}$ satisfies

[TABLE]

Recall that $\rho({\mathcal{S}}):=\max|l_{j}|$ where $l_{j}$ are the eigenvalues of ${\mathcal{S}}.$ Computing the eigenvalues of ${\mathcal{S}}$ amounts to finding the roots of its characteristic polynomial and verifying condition (2.8). Here the characteristic polynomial is a fourth-order polynomial given by

[TABLE]

where the real coefficients $p_{j},j=1,\ldots,4,$ read

[TABLE]

[TABLE]

We can check $\rho({\mathcal{S}})<1$ avoiding the computation of the $\rho({\mathcal{S}})$ by verifying conditions on the parameters $p_{j},j=1,\ldots,4,$ implied by the Schur-Cohn criterion. The strategy is the following (cf. [Jur88]): Define the transpose $P^{\#}$ of $P$ as

[TABLE]

define the $4\times 4$ Schur-Cohn matrix associated to $P$ by

[TABLE]

where $Q^{T}$ is the transpose of a matrix $Q,$ i.e. the matrix with entries $Q^{T}_{ij}=Q_{ji}$; the polynomial P has all the roots inside the unit disk iff $\Delta_{4}(P,P^{\#})$ is positive definite, which corresponds to

[TABLE]

In order to decide about (2.12) we can use the connection it has with the Schur coefficients $(\nu_{k})_{k=0,\ldots,3}$ of the pair $(P,P^{\#});$ in general for the pair $(P,Q)$ we construct the sequence $(P_{k},Q_{k})_{k=0,1,\ldots,}$ as $P_{0}=P,Q_{0}=Q$ and

[TABLE]

[TABLE]

and take

[TABLE]

Then

[TABLE]

therefore condition (2.12) holds iff the Schur coefficients of the pair $(P,P^{\#})$ satisfy

[TABLE]

In particular the Schur coefficients read

[TABLE]

[TABLE]

and thus (2.13) becomes

[TABLE]

and simplifying the last condition and using $p_{4}>0$ we get

[TABLE]

The Schur-Cohn criterion simplifies to (cf. [Jur91])

[TABLE]

An alternative condition for (SCJ3) reads [Ela05, Ex 5.1, p. 255]

[TABLE]

3. Stability conditions for two-step Maruyama methods to the scalar test equation.

Using the definition of the real coefficients (2.10) and (2.11) and the general conditions (2.17) we can argue when a two-step Maruyama method is asymptotically mean-square stable. In all the following we take $\alpha_{0}=\gamma_{1}=1$ and using (2.5) and (2.6) rewrite the complex coefficients $a,b,c,d$ for the standard schemes

[TABLE]

and $b^{*},d^{*}$ for the improved ones

[TABLE]

where also we have used

[TABLE]

*Theorem 3.1  *** The two-step stochastic linear difference equation (2.4) is asymptotically mean-square stable iff

[TABLE]

[TABLE]

and

[TABLE]

Moreover, if (3.4) holds along with

[TABLE]

and

[TABLE]

then (3.5) is also true. For the improved version we take $b^{*}$ and $d^{*}$ in place of $b$ and $d$ respectively. $\Box$

Proof of Theorem 3.

We rewrite the coefficients $p_{i},i=1,...,4,$ by (2.10) and (2.11)

[TABLE]

[TABLE]

We need to check conditions (2.17) to conclude about the stability of the method. Condition (SCJ1) implies $|c|^{2}(|c|^{2}+|d|^{2})<1.$ Note that

[TABLE]

and

[TABLE]

which give

[TABLE]

We also have

[TABLE]

due to (3.9) and

[TABLE]

so (SCJ2) holds when (3.4) holds. Condition (SCJ3) is (3.5).

Furthermore, (3.4), (3.7) and (3.6) imply

[TABLE]

or

[TABLE]

which combined with (3.10) implies

[TABLE]

Denote the left-hand side of (3.5) by $|L|$ and the right side by $R;$ thus we have to show that $|L|<R.$ Using (SCJ2) and (3.10) we get

[TABLE]

by (3.12). Therefore $L>-R.$ It remains to prove $L<R.$

[TABLE]

by (3.8), (3.12) and (SCJ2). The above is negative if the last term is negative, or equivalently if $p_{2}+3p_{4}+p_{3}-p_{1}<1.$ We have

[TABLE]

where we used (3.6), (3.7) and (3.11). ∎

3.1. Two-step Adams-Bashforth and Adams-Moulton Maruyama methods.

In this case $\gamma_{2}=0$ and thus $d=0$ so the recurrences (2.4) simplify to

[TABLE]

for the standard schemes and to

[TABLE]

for the improved ones.

*Proposition 3.2  *** The two-step stochastic linear difference equation (3.13) is asymptotically mean-square stable iff

[TABLE]

[TABLE]

and

[TABLE]

whereas the two-step stochastic linear difference equation (3.14) is asymptotically mean-square stable iff conditions (3.3),(3.4) and (3.5) hold or conditions (3.6),(3.4) and (3.7) hold where $b^{*}$ and $d^{*}$ are replacing $b$ and $d$ respectively. $\Box$

Proof of Proposition 3.1.

We show the first case since the second one is a direct application of Theorem 3. In the case of (3.13) the coefficients read

[TABLE]

[TABLE]

We apply Theorem 3 when $d=0.$ Conditions (3.3) or (3.6) are equivalent to $|c|<1.$ Condition (3.4) is just the right-side of (3.15) and (3.16). Finally (3.7) shrinks to $\Re(a^{2}\overline{c})\geq-|a|^{2}|c|^{2}.$ ∎

*Remark 3.3  *** Consider the case $a,b,c\in{\mathbb{R}}.$ Then conditions (3.15), (3.16) and (3.17) read (see also [TS14, Cor. 6])

[TABLE]

[TABLE]

$\Box$

3.2. Schemes for hereditary systems.

Hereditary systems are used to model processes in a variety of fields such as physics, biology, economy, just to name a few, (cf. [KM92]). Due to their applications, we present them in a separate subsection. The following stochastic difference equation was proposed in [Sha97],

[TABLE]

where necessary and sufficient conditions were given concerning their asymptotic mean-square stability of the zero solution. By taking the trivial case $l=0$ of this delay system with $k=2$ this falls in our setting (2.4) with $b=0,$ that is,

[TABLE]

for the standard schemes and to

[TABLE]

for the improved ones.

*Proposition 3.4  *** The two-step stochastic linear difference equation (3.23) is asymptotically mean-square stable iff

[TABLE]

and

[TABLE]

whereas the two-step stochastic linear difference equation (3.24) is asymptotically mean-square stable iff conditions (3.3),(3.4) and (3.5) hold or conditions (3.6),(3.4) and (3.7) hold where $b^{*}$ and $d^{*}$ are replacing $b$ and $d$ respectively. $\Box$

Proof of Proposition 3.2.

We show the first case since the second one is a direct application of Theorem 3. In the case of (3.23) the coefficients read

[TABLE]

[TABLE]

We apply Theorem 3 when $b=0.$ Conditions (3.3) or (3.6) are equivalent to $|c|^{2}+|d|^{2}<1.$ Condition (3.4) is just the right-side of (3.25). Finally (3.7) shrinks to $\Re(a^{2}\overline{c})\geq-|a|^{2}|c|^{2}.$ ∎

*Remark 3.5  *** Consider the case $a,c,d\in{\mathbb{R}}.$ Then conditions (3.25) and (3.26) read (see also [TS14, Cor. 5])

[TABLE]

[TABLE]

$\Box$

4. Linear MS-stability.

Recall the scalar linear test-equation (2.1)

[TABLE]

where $\lambda,\mu,X_{0}\in{\mathbb{C}}.$ Its zero solution is asymptotically mean-square stable iff $\Re(\lambda)+|\mu|^{2}/2<0;$ in the case $\mu=0$ the above condition reduces to the notion of A-stability. The set

[TABLE]

is called the mean-square (MS-)stability domain of the stochastic equation (2.1). In an analogous manner the MS-stability domain of a two-step stochastic method (SM) for a given step size $h>0$ is defined as

[TABLE]

In case $\lambda,\mu\in{\mathbb{R}}$ we have the notions of the stability regions

[TABLE]

for the sde and

[TABLE]

for the method. A stochastic method is said to be MS-stable if

[TABLE]

The inverse relation

[TABLE]

means that the method is unstable whenever the test-equation is unstable. In this case the notion of conditional MS-stability comes to play, where one has to determine a step size $h_{0}$ such that for a given pair of $(\lambda,\mu)$ in the stability domain or region of the sde the method is mean-square stable for all $h<h_{0}.$

4.1. MS-stability of Adams-Bashforth Maruyama scheme.

The coefficients of the AB2 scheme, see Table 2, read

[TABLE]

and for the improved AB2I

[TABLE]

First we take $(\lambda,\mu)\in{\mathcal{S}}_{AB2}$ where

[TABLE]

Conditions (3.15) give

[TABLE]

Now, inspecting the second inequality further we conclude that

[TABLE]

when

[TABLE]

which implies $\Re(x)<0,$ that is $\Re(\lambda)<0,$ when $|x|^{2}<\frac{4}{3}.$ Conditions (3.16) give

[TABLE]

when

[TABLE]

which holds for any $0<|x|<2$ with $\Re(x)<0.$ Moreover, condition (3.17) reads,

[TABLE]

or equivalently

[TABLE]

which implies $\Re(x)<0$ when $|x|^{2}<20/9.$ Conditions (3.15), (3.16) and (3.17) hold when

[TABLE]

Therefore we get

[TABLE]

for any $h>0,$ which means that AB2 is unstable whenever the test-equation is unstable. Now, given $(\lambda,\mu)\in{\mathbb{C}}\times{\mathbb{C}}$ we want to find $h_{0}>0$ such that ${\mathcal{S}}_{SDE}\subset{\mathcal{S}}_{AB2}(h)$ for any $h<h_{0}.$ Since we chose the parameters in the stability domain ${\mathcal{S}}_{SDE}$ we have that $|\mu|^{2}<-2\Re(\lambda).$ The relation $|x|<1$ gives

[TABLE]

Moreover, by (3.16) we need to show that

[TABLE]

which holds when

[TABLE]

or in terms of $h$ for

[TABLE]

So given $(\lambda,\mu)\in{\mathcal{S}}_{SDE}$ the method AB2 is conditionally MS-stable for any $h<h_{0}$ where

[TABLE]

In case the parameters are real conditions (3.15), (3.16) and (3.17) shrink to (3.20) and (3.21) respectively by Remark 3.1. The asymptotic region reads

[TABLE]

and ${\mathcal{R}}_{AB2}(h)\subset{\mathcal{R}}_{SDE}$ for any $h>0.$ Given $(\lambda,\mu)\in{\mathcal{R}}_{SDE}$ and $h>0$ the method AB2 is conditionally MS-stable for all $h<h_{0}$ where

[TABLE]

In Figure 1 we represent the stability regions of the AB2 and AB2I scheme respectively in the $(x,Y)$-plane where $x=\lambda h$ and $Y=\mu^{2}h,$ where also the stability region of the SDE is shown (it corresponds to the region $0<Y<-2x,$ that is the light-shaded triangle.)

4.2. MS-stability of Adams-Moulton Maruyama scheme.

The coefficients of the AM2 scheme, see Table 2, read

[TABLE]

and for the improved AM2I

[TABLE]

First we take $(\lambda,\mu)\in{\mathcal{S}}_{AM2}$ where

[TABLE]

Conditions (3.15) give

[TABLE]

The first inequality is satisfied by those $x$ with $5\Re(x)<|x|^{2}+6$ and the second inequality holds when $-6<\Re(x)<0.$ Therefore (4.7) holds iff

[TABLE]

which imply $\Re(\lambda)<0.$ Conditions (3.16) give

[TABLE]

which is smaller than $-2\Re(x)$ for any $\Re(x)<0.$ Moreover, condition (3.17) reads,

[TABLE]

or equivalently

[TABLE]

which implies $\Re(x)<0$ and $|x|^{2}<9/2.$ Therefore we get

[TABLE]

for any $h>0,$ which means that AM2 is unstable whenever the test-equation is unstable. Now, given $(\lambda,\mu)\in{\mathbb{C}}\times{\mathbb{C}}$ we want to find $h_{0}>0$ such that ${\mathcal{S}}_{SDE}\subset{\mathcal{S}}_{AM2}(h)$ for any $h<h_{0}.$ Since we chose the parameters in the stability domain ${\mathcal{S}}_{SDE}$ we have that $|\mu|^{2}<-2\Re(\lambda).$ Relation (4.8) implies

[TABLE]

Moreover, by (3.16) we need to show that

[TABLE]

which holds for sufficiently small $h_{1}>0$ implying that the method AM2 is conditionally MS-stable for any $h<h_{0}$ where

[TABLE]

In case the parameters are real we need to show (3.20) and (3.21) respectively. The asymptotic region reads

[TABLE]

and ${\mathcal{R}}_{AM2}(h)\subset{\mathcal{R}}_{SDE}$ for any $h>0.$ Given $(\lambda,\mu)\in{\mathcal{R}}_{SDE}$ and $h>0$ the method AM2 is conditionally MS-stable for all $h<h_{0}$ where

[TABLE]

In Figure 2 we represent the stability regions of the AM2 and AM2I scheme respectively in the $(x,Y)$-plane where $x=\lambda h$ and $Y=\mu^{2}h,$ where also the stability region of the SDE is shown (it corresponds to the region $0<Y<-2x,$ that is the light-shaded triangle.)

4.3. Two-step BDF Maruyama scheme.

The coefficients of the BDF2 scheme, see Table 2, read

[TABLE]

and for the improved BDF2I

[TABLE]

In Figure 3 we represent the stability regions of the BDF2 and BDF2I scheme respectively in the $(x,Y)$-plane where $x=\lambda h$ and $Y=\mu^{2}h,$ where also the stability region of the SDE is shown (it corresponds to the region $0<Y<-2x,$ that is the light-shaded triangle.)

5. Experiments.

In this section we make some simple numerical experiments to complement the stability analysis presented in the previous section. We apply the AB2, AM2 and BDF2 two-step Maruyama schemes as well as their improved versions with constant step-size $h,$ to solve the equation

[TABLE]

For the second initial condition in the two-step schemes we apply the $\theta-$Maruyama method which applied to the linear test equation (2.1) reads

[TABLE]

with $\theta=1/2$ and $n=0.$ We also implement the $\theta-$method (with $\theta=1/2$) and the Euler method (with $\theta=0$) for further comparison. We plot the obtained values in a $log_{2}$-scale against time $t.$ The estimated mean-square norm of $X$ is point-wise estimated by each stochastic numerical method $X^{SM}$ in the following way,

[TABLE]

where we have computed $M$ batches of $L$ simulation paths. The total number of paths in the experiments is $M\cdot L=10^{6}.$ For the first experiment, see Figure 4, we have applied all the methods with time-step size $h=1/8,$ so that they are all stable. The considered time interval is $[0,1].$ In the second experiment, see Figure 5, we integrate over $[0,20]$ with $h=1.$ In this case the $AB2,AB2I$ and $AM2,AM2I$ methods are not stable as well as the forward Euler method. The $BDF2,BDF2I$ methods as well as the $\theta-$methods are asymptotically stable in the mean-square sense with the $BDF2$ performing better. Another remark we can make in the one-dimensional case is about the performance of the proposed improved methods with respect to their stability behavior, which seems to follow the rule that we do not gain more w.r.t to stability performance by using higher order schemes (multiple integrals for the approximation of the diffusion coefficient), as one can see from both Figures 4 and 5. Nevertheless, the situation is different in more dimensions as shown n Section 6.

6. Linear system of SDEs and multi-dimensional noise.

Consider the $d$-system of linear test-equations with $m$-dimensional multiplicative noise

[TABLE]

where $F,G$ are $d\times d$ real-valued matrices and assume w.l.o.g. that $X_{0}$ is non-random.

6.1. Stability of two-step methods for linear system of SDEs driven by multi-dimensional noise.

The two-step Maryuama method with an equidistant step-size $h$ and approximations $X_{i}=(X_{1,i},X_{2,i},\ldots,X_{n,i})^{T}$ of the solution of (6.1) read

[TABLE]

and can be represented as

[TABLE]

where the matrices $A,C$ and $B_{r},D_{r}$ are given by

[TABLE]

[TABLE]

and for the improved versions

[TABLE]

for $r=1,\ldots,m.$

Here, the stability or transition matrix ${\mathcal{S}}$ of the two-step method (6.3) applied to linear system of the form (6.1) reads

[TABLE]

with $R=\sum_{r=1}^{m}(A\otimes D_{r})(B_{r}\otimes{\mathbb{I}}_{d})+\sum_{r=1}^{m}(D_{r}\otimes A)({\mathbb{I}}_{d}\otimes B_{r}).$

A result of the type of Theorem 3, that is a conclusion about the asymptotically zero mean-square stability of the two-step method (6.2) applied to the linear system (6.1), is again related with equivalent conditions for the relation $\rho({\mathcal{S}})<1.$ Now the characteristic polynomial of the stability matrix ${\mathcal{S}}$ is of order $4d^{2}.$ The computational effort of the Schur-Cohn test (SCJ) is now bigger, but one can reduce it by halving the dimensions of the matrix, whose positive-definite character needs to be checked at the expense of some easily checked inequalities on linear combinations of the coefficients of the polynomial (c.f. [AJ73]).

6.2. A linear system of SDEs driven by a single noise term.

Consider the system of linear test-equations (6.1) with $d=2,m=1$ and matrices $F,G$ of the following type

[TABLE]

that is

[TABLE]

with a single noise term. The mean-square stability matrix for (6.9) is

[TABLE]

and the zero solution of (6.9) is asymptotically MS-stable iff (cf. [BS12, Lemma 4.1])

[TABLE]

Below we make a simple experiment implementing the two-step Maruyama methods

[TABLE]

where in particular for the AB2/AB2I methods

[TABLE]

[TABLE]

for the AM2/AM2I methods

[TABLE]

[TABLE]

with $Q=({\mathbb{I}}_{2}-\frac{5}{12}hF)^{-1}$ and for the BDF2/BDF2I methods

[TABLE]

[TABLE]

with $Q=({\mathbb{I}}_{2}-\frac{2}{3}hF)^{-1}.$ We choose the values of $\lambda,\sigma,\epsilon$ such that the spectral abscissa $\alpha({\mathcal{S}})$ of the mean-square stability matrix ${\mathcal{S}}$ is negative, that is $\alpha({\mathcal{S}})<0$, and the spectral radius $\rho({\mathcal{S}})<1$ or in other words such that the condition (6.11) holds. In this case, see Figure 6, the improved versions AM2I and BDF2I are stable whereas AM2 and BDF2 are not.

6.3. A linear system of SDEs driven by two noise terms.

Consider the system of linear test-equations (6.1) with $d=2,m=2$ and matrices $F,G_{1},G_{2}$ of the following type

[TABLE]

that is

[TABLE]

with two commutative noise terms. The mean-square stability matrix for (6.20) is

[TABLE]

and the zero solution of (6.20) is asymptotically MS-stable iff (cf. [BS12, Lemma 4.1])

[TABLE]

Below we make a simple experiment implementing the two-step Maruyama methods

[TABLE]

where for all methods $A$ and $B$ are as in (6.23) and $B_{r},D_{r}$ correspond now to the matrices $G_{r};$ for instance for the AB2/AB2I methods we have

[TABLE]

We choose the values of $\lambda,\sigma,\epsilon$ in a way that the condition (6.22) holds and compute the MS-norm of $X^{(1)},$ just as in [BS12], by

[TABLE]

In this case, see Figure 7, the improved versions AB2I, AM2I and BDF2I are stable whereas AB2, AM2 and BDF2 are not.

Of course, by lowering the step-size the numerical methods become more stable. In the following, we sequentially halve the step-size and confirm the conjecture above. In all cases though, we conclude again that AB2, AM2 and BDF2 are less stable than their improved counterparts, see Figures 8,9 and for a clearer view Figures 10,11 and 12.

Acknowledgments

The author would like to thank Prof. Evelyn Buckwar for fruitful discussions around the subject. It took place during the author’s visit at the Institute for Stochastics, in Linz, Austria in the last four months of 2016.