Mean-square stability of numerical schemes for stochastic differential systems

Yoshihiro SAITO (Gifu Shotoku Gakuen University) and Taketomo MITSUI (Nagoya University)

Abstract

Stochastic differential equations (SDEs) represent physical phenomena dominated by stochastic processes. Similar to deterministic ordinary differential equations (ODEs), various numerical schemes are proposed for SDEs. Stability analysis is significant for numerical SDEs as well, however a few results have been known. We have proposed the mean-square stability of numerical schemes for a scalar SDE, that is, the numerical stability with respect to the mean-square norm. We studied it, however, only for scalar SDEs because of difficulty and complexity in SDE systems. Trying to make a breakthrough, in the present note we will consider a 2-dimensional linear system with one multiplicative noise and give stability criteria under several notions of the matrix norm.

1 Introduction

We ([7]) proposed the numerical mean-square stability (MS-stability) for a scalar stochastic differential equation (SDE) with one multiplicative noise. However we studied it for only scalar SDEs. Komori and Mitsui [4, 5] analyzed numerical MS-stability for a 2-dimensional SDE in a special case, that is, simultaneously diagonalizable case. In this note we will try to analyze numerical MS-stability of the Euler-Maruyama scheme for more general 2dimensional SDE systems.

Consider the SDE of Ito-type given by

dX(t) = f(t, X) dt + g(t, X) dW (t)  (1)

with f(0, t) = g(0, t) = 0 so that the steady state X(t) = 0 is the equilibrium solution. The Euler-Maruyama scheme for the discrete approximate solution \{\overline{X}_{n}\} over the step-points {t_{n}} is given by

\overline{X}_{n+1} = \overline{X}_{n} + f(t_{n}, \overline{X}_{n})h + g(t_{n}, \overline{X}_{n})ΔW_{n}

where h and ΔW_{n} stand for the step-size and the increment of the Wiener process, respectively. Then we can give the definition of the MS-stability.

Definition 1 Steady solution X(t) ≡ 0 is asymptotically stable in mean-square if the estimations

∀ε > 0, ∃δ > 0; E (∥X(t)∥^{2}) < ε for all t ≥ 0 and ∥X_{0}∥ < δ

and

∃δ_{0};  lim_{t→∞} E (∥X(t)∥^{2}) = 0 for all ∥X_{0}∥ < δ_{0}

hold.

Here the norm ∥x∥ stands for the Euclidean norm of a vector x ∈ ℝ^{2}.

The concept of numerical stability means whether a numerical solution can keep a similar asymptotic property as n tends to infinity when it is applied to the asymptotically stable SDE in mean-square. We will consider a general type of linear SDE systems, for we take the standpoint of linear stability analysis. In the next section we describe criteria of MSstability for the SDE system. Section 3 shows the conditions of numerical MS-stability of the Euler-Maruyama scheme corresponding to Section 2. In Section 4 we will show the numerical experiments confirming our stability analysis in Section 3. Finally we will give our conclusions and future aspects.

2 Criteria of MS-stability

To carry out a linear stability analysis, we will restrict the SDE (1) to an Ito-type 2dimensional linear SDE system with one multiplicative noise, which has the form

\{ \begin{matrix}dX(t) & = & DX(t) dt + BX(t) dW (t), \\ X(0) & = & 1. \end{matrix}  (2)

Here the real constant matrices D and B are given by

D = [ \begin{matrix}λ_{1} & 0 \\ 0 & λ_{2} \end{matrix} ] and B = [ \begin{matrix}α_{1} & β_{1} \\ β_{2} & α_{2} \end{matrix} ] .

Komori and Mitsui [4, 5] analyzed MS-stability for SDE system (2) with β_{1} = 0 and β_{2} = 0, that is, for the simultaneously diagonalizable system. We will consider more general SDE system, namely β_{1} ≠ 0 and β_{2} ≠ 0. First we will introduce the conventional and the logarithmic norms of matrices for the analysis.

Definition 2 Corresponding to the vector norms l^{1}, l^{2} and l^{∞} in ℝ^{n}, we define the subordinate matrix norms of square n × n matrix A = (a_{ij}) by

∥A∥_{1} = max_{j} {∑_{i=1}^{n} |a_{ij}|} , ∥A∥_{∞} = max_{i} {∑_{j=1}^{n} |a_{ij}|} ,
∥A∥_{2} = {maximum eigenvalue of A^{T} A}^{1/2} .

The following can be found in [1, 6].

Definition 3 Logarithmic matrix norm μ_{p}[A] (p = 1, 2, ∞) is defined by

μ_{p}[A] = lim_{δ→0+}(∥I + δA∥_{p} − 1)/δ

where I is the unit matrix and h ∈ ℝ.

The following identities are well known to evaluate the logarithmic norms.

μ_{1}[A] = max_{j} {a_{jj} + ∑_{i≠j} |a_{ij}|} , μ_{∞}[A] = max_{i} {a_{ii} + ∑_{j≠i} |a_{ij}|} ,
μ_{2}[A] = maximum eigenvalue of (A + A^{T})/2.

Let P (t) = E(X(t)X(t)^{T}) be the 2 × 2 matrix-valued second moment of the solution of (2). Then P (t) obeys the initial value problem of the following matrix ordinary differential equation (ODE)

\frac{dP}{dt} = DP + P D^{T} + BP B^{T} (t > 0),  (3)

with P (0) = X_{0}X_{0}^{T}. By virtue of the symmetry of the matrix P we have its governing ODEs of 3-dimension

\frac{dY}{dt} = MY  (4)

where

Y (t) = (Y^{1}(t), Y^{2}(t), Y^{3}(t)), Y^{1}(t) = E(X^{1}(t))^{2},
Y^{2}(t) = E(X^{2}(t))^{2}, Y^{3}(t) = E(X^{1}(t)X^{2}(t)).

We can readily obtain the following lemma owing to the logarithmic matrix norm μ_{p}.

Lemma 1 The linear test system with the unit initial value is asymptotically MS-stable w.r.t. logarithmic norm μ_{p} iff

μ_{p}[M] < 0

We will study the MS-stability w.r.t. μ_{∞} for the test system (2). In fact, a direct calculation brings the expression of the matrix in (4) as

M = [ \begin{matrix}2λ_{1} + α_{1}^{2} & β_{1}^{2} & 2α_{1}β_{1} \\ β_{2}^{2} & 2λ_{2} + α_{2}^{2} & 2α_{2}β_{2} \\ α_{1}β_{2} & α_{2}β_{1} & λ_{1} + λ_{2} + α_{1}α_{2} + β_{1}β_{2} \end{matrix} ] .  (5)

Then we reach the following criterion.

Theorem 1 The system (2) is MS-stable w.r.t. μ_{∞} if the estimation

max {2λ_{1} + (|α_{1}| + |β_{1}|)^{2}, 2λ_{2} + (|α_{2}| + |β_{2}|)^{2}} < 0  (6)

holds.

Proof. Eq. (5) implies that μ_{∞}[M] equals the maximum of the following three quantities.

2λ_{1} + α_{1}^{2} + β_{1}^{2} + 2|α_{1}β_{1}| = 2λ_{1} + (|α_{1}| + |β_{1}|)^{2},
2λ_{2} + α_{2}^{2} + β_{2}^{2} + 2|α_{2}β_{2}| = 2λ_{2} + (|α_{2}| + |β_{2}|)^{2},
λ_{1} + λ_{2} + α_{1}α_{2} + β_{1}β_{2} + |α_{1}β_{2}| + |α_{2}β_{1}|

We see the following inequality straightforwardly.

λ_{1} + λ_{2} + α_{1}α_{2} + β_{1}β_{2} + |α_{1}β_{2}| + |α_{2}β_{1}|
≤ λ_{1} + λ_{2} + |α_{1}α_{2}| + |β_{1}β_{2}| + |α_{1}β_{2}| + |α_{2}β_{1}|
= λ_{1} + λ_{2} + (|α_{1}| + |β_{1}|)(|α_{2}| + |β_{2}|)
≤ \frac{2λ_{1} + (|α_{1}| + |β_{1}|)^{2}}{2} + \frac{2λ_{2} + (|α_{2}| + |β_{2}|)^{2}}{2}

Thus we have (6).

Hereafter we assume λ_{1} < λ_{2} < 0, for it is a natural condition of asymptotic stability of ODEs. Theorem 1 has several consequencies as its corollaries.

Corollary 1 The singly anti-diagonal case (SAD case) of the diffusion matrix implies

B = [ \begin{matrix}0 & β \\ β & 0 \end{matrix} ] ,

which yields

M = [ \begin{matrix}2λ_{1} & β^{2} & 0 \\ β^{2} & 2λ_{2} & 0 \\ 0 & 0 & λ_{1} + λ_{2} + β^{2} \end{matrix} ] .

Then we have the stability criterion w.r.t. μ_{∞} by

max{2λ_{1} + β^{2}, 2λ_{2} + β^{2}} < 0.  (7)

Remark. In SAD case, since the matrix M has a special form, the condition (7) turns out to be that w.r.t. μ_{1} as well.

Note that the condition represented by μ_{∞} is a sufficient condition for the convergence to the zero solution. We will show this through the following example of SAD case.

Example 1 The combination with

D = [ \begin{matrix}−100 & 0 \\ 0 & −1 \end{matrix} ] and B = [ \begin{matrix}0 & 2 \\ 2 & 0 \end{matrix} ]

yields

M = [ \begin{matrix}−200 & 4 & 0 \\ 4 & −2 & 0 \\ 0 & 0 & −97 \end{matrix} ] ,

whose logarithmic norms are

μ_{∞}(M) = 2 > 0 but μ_{2}(M) = −101 + \sqrt{9817} < 0.

Corollary 2 The singly diagonal and anti-diagonal case (SDAD case) of the diffusion matrix implies

B = [ \begin{matrix}α & β \\ β & α \end{matrix} ] ,

which yields

M = [ \begin{matrix}2λ_{1} + α^{2} & β^{2} & 2αβ \\ β^{2} & 2λ_{2} + α^{2} & 2αβ \\ αβ & αβ & λ_{1} + λ_{2} + α^{2} + β^{2} \end{matrix} ] .

Therefore the SDAD case brings the stability criterion w.r.t. μ_{∞} as

max {2λ_{1} + (|α| + |β|)^{2}, 2λ_{2} + (|α| + |β|)^{2}} < 0.  (8)

3 MS-stability of Euler-Maruyama scheme

We now ask what conditions must be imposed in order that the numerical solution \{\overline{X}_{n}\} of (2) generated by a numerical scheme satisfies

\overline{Y}_{n} = E|\overline{X}_{n}|^{2} → 0 as n → ∞.  (9)

In [7], we showed the (scalar) numerical stability factor |R(\overline{h}, k)| of the Euler-Maruyama scheme as

|R(\overline{h}, k)| = |1 + \overline{h}^{2}| + |k\overline{h}|

with \overline{h} = hλ and k = μ^{2}/λ when the scheme is applied to the scalar test equation

dX(t) = λX dt + μX dW (t) (λ, μ ∈ ℂ).

The region R_{EM} defined by

R_{EM} = \{(h̄, k); |R(\overline{h}, k)| < 1 holds\}

is called the MS-stability region of the Euler-Maruyama scheme in the scalar case. The region is displayed in Fig. 1 in the case of λ, μ ∈ ℝ. Through our analysis we will show that the stability factor R(\overline{h}, k) as well as the MS-stability region R_{EM} is still efficient for the linear system (2).

When we apply a numerical scheme to (2) and calculate the components of the second moment of \overline{X}_{n}, we obtain a one-step difference equation of the form

\overline{Y}_{n+1} = \overline{M} \overline{Y}_{n}  (10)

Figure 1: MS-stability region of the Euler-Maruyama scheme

where

\overline{Y}_{n} = (\overline{Y} _{n}^{1}, \overline{Y} _{n}^{2}, \overline{Y} _{n}^{3}), \overline{Y} _{n}^{1} = E(\overline{X}_{n}^{1})^{2}, \overline{Y} _{n}^{2} = E(\overline{X}_{n}^{2})^{2}, \overline{Y} _{n}^{3} = E(\overline{X}_{n}^{1}\overline{X}_{n}^{2}).  (11)

We shall call \overline{M} the stability matrix of the scheme. Note that \overline{Y}_{n} → 0 as n → ∞ if

∥\overline{M}∥_{p} < 1.  (12)

Therefore we introduce

Definition 4 The numerical scheme is said to be MS-stable w.r.t. ∥·∥_{p} if it has \overline{M} satisfying ∥\overline{M}∥_{p} < 1.

We will calculate the stability matrix \overline{M} and an MS-stability criterion w.r.t. ∥ · ∥_{∞} of the Euler-Maruyama scheme for the system (2). In the following, let the symbol r(x) stand for 1 + x.

Theorem 2 For the system (2) we obtain

\overline{M} = [ \begin{matrix}r^{2}(λ_{1}h) + α_{1}^{2}h & β_{1}^{2}h & 2α_{1}β_{1}h \\ β_{2}^{2}h & r^{2}(λ_{2}h) + α_{2}^{2}h & 2α_{2}β_{2}h \\ α_{1}β_{2}h & α_{2}β_{1}h & r(λ_{1}h)r(λ_{2}h) + (α_{1}α_{2} + β_{1}β_{2})h \end{matrix} ] ,

which yields the numerical stability criterion w.r.t. ∥ · ∥_{∞} as

max{(1 + λ_{1}h)^{2} + (|α_{1}| + |β_{1}|)^{2}h, (1 + λ_{2}h)^{2} + (|α_{2}| + |β_{2}|)^{2}h} < 1.  (13)

Proof runs similarly as that for Theorem 1. That is, the inequality

|r(λ_{1}h)r(λ_{2}h) + (α_{1}α_{2} + β_{1}β_{2})h| + |α_{1}β_{2}|h + |α_{2}β_{1}|h
≤ |r(λ_{1}h)r(λ_{2}h)| + |α_{1}α_{2}|h + |β_{1}β_{2}|h + |α_{1}β_{2}|h + |α_{2}β_{1}|h
= |r(λ_{1}h)r(λ_{2}h)| + (|α_{1}| + |β_{1}|)(|α_{2}| + |β_{2}|)h
≤ \frac{|r(λ_{1}h)|^{2} + |r(λ_{2}h)|^{2}}{2} + \frac{(|α_{1}| + |β_{1}|)^{2} + (|α_{2}| + |β_{2}|)^{2}}{2} h
= \frac{|r(λ_{1}h)|^{2} + (|α_{1}| + |β_{1}|)^{2}h}{2} + \frac{|r(λ_{2}h)|^{2} + (|α_{2}| + |β_{2}|)^{2}h}{2}

is employed to derive the above result.

Observing the left-hand side of the MS-stability condition (13), we conclude that we can check the numerical MS-stability whether the pair (\overline{h}, k) ≡ (λh, (|α|^{2} + |β|^{2})/λ) satisfies |R(\overline{h}, k)| < 1 for both (λ_{1}, α_{1}, β_{1}) and (λ_{2}, α_{2}, β_{2}). Namely, we should check whether the inclusions

(\overline{h}_{1}, k_{1}) = (λ_{1}h, (|α_{1}|^{2} + |β_{1}|^{2})/λ_{1}), (\overline{h}_{2}, k_{2}) = (λ_{2}h, (|α_{2}|^{2} + |β_{2}|^{2})/λ_{2}) ∈ R_{EM}

hold for the Euler-Maruyama scheme applied to the system (2).

4 Numerical experiments

We will confirm our MS-stability analysis of the Euler-Maruyama scheme through numerical experiments. We will describe four examples as follows. All the examples have the same initial condition as X(0) = (1, 1). In each calculation, 10, 000 sample paths are generated and the values \overline{Y} _{n}^{1}, \overline{Y} _{n}^{2} and \overline{Y} _{n}^{3} in (11) are plotted versus the time-axis in the figures.

Example 2 Simultaneously diagonalizable case described in [4, 5].

dX = [ \begin{matrix}−200 & 0 \\ 0 & −100 \end{matrix} ] X dt + [ \begin{matrix}10 & 0 \\ 0 & 10 \end{matrix} ] X dW (t)  (14)

Since λ_{1} = −200, λ_{2} = −100, α_{1} = α_{2} = 10 and β_{1} = β_{2} = 0, we can discriminate the following cases. h = 0.005, (\overline{h}, k) = (−1, −0.5), (−0.5, −1) : stable h = 0.01, (\overline{h}, k) = (−2, −0.5), (−1, −1) : unstable h = 0.02, (\overline{h}, k) = (−4, −0.5), (−2, −1) : unstable h = 0.05, (\overline{h}, k) = (−10, −0.5), (−5, −1) : unstable

Numerical results are in Fig. 2.

Example 3 SAD case given in Cor. 1.

dX = [ \begin{matrix}−200 & 0 \\ 0 & −100 \end{matrix} ] X dt + [ \begin{matrix}0 & 10 \\ 10 & 0 \end{matrix} ] X dW (t)

Figure 2: Example 2. Upper left: h = 0.005, upper right: h = 0.01, lower left: h = 0.02 and lower right: h = 0.05.

Since λ_{1} = −200, λ_{2} = −100, α = 0 and β = 10, we can have the following cases similar to Example 2. h = 0.005, (\overline{h}, k) = (−1, −0.5), (−0.5, −1) : stable h = 0.01, (\overline{h}, k) = (−2, −0.5), (−1, −1) : unstable h = 0.02, (\overline{h}, k) = (−4, −0.5), (−2, −1) : unstable h = 0.05, (\overline{h}, k) = (−10, −0.5), (−5, −1) : unstable

Numerical results are in Fig. 3.

Example 4 SDAD case given in Cor. 2.

dX = [ \begin{matrix}−200 & 0 \\ 0 & −100 \end{matrix} ] X dt + [ \begin{matrix}10 & 2 \\ 2 & 10 \end{matrix} ] X dW (t)

Since λ_{1} = −200, λ = −100, α = 10 and β = 2 in (8), we can discriminate the following two cases. h = 0.005, (\overline{h}, k) = (−1, −0.72), (−0.5, −1.44) : stable h = 0.01, (\overline{h}, k) = (−2, −0.72), (−1, −1.44) : unstable

Numerical results are in Fig. 4.

Figure 3: Example 3. Upper left: h = 0.005, upper right: h = 0.01, lower left: h = 0.02 and lower right: h = 0.05.

Example 5 Another SDAD case.

dX = [ \begin{matrix}−200 & 0 \\ 0 & −100 \end{matrix} ] X dt + [ \begin{matrix}2 & 10 \\ 10 & 2 \end{matrix} ] X dW (t)

Since λ_{1} = −200, λ = −100, α = 2 and β = 10 in (8), discrimination in the numerical stability turns out to be same as in Example 4. h = 0.005, (\overline{h}, k) = (−1, −0.72), (−0.5, −1.44) : stable h = 0.01, (\overline{h}, k) = (−2, −0.72), (−1, −1.44) : unstable

Numerical results are in Fig. 5.

Figure 4: Example 4. Left: h = 0.005 and right: h = 0.01.

Example 6 Another SDAD case.

dX = [ \begin{matrix}−200 & 0 \\ 0 & −100 \end{matrix} ] X dt + [ \begin{matrix}5 & 10 \\ 10 & 5 \end{matrix} ] X dW (t)

Since λ_{1} = −200, λ = −100, α = 5 and β = 10 in (8), discrimination in the numerical stability turns out to be same as in Example 4. h = 0.005, (\overline{h}, k) = (−1, −0.625), (−0.5, −1.25) : stable h = 0.01, (\overline{h}, k) = (−2, −0.625), (−1, −1.25) : unstable

Numerical results are in Fig. 5.

5 Conclusions and Future aspects

In this note we have extended the numerical MS-stability analysis from a scalar SDE with one multiplicative noise into a 2-dimensional SDE system with one multiplicative noise for the Euler-Maruyama scheme, and showed that the MS-stability region for the scalar case is still efficient for the system case. This must suggest a way for linear stability analysis of the

Figure 5: Example 5. Left: h = 0.005 and right: h = 0.01.

SDE system in the MS sense. Therefore we will analyze MS-stability for the matrices D and B with complex numbers and of more dimension. And we will investigate the relationship of the MS-stability conditions in different matrix norms, for example, between ∥ · ∥_{∞} and ∥ · ∥_{2}. Also we plan to MS-stability analysis of other numerical schemes.

Acknowledgements

This work is supported by the Grant-in-Aid for Scientific Research (A) of Japan Society for Promotion of Science (JSPS) (No. 11304004).

References

[1] K. Dekker and J. G. Verwer, Stability of Runge-Kutta methods for stiff nonlinear differential equations, North-Holland, Amsterdam, 1984.

[2] T. C. Gard, Introduction to Stochastic Differential Equations, Marcel Dekker, New York, 1988.

[3] P. E. Kloeden and E. Platen, The Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 1992.

[4] Y. Komori and T. Mitsui, Stable ROW-type weak scheme for stochastic differential equations, Monte Carlo Methods and Applic., 1(1995), 279-300.

[5] Y. Komori, Y. Saito and T. Mitsui, Some issues in discrete approximate solution for stochastic differential equations, Computers Math. Applc., 28 (1994), no.10-12, 269– 278.

[6] J. D. Lambert, Numerical Methods for Ordinary Differential Systems, Wiley, New York, 1991.

[7] Y. Saito and T. Mitsui, Stability analysis of numerical schemes for stochastic differential equations, SIAM J. Numer. Anal., 33 (1996), 2254–2267.