Mean-square stability of Milstein method for linear hybrid stochastic delay integro-differential equations

Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

Contents lists available at ScienceDirect

Nonlinear Analysis: Hybrid Systems journal homepage: www.elsevier.com/locate/nahs

Mean-square stability of Milstein method for linear hybrid stochastic delay integro-differential equations A. Rathinasamy a,∗ , K. Balachandran b a

Department of Mathematics and Computer Applications, PSG College of Technology, Coimbatore 641004, India

b

Department of Mathematics, Bharathiar University, Coimbatore 641046, India

article

info

Article history: Received 5 May 2008 Accepted 19 September 2008 Keywords: Stochastic delay integro-differential equations Mean-square stability Numerical methods Markovian switching

a b s t r a c t In this paper we study the mean-square (MS) stability of the Milstein method for linear stochastic delay integro-differential equations (SDIDE) with Markovian switching by extending the techniques of [Z. Wang, C. Zhang, An analysis of stability of Milstein method for stochastic differential equations with delay, Computers and Mathematics with Applications 51 (2006) 1445–1452; L. Ronghua, H. Yingmin, Convergence and stability of numerical solutions to SDDEs with Markovian switching, Applied Mathematics and Computation 175 (2006) 1080–1091]. It is established that the Milstein method is MSstable for linear stochastic delay differential equations (Wang and Zhang (2006); in the above reference). Here we prove that it is MS-stable for linear SDIDE with Markovian switching also under suitable conditions on the integral term. A numerical example is provided to illustrate the theoretical results. © 2008 Elsevier Ltd. All rights reserved.

1. Introduction The hybrid systems driven by continuous-time Markov chains have been used to model many practical systems [e.g. hybrid queueing systems [1]] where they may experience abrupt changes in their structure and parameters caused by phenomena such as component failures or repairs, changing subsystem interconnections, and abrupt environmental disturbances. The hybrid systems combine a part of the state that takes values continuously and another part of the state that takes discrete values. One of the important classes of hybrid systems is the stochastic delay differential equations (SDDEs) with Markovian switching. The natural extension of this system is stochastic delay integro-differential equations (SDIDE) with Markovian switching. Explicit solutions can hardly be obtained for the SDIDEs with Markovian switching. Thus, it is necessary to develop appropriate numerical methods and to study the properties of these approximate schemes. Stability of numerical schemes for SDDEs is essential to avoid a possible explosion of numerical solutions. The convergence and stability properties of the numerical methods for the stochastic ordinary differential equations have been studied by many authors, for instance, Higham [2], Kloeden and Platen [3], and Milstein [4]. The stability of the Euler–Maruyama numerical solutions for stochastic differential delay equations is studied by Cao et al. [5]. The convergence and stability of the semi-implicit Euler method for a linear stochastic differential delay equation is discussed by Liu et al. [6]. The stability of the Milstein method for stochastic differential delay equations is studied by Wang and Zhang [7]. The theoretical study of stochastic differential delay equations with Markovian switching is considered by several authors [8– 10]. Rathinasamy and Balachandran [11] analyzed the numerical MS-stability of second-order Runge–Kutta schemes for multi-dimensional stochastic differential systems with one multiplicative noise. The convergence and stability of the Euler

∗

Corresponding author. E-mail addresses: [email protected] (A. Rathinasamy), [email protected] (K. Balachandran).

1751-570X/$ – see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.nahs.2008.09.015

A. Rathinasamy, K. Balachandran / Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

1257

method for a class of linear stochastic differential delay equations with Markovian switching is studied by Ronghua and Yingmin [12]. Koto [13] investigated the stability of the numerical methods for the integro-differential equation whereas Shaikhet and Rorberts [14] studied the stability of the numerical methods for the linear stochastic functional differential equations. Ding et al. [15] studied the convergence and stability of the semi-implicit Euler method for linear stochastic delay integro-differential equations. In this paper we consider the linear n-dimensional stochastic delay integro-differential equation with Markovian switching of the form



dx(t ) = A(r (t ))x(t ) + B(r (t ))x(t − τ ) + C (r (t ))

+ [D(r (t ))x(t ) + E (r (t ))x(t − τ )] dw(t ), x(t ) = ψ(t ), t ∈ [−τ , 0],

t

Z

t −τ





 x(s)ds dt 

t ≥ 0,

(1)

 

where r (t ) is a Markov chain taking values in S = {1, 2, . . . , N } and A(.), B(.), C (.), D(.) and E (.) ∈ R, {w(t )} is a standard one-dimensional Brownian motion and τ is a positive fixed delay. In Section 2 we will introduce some necessary notations and assumptions. Further we shall discuss the exponential stability in mean squares for Eq. (1) in the same section. In Section 3, the Milstein method will be used to produce the numerical solutions. Furthermore, our main result will be shown and proved in Section 4 and a numerical example is given in Section 5. 2. Notation and definitions Let (Ω , F , {Ft }t ≥0 , P ) be a complete probability space with a filtration (Ft )t ≥0 , which satisfies the usual conditions (i.e., it is increasing and right continuous while F0 contains all P − nullsets). Let w(t ), t ≥ 0 in Eq. (1) be Ft -adapted and independent of F0 . Moreover, |.| is the Euclidean norm in Rn and the matrix norm of matrix G is defined by kGk = sup {|Gx| : |x| = 1}. We assume the initial segment ψ(t ), t ∈ [−τ , 0] to be F0 -measurable and right continuous, E kψk2 < ∞, where kψk = sup−τ ≤t ≤0 |ψ(t )|. Let r (t ), t ≥ 0, be a right-continuous
Markov chain on the probability space taking values in a finite state space S = {1, 2, . . . , N } with generator Γ = γij N ×N given by P {r (t + ∆) = j|r (t ) = i} =



γij ∆ + o (∆) 1 + γii ∆ + o (∆)

if i 6= j, if i = j,

where ∆ > 0. Here γij ≥ 0 is the transition rate from i to j if i 6= j while γii = − j6=i γij . We assume that the Markov chain r (.) is independent of the Brownian motion w(.). It is known that almost every sample path of r (t ) is a right-continuous step function with a finite number of simple jumps in any finite subinterval of R+ . As for r (t ), the following lemma is satisfied (see [12]).

P

Lemma 2.1. Given h > 0, then probability matrix P (h) = Pij (h)


N ×N



rnh = r (nh) , n = 0, 1, 2, . . .



is a discrete Markov chain with the one-step transition

= ehΓ .

Since the γij are independent of x, the paths of Markov chain r can be generated independent of x and, in fact, before computing x. We impose the following standing hypotheses for the functions f : Rn × Rn × R+ × S → Rn , g : Rn × Rn × R+ × S → Rn×m . (H1 ): Both f and g satisfy the local Lipschitz condition. That is, for each k = 1, 2, . . . , there is an hk > 0 such that

|f (x, y, t , i) − f (¯x, y¯ , t , i)| + |g (x, y, t , i) − g (¯x, y¯ , t , i)| ≤ hk (|x − x¯ | + |y − y¯ |) for all t ≥ 0 and i ∈ S and those x, y, x¯ , y¯ ∈ Rn with |x| ∨ |y| ∨ |¯x| ∨ |¯y| ≤ k. (H2 ): For every i ∈ S, there are constants αi ∈ R and ρi , σi , βi ≥ 0 such that 2xT f (x, 0, t , i) ≤ αi |x|2 ,

|f (x, 0, t , i) − f (x, y, t , i)| ≤ ρi |y| , |g (x, y, t , i)| ≤ σi |x| + βi |y|2 2

2

for all x, y ∈ Rn and t ≥ 0. In the scalar test Eq. (1), we can match in the following form: f (x(t ), x(t − τ ), t , r (t )) = A(r (t ))x(t ) + B(r (t ))x(t − τ ) + C (r (t ))

Z

t

x(s)ds t −τ

and g (x(t ), x(t − τ ), t , r (t )) = D(r (t ))x(t ) + E (r (t ))x(t − τ ).

1258

A. Rathinasamy, K. Balachandran / Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

If we take αi = 2A0 (i), ρi = |B(i)|+τ |C (i)|, σi = |D(i)| (|D(i)| + |E (i)|) and βi = |E (i)| (|D(i)| + |E (i)|), then the hypotheses (H1 ) and (H2 ) are satisfied naturally for Eq. (1). Then, for any given initial data ψ ∈ CFb 0 ([−τ , 0] : Rn ), Eq. (1) has a unique continuous solution denoted by x(t ; ψ) on t ≥ −τ . It is also easy to see that the equation admits the trivial solution x(t ; 0) ≡ 0. By means of Lyapunov’s function method [9], we obtain the following lemma: Theorem 2.1. Define

A = diag (− (α1 + ρ1 + σ1 ) , . . . , − (αN + ρN + σN )) − Γ .

(2)

Assume that A is a non-singular M-matrix and



1

(ρ1 + β1 )

1

,...,

T

(ρN + βN )

 A−1 1E,

(3)

E = (1, 1, . . . , 1)T . Then the trivial solution of Eq. (1) is exponentially stable in mean squares. where 1 3. Milstein method The general nonlinear multidelay systems without Markovian switching is dx(t ) = b

(0)

(x(t ), x(t − τ ), t )dt +

x(t ) = ψ(t ),

d X

(j)

(j)

b (x(t ), x(t − τ ), t )dw (t ),

t ≥ 0,

j =1

t ∈ [−τ , 0].

  

(4)

 

Küchler and Planten [16] presented the so-called order 1 strong Taylor approximation formula for (4). Applying the formula to the linear one-delay integro-differential system with Markovian switching, then we get the Milstein method to Eq. (1) leading to a numerical process of the following type: y¯ n+1 = y¯ n +

A

rnh




y¯ n + B

rnh




rnh

y¯ n−m + C

m
X

h

! y¯ n−k

h + D rnh y¯ n + E rnh y¯ n−m ∆wn










k=1

+D

rnh




D

rnh




y¯ n + E

rnh




y¯ n−m I1 + E rnh







D rnh y¯ n−m + E rnh y¯ n−2m I2 ,










(5)

where h > 0 is a stepsize which satisfies τ = mh for some positive integer m and tn = nh, = r (tn ) ∈ S, y¯ n is an approximation to x(tn ) = xn , if tn ≤ 0, we have y¯ n = ψ(tn ). Moreover, the increments ∆wn := w(tn+1 ) − w(tn ) are independent N (0, h)-distributed Gaussian random variables. Further, we assume that y¯ n is Ftn -measurable at the mesh points tn . Let I1 and I2 denote the two double integrals defined, respectively, by rnh

Z

tn+1

Z

s

dw(t )dw(s) =

I1 = tn

  (∆wn )2 − h

tn

2

and

Z

tn+1

Z

s

dw(t − τ )dw(s).

I2 = tn

tn

We will refer to the numerical scheme (5) as the Milstein method since the scheme is just the Milstein method when applied to a system without delay and Markovian switching. The convergence of method (5) without Markovian switching is discussed in [16]. Küchler and Planten [16] proved that the order 1 strong Taylor approximation formula converges strongly with order 1 whenever the coefficients b( j) (j = 1, 2, . . . , d) of system (4) without Markovian switching are homogeneous and satisfy both the generalized Lipschitz condition and the generalized growth condition. A simple check shows that the Milstein method (5) satisfies these conditions and hence it is strongly convergent of order 1. The following lemma [7] will be key to the proof. Lemma 3.1. The double integrals I1 and I2 have the same expectation and variation and satisfy E [I1 ] = E [I2 ] = 0,

E [I12 ] = E [I22 ] =

h2 2

,

E [ I1 I2 ] = 0 .

4. Mean-square stability of numerical scheme We will investigate the mean-square stability of the Milstein method in this section.

A. Rathinasamy, K. Balachandran / Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

1259

Definition 4.1. Under the conditions (3) and for any i ∈ S, Ai + |Bi | + τ |Ci | +

1 2

(|Di | + |Ei |)2 < 0,

(6)

where Ai = A (i), Bi = B (i), Ci = C (i), Di = D (i) and Ei = E (i), a numerical method is said to be mean square stable (MS-stable), if there exists a h0 (Ai , . . .) > 0 such that any application of the method to problem (1) generates a numerical approximation {¯yn }, which satisfies lim E |¯yn |2 = 0

(7)

n→∞

for all h ∈ (0, h0 (Ai , . . .)) with h = τ /m for an integer m. Now, we will prove the main theorem of this paper. Theorem 4.1. Assume the conditions (3) and (6) are satisfied. Then, the Milstein method applied to Eq. (1) is MS-stable. Proof. By re-arranging the right-hand side of (5), we have m X

y¯ n+1 = y¯ n + A rnh h + D rnh ∆wn y¯ n + C rnh h2













y¯ n−k + B rnh h + E rnh ∆wn y¯ n−m










k=1

+ D rnh




D rnh y¯ n + E rnh y¯ n−m I1 + E rnh













D rnh y¯ n−m + E rnh y¯ n−2m I2 .










By Lemma 2.1, the generation of rnh occurs before computing y¯ n+1 , then rnh is known. Since rnh ∈ S, so the above expression can be rewritten to y¯ n+1 = (1 + Ai h + Di ∆wn ) y¯ n + Ci h2

m X

y¯ n−k + (Bi h + Ei ∆wn ) y¯ n−m

k=1

+ Di (Di y¯ n + Ei y¯ n−m ) I1 + Ei (Di y¯ n−m + Ei y¯ n−2m ) I2 , where Ai = A (i), Bi = B (i), Ci = C (i), Di = D (i) and Ei = E (i) for any i ∈ S. Squaring both sides of the above equality, we have

" y2n+1

¯

= (1 + Ai h + Di ∆wn ) ¯ + Ci h 2

y2n

2

m X

#2 y¯ n−k

+ (Bi h + Ei ∆wn )2 y¯ 2n−m + D2i (Di y¯ n + Ei y¯ n−m )2 I12

k=1

+ Ei2 (Di y¯ n−m + Ei y¯ n−2m )2 I22 + 2 (1 + Ai h + Di ∆wn ) Ci h2

m X

y¯ n−k y¯ n

k=1

+ 2 (1 + Ai h + Di ∆wn ) (Bi h + Ei ∆wn ) y¯ n−m y¯ n + 2 (1 + Ai h + Di ∆wn ) Di (Di y¯ n + Ei y¯ n−m ) I1 y¯ n m X + 2 (1 + Ai h + Di ∆wn ) Ei (Di y¯ n−m + Ei y¯ n−2m ) I2 y¯ n + 2Ci h2 (Bi h + Ei ∆wn ) y¯ n−k y¯ n−m k=1

+ 2Ci h2 Di (Di y¯ n + Ei y¯ n−m ) I1

m X

y¯ n−k + 2Ci h2 Ei (Di y¯ n−m + Ei y¯ n−2m ) I2

k=1

m X

y¯ n−k

k=1

+ 2 (Bi h + Ei ∆wn ) y¯ n−m Di (Di y¯ n + Ei y¯ n−m ) I1 + 2 (Bi h + Ei ∆wn ) y¯ n−m Ei (Di y¯ n−m + Ei y¯ n−2m ) I2 + 2Di (Di y¯ n + Ei y¯ n−m ) Ei (Di y¯ n−m + Ei y¯ n−2m ) I1 I2 . Pm It follows from that 2xy ≤ x2 + y2 (∀x, y ∈ R), k=1 y¯ n−k ≤ m max(1≤k≤m) y¯ n−k and mh = τ , we get y¯ 2n+1 ≤ (1 + Ai h + Di ∆wn )2 y¯ 2n + Ci2 τ 2 h2 max y¯ 2n−k + (Bi h + Ei ∆wn )2 y¯ 2n−m (1≤k≤m)

D2i



 

+ + |Ei Di | y¯ 2n + Ei2 + |Ei Di | y¯ 2n−m I12 + Ei2 D2i + |Ei Di | y¯ 2n−m + Ei2 + |Ei Di | y¯ 2n−2m I22

+ |1 + Ai h| |Bi | h y¯ 2n + y¯ 2n−m + |Di Ei | (∆wn )2 y¯ 2n + y¯ 2n−m   + 2 [(1 + Ai h) Ei + Bi Di h] ∆wn y¯ n−m y¯ n + |1 + Ai h| |Ci | τ h max y¯ 2n−k + y¯ 2n (1≤k≤m)   + |Ci Bi | τ h2 max y¯ 2n−k + y¯ 2n−m + 2Ci Ei τ h∆wn max y¯ n−k y¯ n−m + 2Di Ci τ h∆wn max y¯ n−k y¯ n 

D2i

(1≤k≤m)

(1≤k≤m)

(1≤k≤m)

+ 2 (1 + Ai h + Di ∆wn ) Di (Di y¯ n + Ei y¯ n−m ) I1 y¯ n + 2 (1 + Ai h + Di ∆wn ) Ei (Di y¯ n−m + Ei y¯ n−2m ) I2 y¯ n + 2Ci τ hDi (Di y¯ n + Ei y¯ n−m ) I1 max y¯ n−k + 2Ci τ hEi (Di y¯ n−m + Ei y¯ n−2m ) I2 max y¯ n−k (1≤k≤m)

(1≤k≤m)

1260

A. Rathinasamy, K. Balachandran / Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

+ 2 (Bi h + Ei ∆wn ) y¯ n−m Di (Di y¯ n + Ei y¯ n−m ) I1 + 2 (Bi h + Ei ∆wn ) y¯ n−m Ei (Di y¯ n−m + Ei y¯ n−2m ) I2 + 2Di (Di y¯ n + Ei y¯ n−m ) Ei (Di y¯ n−m + Ei y¯ n−2m ) I1 I2 . (8)   2 Note that E (∆wn ) = 0, E (∆wn ) = h and any two of ∆wn , I1 and I2 are independent. Furthermore, y¯ n ,y¯ n−1 , . . . , y¯ n−m and y¯ n−2m are all Ftn -measurable. Hence, E (∆wn y¯ n y¯ n−m ) = E (¯yn y¯ n−m E (∆wn |Ftn )) = 0, E (∆wn )2 y¯ 2j = E y¯ 2j E ((∆wn )2 |Ftn ) = hE (¯y2j ),








j ∈ {n, n − m} .

Similarly, it can be derived that E I12 y¯ 2j =





E I22 y¯ 2k =





h2 2 h2 2

E y¯ 2j ,

j ∈ { n, n − m } ,

E y¯ 2k ,

k ∈ {n − m, n − 2m}

   

and

h

i

i

E y¯ in1 y¯ n2−m ∆wnj1 I1 = 0,

i1 , i2 ∈ {0, 1, 2} , j1 ∈ {0, 1} .

With a recursive calculation and Lemma 3.1, we obtain

h

i

i

j

i

E y¯ in1 y¯ n2−m y¯ n3−2m ∆wnj1 I12 I2 = 0,

i1 , i2 ∈ {0, 1, 2} , i3 , j1 , j2 ∈ {0, 1} ,

and j1 + j2 ≤ 1. Let Yn = E |yn |2 and taking expectations on both sides of inequality (8). Then, it holds that Yn+1 ≤ P (Ai , . . .) Yn + Q (Ai , . . .) Yn−m + R (Ai , . . .) Yn−2m + U (Ai , . . .) max Yn−k , 1≤k≤m

where P (Ai , . . .) = (1 + Ai h)2 + D2i h + |1 + Ai h| |Bi | h + |Di Ei | h + Q (Ai , . . .) = B2i h2 + Ei2 h + |1 + Ai h| |Bi | h + |Di Ei | h + R (Ai , . . .) =

h2

h2 2

h2 2

D2i D2i + |Di Ei | + |1 + Ai h| |Ci | τ h,




D2i Ei2 + |Di Ei | +




h2 2

Ei2 D2i + |Di Ei | + |Bi Ci | τ h2 ,




+ |Di Ei | , 2 U (Ai , . . .) = Ci2 τ 2 h2 + |1 + Ai h| |Ci | τ h + |Bi Ci | τ h2 . Ei2

Ei2




This implies that Yn+1 ≤ [P (Ai , . . .) + Q (Ai , . . .) + R (Ai , . . .) + U (Ai , . . .)] max {Yn−k , Yn−2m } . 0≤k≤m

By the above inequality we conclude that Yn → 0 (n → ∞) if P (Ai , . . .) + Q (Ai , . . .) + R (Ai , . . .) + U (Ai , . . .) ≤ 1, that is,



(1 + Ai h)2 + [|Bi | + τ |Ci |]2 h2 + (|Di | + |Ei |)2 h + 2 |1 + Ai h| [|Bi | + τ |Ci |] h +

h2 2

D2i

+

Ei2




(|Di | + |Ei |)

2



< 1.

Write

(

)   − 2Ai + 2 |Bi | + 2 |Ci | τ + (|Di | + |Ei |)2
h1 (Ai , . . .) = min i (|Ai | + |Bi | + |Ci | τ )2 + 21 D2i + Ei2 (|Di | + |Ei |)2 and



h2 (Ai , . . .) = min min i



1

|Ai |



 , min {T (i)} , i

(9)

A. Rathinasamy, K. Balachandran / Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

Fig. 1. Numerical simulation of case 1 with h =

1 26

and h =

1 25

1261

.

where T (i) =

  − 2Ai + 2 |Bi | + 2 |Ci | τ + (|Di | + |Ei |)2
. (Ai + |Bi | + |Ci | τ )2 + 12 D2i + Ei2 (|Di | + |Ei |)2

It follows from that condition (6) that h1 (Ai , . . .) > 0 and h2 (Ai , . . .) > 0. If h ∈ (0, h1 (Ai , . . .)), then we have

    
1 2 1 Di + Ei2 (|Di | + |Ei |)2 h2 + 2Ai + 2 |Bi | + 2 |Ci | τ + (|Di | + |Ei |)2 h < 0, (|Ai | + |Bi | + |Ci | τ )2 + 2

2

which implies that inequality (9) holds. If h ∈ (0, h2 (Ai , . . .)), then we have 1 + Ai h > 0 and

    
1 2 1 2 2 2 2 2 Di + Ei (|Di | + |Ei |) h + 2Ai + 2 |Bi | + 2 |Ci | τ + (|Di | + |Ei |) h < 0, (Ai + |Bi | + |Ci | τ ) + 2

2

which implies that inequality (9) still holds. Let h0 (Ai , . . .) = max {h1 (Ai , . . .) , h2 (Ai , . . .)} .

(10)

Then, we can conclude that inequality (9) holds whenever h ∈ (0, h0 (Ai , . . .)). This completes proof.



5. Numerical example We shall discuss an example to illustrate our theory. Let w(t ) be a scalar Brownian motion. Let r (t ) be a right continuous Markov chain taking values in S = {1, 2} with the generator

Γ = γij


2×2

 −1 = γ

1

−γ



.

Of course w(t ) and r (t ) are assumed to be independent. Consider one-dimensional linear stochastic delay integrodifferential equations with Markovian switching



dx(t ) = A (r (t )) x(t ) + B (r (t )) x(t − 1) + C (r (t ))

t

Z



x(s)ds dt + E (r (t )) x(t − 1)dw(t )

(11)

t −1

on t ≥ 0 with initial segment ψ(t ) = t + 1 for t ∈ [−1, 0] as a function, r (0) = 1. Case 1. Let A (1) = −2, B (1) = 1, C (1) = 12 , E (1) = 52 and A (2) = −9, B (2) = 4, C (2) = 3, E (2) = regard Eq. (11) as the result of the following two equations:



dx(t ) = −2x(t ) + x(t − 1) +

1

t

Z

2



x(s)ds dt + t −1

2 5

x(t − 1)dw(t )

4 . 5

It is interesting to

(12)

and



dx(t ) = −9x(t ) + 4x(t − 1) + 3

Z

t



x(s)ds dt + t −1

4 5

x(t − 1)dw(t )

(13)

1262

A. Rathinasamy, K. Balachandran / Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

Fig. 2. Numerical simulation of case 1 with h =

Fig. 3. Numerical simulation of case 1 with h =

1 24

and h =

1 22

and h =

1 23

.

1 . 2

switching from one to the other according to the movement of the Markov chain r (t ). It is known that Eqs. (12) and (13) are almost surely exponentially stable: however, as the result of Markovian switching, the overall behaviour, that is, Eq. (11), will be exponentially stable if 0 < γ < 2.354. So we take γ = 1. In the following tests, we show the influence of stepsize h on MS-stability of the Milstein method. The data used in all P100 1 figures are obtained by the mean square of data by 100 trajectories, that is, ωr : 1 ≤ r ≤ 100, yn = 100 yn (ωr )|2 . r =1 |¯ In all figures tn denotes the mesh point. In this case h0 (Ai , Bi , Ci , Di , Ei ) =

1 . 9

Thus, when a method of form (5) with

stepsize h : 0 < h < is used to solve the above system, the corresponding numerical solution is mean-square stable by Theorem 4.1. By the numerical tests, we can intuitively see the stepsize’s influence on the stability of the method. In the Milstein method by taking stepsizes h = 216 , h = 215 , h = 214 and h = 213 , respectively, we obtain four groups of numerical solutions of Eq. (11) on interval [0, 15], which are displayed in Figs. 1 and 2. One may
be interested to know about the stability of the numerical methods when stepsizes are outside the stability range 0, 19 . Here, we make an insight into the problem. Taking stepsizes h = 212 , h = 12 , respectively, we obtain two groups of numerical solutions of Eq. (11) on interval [0, 15], which are displayed in Fig. 3. 1 9

Case 2. A (1) = −2, B (1) = 1, C (1) = 0.5, E (1) = 0.4 and A (2) = 0.5, B (2) = 0.2, C (2) = 0.1, E (2) = 0.8. For these values, neither the condition of Theorem 2.1 nor the condition of Theorem 4.1 are satisfied. To carry out the numerical 1 simulation we choose the stepsize h = 1024 . The computer simulation result is shown in Fig. 4. Clearly, the Milstein method reveals the unstable property of the solution.

A. Rathinasamy, K. Balachandran / Nonlinear Analysis: Hybrid Systems 2 (2008) 1256–1263

Fig. 4. Numerical simulation of case 2 with h =

1 210

1263

.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]

L. Abolnikov, J.H. Dshalalow, A Treerattrakoon, On a dual hybrid queueing system, Nonlinear Analysis: Hybrid Systems 2 (2008) 96–109. D.J. Higham, Mean-square and asymptotic stability of the stochastic theta method, SIAM Journal of Numerical Analysis 38 (2000) 753–769. P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 1992. G.N. Milstein, Numerical Integration of Stochastic Differential Equations, Kluwer, Dordrecht, 1995. W.R. Cao, M.Z. Liu, Z.C. Fan, MS-stability of the Euler–Maruyama method for stochastic differential delay equations, Applied Mathematics and Computation 159 (2004) 127–135. M.Z. Liu, W.R. Cao, Z.C. Fan, Convergence and stability of the semi-implicit Euler method for a linear stochastic differential delay equation, Journal of Computational and Applied Mathematics 170 (2004) 255–268. Z. Wang, C. Zhang, An analysis of stability of Milstein method for stochastic differential equations with delay, Computers and Mathematics with Applications 51 (2006) 1445–1452. X. Mao, A. Matasov, A.B. Piunovkiy, Stochastic differential delay equations with Markovian switching, Bernoulli 6 (2000) 73–90. X. Mao, Robustness of stability of stochastic differential delay equations with Markovian switching, SACTA 3 (2000) 48–61. X. Mao, C. Yuan, Stochastic Differential Equations with Markovian Switching, Imperial College Press, London, 2006. A. Rathinasamy, K. Balachandran, Mean-square stability of second-order Runge–Kutta methods for multi-dimensional linear stochastic differential systems, Journal of Computational and Applied Mathematics 219 (2008) 170–197. L. Ronghua, H. Yingmin, Convergence and stability of numerical solutions to SDDEs with Markovian switching, Applied Mathematics and Computation 175 (2006) 1080–1091. T. Koto, Stability of θ -methods for delay integro-differential equations, Journal of Computational and Applied Mathematics 161 (2003) 393–404. L.E. Shaikhet, J.A. Rorberts, Reliability of difference analogues to preserve stability properties of stochastic Volterra integro-differential equations, Advances in Difference Equations (2006). X. Ding, K. Wu, M. Liu, Convergence and stability of the semi-implicit Euler method for linear stochastic delay integro-differential equations, International Journal of Computer Mathematics 83 (2006) 753–763. U. Küchler, E. Platen, Strong discrete time approximation of stochastic differential equations with time delay, Mathematics and Computer Simulation 54 (2000) 189–205.