- Email: [email protected]
Convergence and stability of the semi-tamed Euler scheme for stochastic differential equations with non-Lipschitz continuous coefficients
Applied Mathematics and Computation 228 (2014) 240–250
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Convergence and stability of the semi-tamed Euler scheme for stochastic differential equations with non-Lipschitz continuous coefficients Xiaofeng Zong 1, Fuke Wu ⇑,2, Chengming Huang 3 School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan, Hubei 430074, PR China
a r t i c l e
i n f o
Keywords: Stochastic differential equations Semi-tamed Euler scheme Strong convergence Mean square stability
a b s t r a c t Recently, explicit tamed schemes were proposed to approximate the SDEs with the non-Lipschitz continuous coefficients. This work proposes a semi-tamed Euler scheme, which is also explicit, to solve the SDEs with the drift coefficient equipped with the Lipschitz continuous part and non-Lipschitz continuous part. It is shown that the semitamed Euler converges strongly with the standard order one-half to the exact solution of the SDE. We also investigate the stability inheritance of the semi-tamed Euler schemes and reveal that this scheme does have advantage in reproducing the exponential mean square stability of the exact solution. Numerical experiments confirm the theoretical analysis. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction In this work, we study the numerical solution to the stochastic differential equations (SDEs) of the form
dxt ¼ lðxt Þdt þ rðxt ÞdW t ;
xð0Þ ¼ n;
ð1:1Þ d
dm
where W t is a m-dimensional Brownian motion, r : R # R is global Lipschitz continuous, lðxÞ ¼ f ðxÞ þ gðxÞ is globally one-side Lipschitz continuous, the function f : Rd # Rd is the Lipschitz continuous part of l, and g : Rd # Rd is non-Lipschitz continuous part. The SDE of this form does involve many stochastic models in the real world, such as stochastic Duffing–van der Pol oscillator [1,2], stochastic Lorenz equation [3,2], experimental psychology model [2], stochastic Ginzburg–Landau equation [4,2], stochastic Lotka–Volterra equations [5,4,2] and volatility processes [4,2], to name a few. Numerical analysis is an important tool in studying stochastic models, since most SDEs cannot be solved explicitly. In judging the quality of a numerical scheme, it is necessary to examine its convergence and stability. Since the drift coefficient is non-Lipschitz continuous, the classic explicit Euler–Maruyama (EM) method, investigated in Kloeden and Platen [4], Maruyama [6] and Milstein [7] for approximating the SDEs with globally Lipschitz continuous coefficients, may not converge in the strong mean square sense to the exact solution. Hutzenthaler et al. [8,9] shown that absolute moments of the explicit EM approximation for a SDE with a superlinearly growing and globally one-sided Lipschitz ⇑ Corresponding author. E-mail addresses: [email protected] (X. Zong), [email protected] (F. Wu), [email protected] (C. Huang). The research of this author was supported in part by the National Science Foundation of China (Grant No. 11001091). 2 The research of this author was supported in part by the National Science Foundation of China (Grant No. 11001091) and the Program for New Century Excellent Talents in University. 3 The research of this author was supported in part by the National Science Foundation of China (Grant Nos. 91130003 and 11371157) and the Fundamental Research Funds for the Central Universities (Grant No. 2013TS137). 1
0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.11.100
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
241
continuous drift coefficient, diverge to infinity at a finite time point. Although the implicit method, such as the backward Euler, split-step backward Euler methods, holds the strong convergence [10,11], it requires additional computation effort to solve a implicit system. Recently, Hutzenthaler et al. [12] proposed a new explicit approximation, called drift-tamed Euler method
Y Nnþ1 ¼ Y Nn þ
lðY Nn ÞT=N þ rðY Nn ÞDW Nn ; 1 þ klðY Nn ÞkT=N
ð1:2Þ
where DW Nn ¼ W ðnþ1ÞT=N W nT=N . It is shown that if the drift coefficient function is globally one-side Lipschitz continuous and has an at most polynomially growing derivative, the drift-tamed Euler method (1.2) converge strongly to the exact solution with the standard convergence order 0:5. Following the ideas, Wang and Gan [13] proposed the tamed Milstein method and proved it converge strongly to the exact solution with the convergence order 1. For the SDEs with non-Lipschitz coefficients, Hutzenthaler et al. [2] also proposed a series of the numerical approximations, such as partially drift-implicit approximation, linear implicit scheme, drift-truncated Euler and increment-tamed Euler, etc. Sauer and Stannat [14] proved strong convergence of the finite differences approximation in space for stochastic reaction diffusion equations with multiplicative noise under a one-sided Lipschitz condition using the tamed idea. Numerical stability for SDEs is primarily concerned with ascertaining for what values of stepsize does a particular numerical method replicate the stability properties of the exact solution. Under the globally Lipschitz condition, the explicit Euler– Maruyama scheme can reproduce the exponential mean square stability of the exact solution (see [15–19]), but for the SDEs (1.1) with non-Lipschitz coefficients, Higham et al. [20] gave a counterexample to show that the explicit Euler scheme may be explosive even though the exact solution is stable. Although the implicit scheme can share the stability of the exact solution (see [17,15,21]), it also requires additional computation effort to solve a implicit system. Therefore, it is necessary to seek for a new numerical scheme, which is explicit and can replicate the stability properties of the exact solution for the SDEs with non-Lipschitz continuous coefficients. Inspired by the literature [22,12,2], we derive a non-Lipschitz term tamed approximation for SDE (1.1):
Y Nnþ1 ¼ Y Nn þ f ðY Nn ÞT=N þ
gðY Nn ÞT=N 1 þ kgðY Nn ÞkT=N
þ rðY Nn ÞDW Nn :
ð1:3Þ
We refer to this approximation as semi-tamed Euler scheme. The increment function defined in [2] is /ðx; h; yÞ ¼ gðxÞ f ðxÞh þ 1þkgðxÞkh h þ rðxÞy. We will firstly show that the semi-tamed Euler scheme converges strongly to the exact solution
with the standard convergence order 0:5. Then we prove that the semi-tamed Euler scheme (1.3) can not only reproduce the exponential mean square stability of the exact solution, but also work better than the tamed Euler (1.2), backward Euler and split-step backward Euler schemes at stability preservation under a stepsize restriction. The rest of the paper is organized as follows. Section 2 begins with notations and estimation of the moment. Section 3 aims to show the strong convergence of the semi-tamed Euler scheme. Section 4 devotes to investigating the exponential mean square stability of the semi-tamed scheme for SDEs. Section 6 provides some numerical experiments for demonstration. 2. Estimation of the pth moments Throughout the whole article, unless otherwise specified, we use the following notations. Let T 2 ð0; 1Þ be a fixed real number and N 2 N be the step number of the uniform mesh with the stepsize T=N. Let ðX; F; PÞ be a complete probability space with a filtration fFt gtP0 satisfying the usual conditions and W t be a m-dimensional Brownian motion defined on this 1=2
probability space. Moreover, we use the notation kxk ¼ ðjx1 j2 þ þ jxk j2 Þ ; hx; yi ¼ x1 y1 þ x2 y2 þ þ xk yk for all x; y 2 Rk ; k 2 N, and kAk :¼ supx2Rl ; kxk61 kAxk for all A 2 Rkl ; k; l 2 N. a _ b represents maxfa; bg and a ^ b denotes minfa; bg. In this work, we also make the following assumptions. Assumption 2.1. Let gðxÞ be a continuously differentiable function and there exist positive constants K > 1 and c > 0, such that for any x; y 2 Rd
kf ðxÞ f ðyÞk _ krðxÞ rðyÞk 6 Kkx yk;
ð2:1Þ
hx y; lðxÞ lðyÞi 6 Kkx yk2 ;
ð2:2Þ
kg 0 ðxÞk 6 Kkxkc :
ð2:3Þ
Remark 2.1. There many stochastic models hold Assumption 2.1, such as stochastic Ginzburg–Landau equation, stochastic Verhulst equation, volatility processes (see [2]). Note that conditions (2.1) and (2.2) imply that SDE (1.1) admits the pth moment bounded solution in any finite time T (see [23]).
242
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
In order to prove the uniform boundedness of pth moments of the numerical solution (1.3), we follow the ideas in [12] to introduce the dominating stochastic processes DNn and appropriate sub-events XNn
DNn
! n1 X N 2 N ¼ ðk þ knkÞ exp k þ sup ½kkDW k k þ ak
ð2:4Þ
06u6n k¼u
and
( N n
) j sup DNk ð 06k6n1
x2X
X ¼
xÞ 6 N
2c
; sup kDW Nk k 06k6n1
61 ;
ð2:5Þ
where
k ¼ ð1 þ 3KT þ Tkf ð0Þk þ krð0Þk þ KÞ
4
and
aNn ¼ 1fkY Nn kP1g
* Y Nn þ f ðY Nn ÞT=N kY Nn k
;
rðY Nn Þ kY Nn k
+
DW Nn :
ð2:6Þ
The following lemmas are needed for proving the uniform boundedness of the pth moment. N N Lemma 2.1. Let Y N n ; Dn and Xn be given by (1.3), (2.4) and (2.5). Then
1XNn Y Nn 6 DNn
ð2:7Þ
for all n 2 f0; 1; 2; . . . ; Ng. Proof. Note that kDW Nn k 6 1 on XNnþ1 for all n 2 f0; 1; . . . ; N 1g and all N 2 N. Then Assumption 2.1 implies
kY Nnþ1 k 6 kY Nn k þ kf ðY Nn ÞkT=N þ kgðY Nn ÞkT=N þ krðY Nn ÞkkDW Nn k 6 1 þ Tkf ðY Nn Þ f ð0Þk þ Tkf ð0Þk þ TkgðY Nn Þk þ krðY Nn Þ rð0Þk þ krð0Þk 6 1 þ 3KT þ Tkf ð0Þk þ krð0Þk þ K 6 k N nþ1
ð2:8Þ
kY Nn k
on X \ fx 2 X : 6 1g for all n 2 f0; 1; 2; . . . ; Ng. Using the Cauchy–Schwarz inequality and the estimate ab 6 for all a; b 2 R, we have
kY Nnþ1 k2 ¼ kY Nn þ f ðY Nn ÞT=N þ
a2 2
2
þ b2
gðY Nn ÞT=N
þ rðY Nn ÞDW Nn k2 1 þ kgðY Nn ÞkT=N gðY N ÞT=N 2 2hY Nn ; gðY Nn ÞT=Ni N 2 N 2 2 2 N N 2 N N n ¼ kY n k þ kf ðY n Þk T =N þ þ krðY n ÞDW n k þ 2hY n ; f ðY n ÞT=Ni þ 1 þ kgðY Nn ÞkT=N 1 þ kgðY Nn ÞkT=N þ
2hf ðY Nn ÞT=N; gðY Nn ÞT=Ni 1þ
kgðY Nn ÞkT=N
þ 2hY Nn þ f ðY Nn ÞT=N; rðY Nn ÞDW Nn i þ
hgðY Nn ÞT=N; rðY Nn ÞDW Nn i 1 þ kgðY Nn ÞkT=N
6 kY Nn k2 þ 2kf ðY Nn Þk2 T 2 =N2 þ 3kgðY Nn ÞT=Nk2 þ 2krðY Nn ÞDW Nn k2 þ 2hY Nn ; f ðY Nn ÞT=Ni þ 2hY Nn þ f ðY Nn ÞT=N; rðY Nn ÞDW Nn i þ
2hY Nn ; lðY Nn ÞT=Ni 1 þ kgðY Nn ÞkT=N
2hY Nn ; f ðY Nn ÞT=Ni 1 þ kgðY Nn ÞkT=N
on X for all n 2 f0; 1; . . . ; N 1g and all N 2 N. The globally Lipschitz condition of f and condition of l give that 2
2
ð2:9Þ
:
r as well as the one side Lipschitz
2
kf ðxÞk2 6 ðkf ðxÞ f ð0Þk þ kf ð0ÞkÞ 6 ðKkxk þ kf ð0ÞkÞ 6 ðK þ kf ð0ÞkÞ kxk2 ;
ð2:10Þ
pffiffiffi jhx; f ðxÞij 6 kxkkf ðxÞ f ð0Þk þ kxkkf ð0Þk 6 Kkxk2 þ kf ð0Þkkxk 6 ðK þ kf ð0ÞkÞkxk2 6 kkxk2 ;
ð2:11Þ
krðxÞk2 6 ðkrðxÞ rð0Þk þ krð0ÞkÞ2 6 ðKkxk þ krð0ÞkÞ2 6 ðK þ krð0ÞkÞ2 kxk2 ; 6 kkxk2 ;
ð2:12Þ
pffiffiffi hx; lðxÞi 6 Kkxk2 þ kxkklð0Þk 6 kkxk2
ð2:13Þ
d
for all x 2 R with kxk P 1 and the polynomial bound on g 0 yields that
pffiffiffi kgðxÞk2 6 2K 2 kxk2ðcþ1Þ 6 N kkxk2
ð2:14Þ
243
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
for all x 2 Rd with 1 6 kxk 6 N 1=ð2cÞ and all N 2 N. Firstly, resorting to (2.11) and (2.13) as well as the definition of aNn , we have
2hY Nn þ f ðY Nn ÞT=N; rðY Nn ÞDW Nn i ¼ 2aNn kY Nn k2 on fx 2 X : 1 6
2hY Nn ; 1þ
kY Nn k
6N
ðY Nn ÞT=Ni kgðY Nn ÞkT=N
l
kgðY Nn ÞkT=N
ð2:15Þ
g for all N 2 N, and
6
2jhY Nn ; f ðY Nn ÞT=Nij 1þ
1=2c
6
pffiffiffi T 2 k kY Nn k2 ; N
ð2:16Þ
pffiffiffi T 2 k kY Nn k2 : N
2
Note that 2ðK þ kf ð0ÞkÞ T 2 < k and 3T 2 þ 6T 6 we obtain
ð2:17Þ pffiffiffi k. Then substituting the inequalities (2.10)–(2.12), (2.14)–(2.16) into (2.9),
pffiffiffi pffiffiffi 2 kY Nnþ1 k2 6 kY Nn k2 ð1 þ 2ðK þ kf ð0ÞkÞ T 2 =N þ 3 kT 2 =N þ 2kkDW Nn k2 þ 6 kT=N þ 2an Þ ! pffiffiffi 2 2ðK þ kf ð0ÞkÞ T 2 þ ð3T 2 þ 6TÞ k ¼ kY Nn k2 1 þ þ 2kkDW Nn k2 þ 2aNn N 2k þ 2kkDW Nn k2 þ 2aNn 6 kY Nn k2 exp N
ð2:18Þ
on fx 2 X : 1 6 kY Nn k 6 N 1=2c g for all N 2 N. Then we can establish (2.7) by the mathematical induction. Repeating the proof Lemma 2.1 in [12] completes this proof. h Lemma 2.2. For all p P 1,
sup N2N; NP4kpT
" !# N 1 X < 1: E exp pk kDW Nk k2
ð2:19Þ
k¼0
Proof. This proof is similar to Lemma 3.3 in [12] with only different k. h Lemma 2.3. Let an be given by (2.6). Then for all p P 1
" !# n1 X < 1: sup supE sup exp pz aNk
z2f1;1g N2N
06n6N
ð2:20Þ
k¼0
P N Proof. Note that the discrete stochastic process z n1 k¼0 ak is an ðFnT=N Þn2f0;1;...;Ng -martingale for every z 2 f1; 1g and every P N N 2 N. Then one can easily verify that the discrete stochastic process expðz n1 k¼0 ak Þ is a positive ðFnT=N Þn2f0;1;...;Ng submartingale for every z 2 f1; 1g and every N 2 N. Using Doob’s maximal inequality (see [24]), we have
! n1 X N sup exp z ak n2f0;1;...;Ng k¼0
Lp ðX;RÞ
! n1 X p N 6 exp z ak p1 k¼0
ð2:21Þ
Lp ðX;RÞ
for all N 2 N; p 2 ð1; 1Þ and all z 2 f1; 1g. Moreover, we have that
" 2 # x þ f ðxÞT=N rðxÞ p2 T krðxÞ ðx þ f ðxÞT=NÞk2 p2 T krðxÞk2 kx þ f ðxÞT=Nk2 E pz1fkxkP1g ; DW Nk ¼ 1fkxkP1g 6 1fkxkP1g 4 kxk kxk N N kxk kxk4 2
6
p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞ2 N
for all x 2 Rd ; k 2 f0; 1; . . . ; N 1g; N 2 N; p 2 ð1; 1Þ and all z 2 f1; 1g. The proof of Lemma 3.4 in [12] gives
!
2 x þ f ðxÞT=N rðxÞ p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞ2 ; 6 exp E exp pz1fkxkP1g ; DW Nk kxk kxk N which also shows that 2
E½expðpzaNk ÞjFkT=N 6 exp
p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞ2 N
!
244
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
for all x 2 Rd ; k 2 f0; 1; . . . ; N 1g; N 2 N; p 2 ½1; 1Þ and all z 2 f1; 1g. Hence, we have
" !# " ! # n1 n2 X X ¼ E exp pz aNk E½expðpzaNN1 ÞjFðN1ÞT=N E exp pz aNk k¼0
k¼0
"
n2 X 6 E exp pz aNk
!#
k¼0
2
p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞÞ2 exp N
!
2
6 6 expðp2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞÞ2 Þ; which together with (2.21) implies
" !# n1 X <1 sup 2 f1; 1gsup E sup exp pz aNk z
N2N
06n6N
k¼0
d
for all x 2 R ; k 2 f0; 1; . . . ; N 1g; N 2 N; p 2 ð1; 1Þ and all z 2 f1; 1g. h The following two lemmas are based on Lemmas 2.2 and 2.3 and can be proved by following the proof Lemma 3.5 and Lemma 3.6 in [12]. Lemma 2.4. Let DN n be given by (2.4). Then for all p 2 ½1; 1Þ
sup N2N; NP8kpT
E sup jDNn jp < 1:
ð2:22Þ
06n6N
Lemma 2.5. Let XNN be given by (2.5). Then for all p 2 ½1; 1Þ c
supðNp P½ðXNN Þ Þ < 1:
ð2:23Þ
N2N
Now, we establish the moment boundedness theorem of the approximation (1.3). Theorem 2.6. Let Y N n be given by (1.3), then for all p 2 ð1; 1Þ
sup sup EkY Nn kp < 1:
ð2:24Þ
N2N 06n6N
Proof. Firstly, we rewrite (1.3) as the form
Y Nn ¼ n þ
n1 n1 X X T=N f ðY Nk Þ þ k¼0
k¼0
T=NgðY Nn Þ 1 þ kgðY Nn ÞkT=N
¼ n þ f ð0ÞnT=N þ rð0ÞW NnT=N þ
þ
n1 X
rðY Nk ÞDW Nk
k¼0
n1 n1 X X T=Nðf ðY Nk Þ f ð0ÞÞ þ k¼0
k¼0
T=NgðY Nn Þ 1þ
kgðY Nn ÞkT=N
þ
n1 X ðrðY Nk Þ rð0ÞÞDW Nk k¼0
for all n 2 f0; 1; . . . ; N 1g and all N 2 N. By the Burkholder–Davis–Gundy type inequality in [12, Lemma 3.8], we have
kY Nn kLp ðX;Rd Þ
X n1 N 6 knkLp ðX;Rd Þ þ kf ð0ÞknT=N þ kr þ T=N½f ðY k Þ f ð0Þ p d k¼0 L ðX;R Þ X X n1 n1 gðY Nk Þ þ T=N þ ½rðY Nk Þ rð0ÞDW Nk p d k¼0 1 þ kgðY Nk ÞkT=NLp ðX;Rd Þ k¼0 L ðX;R Þ !1=2 m n1 X X 2 N 6 knkLp ðX;Rd Þ þ kf ð0ÞknT=N þ p nT=N kri ð0Þk þ T=N½f ðY k Þ f ð0Þ p d i¼1 k¼0 L ðX;R Þ !1=2 n1 m XX þNþp T=Nkri ðY Nk Þ ri ð0Þk2Lp ðX;Rd Þ ; ð0ÞW NnT=N kLp ðX;Rd Þ
k¼0 i¼1
which together with the global Lipschitz condition of f and
r gives
n1 X pffiffiffiffiffiffiffi kY Nn kLp ðX;Rd Þ 6 2ðknkLp ðX;Rd Þ þ kf ð0ÞkT þ p Tmkrð0Þk þ NÞ þ ðK þ 2p2 mK 2 Þ T=NkY Nn kLp ðX;Rd Þ
!
k¼0
for all n 2 f0; 1; . . . ; N 1g and all N 2 N. Resorting to Gronwall’s lemma, we have that for all N 2 N and all p 2 ½2; 1Þ
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
sup n2f0;1;...;Ng
kY Nn kLp ðX;Rd Þ 6 2ðknkLp ðX;Rd Þ þ C 0 þ NÞ exp ðK þ 2p2 mK 2 ÞT ;
245
ð2:25Þ
pffiffiffiffiffiffiffi where C 0 ¼ kf ð0ÞkT þ p Tmkrð0Þk. Note that the inequality above does not prove this theorem due to N 2 N on the righthand side of the inequality (2.25). However, exploiting (2.25) in an appropriate bootstrap argument will enable us to establish Theorem 2.6. Firstly,
sup sup k1ðXN Þc Y Nn kLp ðX;Rd Þ 6 sup sup ðk1ðXN Þc kL2p ðX;Rd Þ kY Nn kL2p ðX;Rd Þ Þ n
N2N 06n6N
n
N2N 06n6N
6 sup Nk1ðXN Þc kL2p ðX;Rd Þ sup sup N 1 kY Nn kL2p ðX;Rd Þ n
N2N
N2N 06n6N
1=2p c 6 sup N2p P½ðXNn Þ sup sup N1 kY Nn kL2p ðX;Rd Þ < 1; N2N
ð2:26Þ
N2N 06n6N
where we used Hölder’s inequality, estimate (2.25) and Lemma 2.5. In addition, By Lemmas 2.1 and 2.4, we obtain
sup sup k1XNn Y Nn kLp ðX;Rd Þ 6 sup sup kDNn kLp ðX;Rd Þ < 1; N2N 06n6N
N2N 06n6N
which together with (2.26) implies the desired assertion.
h
3. Strong convergence of the semi-tamed Euler scheme In order to establish our convergence theorem for the semi-tamed Euler method (1.3), we now introduce appropriate time continuous interpolations of the time discrete numerical approximations (1.3). More accurately, we define the time continuous approximation yN : ½0; T X ! Rd ; N 2 N given by
t nT gðY Nn Þ nT N f ðY Nn Þ þ yNt ¼ Y Nn þ t þ rðY Nn ÞðW t W nT=N Þ N 1 þ kgðY Nn ÞkT=N
ð3:1Þ
for all t 2 ½nT ; ðnþ1ÞT Þ; n 2 f0; 1; . . . ; N 1g and all N 2 N. It is easy to see that yN is Ft -adapted stochastic process. N N Theorem 3.1. Let Assumption 2.1 hold. Then there exists a family C p 2 ð0; 1Þ; p 2 ½1; 1Þ, of real numbers such that
"
#!1=p
E sup kxt yNt kp
6 C p N1=2
ð3:2Þ
t2½0;T
for all N 2 N and all p 2 ½1; 1Þ. Before proving this theorem, we firstly present a lemma, which will be used in the proof of Theorem 3.1. Lemma 3.2. Let Y N n be defined by (1.3). Then we have
sup sup Ekf ðY Nn Þkp < 1; N2N n2f0;...;Ng
sup sup EkgðY Nn Þkp < 1 N2N n2f0;...;Ng
and
sup sup EkrðY Nn Þkp < 1 N2N n2f0;...;Ng
for all p 2 ½1; 1Þ. Proof. This proof can be completed by Assumption 2.1 and Theorem 2.6. h
kgðxÞkgðxÞ ðxÞ ¼ lðxÞ T=N 1þkgðxÞkT=N Proof of Theorem 3.1 Define yNs ¼ Y Nn for s 2 ½nT ; ðnþ1ÞT Þ for n 2 f0; 1; . . . ; N 1g; N 2 N; l and rewrite N N the continuous approximation (3.1) as
yNs ¼ n þ
Z
s
0
l ðyNu Þdu þ
Z
s 0
rðyNu ÞdW u :
Then we obtain from (1.1) and (3.3)
xs yNs ¼
Z 0
s
ðyNu Þdu þ ½ðlðxu Þ l
Z 0
s
½rðxu Þ rðyNu ÞdW u
ð3:3Þ
246
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
for all s 2 ½0; T and N 2 N. Itô’s formula produces
+ Z s* kgðyNu ÞkgðyNu Þ du xu yNu ; 1 þ kgðyNu ÞkT=N 0 0 m Z s m Z s X X i þ2 hxu yNu ; ri ðxu Þ ri ðyNu ÞidW u þ kri ðxu Þ ri ðyNu Þk2 du
kxs yNs k2 ¼ 2
Z
s
hxu yNu ; lðxu Þ lðyNu Þidu þ 2T=N
0
i¼1
¼2
Z
hxu yNu ; lðxu Þ lðyNu Þidu þ 2
0
þ
i¼1
s
m Z s X 0
i¼1
kri ðxu Þ ri ðyNu Þk2 du þ 2
Z 0
hxu yNu ; lðyNu Þ lðyNu Þidu þ 2T=N
m Z s X 0
i¼1
2
0
s
+ Z s* kgðyNu ÞkgðyNu Þ xu yNu ; du 1 þ kgðyNu ÞkT=N 0
i
hxu yNu ; ri ðxu Þ ri ðyNu ÞidW u :
2
By the inequality ða þ bÞ 6 2a2 þ 2b for all a; b 2 R, the Cauchy–Schwarz inequality and Assumption 2.1, we have
kxs yNs k2 6 ð2K þ 2K 2 mÞ þ 2K 2 m
Z 0
s
Z
s 0
Z
kxu yNu k2 du þ 2
kyNu yNu k2 du þ 2
m Z X i¼1
0
s
0
kxu yNu k klðyNu Þ lðyNu Þkdu þ 2T=N
s
Z 0
s
kxu yNu k kgðyNu Þk2 du
i
hxu yNu ; ri ðxu Þ ri ðyNu ÞidW u 2
for all s 2 ½0; T and all N 2 N. The inequality ab 6 a2 =2 þ b =2 for all a; b 2 R hence yields that
Z t t kxu yNu k2 du þ klðyNu Þ lðyNu Þk2 du þ T 2 =N2 kgðyNu Þk4 du 0 0 0 Z t X m Z t i 2 N N 2 N N þ 2K m kyu yu k du þ 2 sup hxu yu ; ri ðxu Þ ri ðyu ÞidW u ; 06s6t i¼1 0 0
sup kxs yNs k2 6 2ðK þ K 2 m þ 1Þ
06s6t
Z
Z
t
which together with the Burkholder–Davis–Gundy type inequality in [12, Lemma 3.7] implies
Z t Z t N 2 sup xs yN k2 k p=2 d 6 2ðK þ K 2 m þ 1Þ kx y k du þ klðyNu Þ lðyNu Þk2Lp ðX;Rd Þ du p d u s u L ðX;R Þ L ðX;R Þ 06s6t 0 0 Z t Z t þ T 2 =N2 kgðyNu Þk4L2p ðX;Rd Þ du þ 2K 2 m kyNu yNu k2Lp ðX;Rd Þ du 0
þp
m Z X i¼1
0
t 0
!1=2
khxu yNu ; ri ðxu Þ ri ðyNu Þik2Lp=2 ðX;RÞ du
ð3:4Þ
for all t 2 ½0; T; N 2 N and p 2 ½4; 1Þ. The last term in (3.4) has been estimated in the proof Theorem 1.1 of [12] as
p
m Z X i¼1
0
t
!1=2
khxu yNu ; ri ðxu Þ ri ðyNu Þik2Lp=2 ðX;RÞ du 2
6
2 Z 1 p2 K 2 m t sup kxs yN k þ kxs yNs k2Lp ðX;Rd Þ ds: s 2 06s6t 2 p d 0 L ðX;R Þ
ð3:5Þ
2
Using the inequality ða þ bÞ 6 2a2 þ 2b for all a; b 2 R and inserting (3.5) into (3.4) produce
2 sup kxs yN k s p 06s6t
L ðX;Rd Þ
!Z Z t t p2 K 2 m 62 K þK mþ1þ kxu yNu k2Lp ðX;Rd Þ du þ T 2 =N2 kgðyNu Þk4L2p ðX;Rd Þ du þ ð2K 2 m 2 0 0 2 Z t Z t 1 N sup kyNu yNu k2Lp ðX;Rd Þ du þ klðyNu Þ lðyNu Þk2Lp ðX;Rd Þ du þ kx y k ; þ p2 K 2 mÞ s s 2 06s6t 0 0 Lp ðX;Rd Þ 2
which implies
2 Z t Z t 1 2 2 2 2 N 2 sup kxs yN k 6 2ðK þ 1 þ p K mÞ kx y k du þ T =N kgðyNu Þk4L2p ðX;Rd Þ du p d u s u L ðX;R Þ 2 06s6t 0 0 Lp ðX;Rd Þ Z t Z t kyNu yNu k2Lp ðX;Rd Þ du þ klðyNu Þ lðyNu Þk2p þ ð2K 2 m þ p2 K 2 mÞ du L ðX;Rd Þ 0 0 for all t 2 ½0; T; N 2 N and p 2 ½4; 1Þ. The Gronwall lemma and the inequality a; b; c 2 ½0; 1Þ give
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi pffiffiffi pffiffiffi a þ b þ c 6 a þ b þ b for all
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
sup kxs yN k s 06s6T
6
pffiffiffiffiffiffi pffiffiffiffiffiffiffi 2T expð2K þ 2 þ 2p2 K 2 mÞ 2p2 K 2 2m sup kyNu yNu kLp ðX;Rd Þ 06u6T
Lp ðX;Rd Þ
247
þT=N sup kgðY Nn Þk2L2p ðX;Rd Þ n2f0;1;...;Ng
þ sup kl 06u6T
ðyNu Þ
ðyNu ÞkLp ðX;Rd Þ
l
! :
ð3:6Þ
Additionally, (3.3) and the Burkholder–Davis–Gundy type inequality show that for all N 2 N and p 2 ½2; 1Þ
sup 06u6T
kyNu
yNu kLp ðX;Rd Þ
Z u T T T N N N 6 sup kf ðyu ÞkLp ðX;Rd Þ þ sup kgðyu ÞkLp ðX;Rd Þ þ sup rðyu ÞdW u p d N 06u6T N 06u6T N 06u6T u L ðX;R Þ pffiffiffiffiffiffiffi T T p Tm N N sup sup kri ðY Nn ÞkLp ðX;Rd Þ ; 6 sup kf ðY n ÞkLp ðX;Rd Þ þ sup kgðY n ÞkLp ðX;Rd Þ þ pffiffiffiffi N 06n6N N 06n6N N i2f1;2;...;mg06n6N
which together with Lemma 3.2 gives
sup N2N
pffiffiffiffi N sup kyNu yNu kLp ðX;Rd Þ < 1
ð3:7Þ
06u6T
for all p 2 ½1; 1Þ. Particularly, we obtain
sup sup kyNu kLp ðX;Rd Þ < 1
ð3:8Þ
N2N 06u6T
for all p 2 ½1; 1Þ. Assumption 2.1 gives that
sup kyNu yNu kLp ðX;Rd Þ ; sup klðyNu Þ lðyNu ÞkLp ðX;Rd Þ 6 K 1 þ sup kyNu kcLp ðX;Rd Þ
06u6T
06u6T
06u6T
which together with (3.8) produces
sup N2N
pffiffiffiffi N sup klðyNu Þ lðyNu ÞkLp ðX;Rd Þ < 1
ð3:9Þ
06u6T
for all N 2 N and p 2 ½1; 1Þ. Combining (3.6), (3.7) and (3.9) and Lemma 3.2 gives (3.2). h 4. Exponential mean square stability of the semi-tamed Euler scheme In this section, we examine whether the semi-tamed scheme can share the exponential mean square stability of the exact solution. Since our aim involves the asymptotic behavior at the larger time, we rewrite the semi-tamed scheme (1.3) as follows:
Y nþ1 ¼ Y n þ f ðY n Þh þ
gðY n Þh þ rðY n ÞDW n ; 1 þ kgðY n Þkh
ð4:1Þ
where h > 0 is the stepsize and DW n ¼ W ðnþ1Þh W nh is the Brownian increment. For the stability purpose, we assume f ð0Þ ¼ gð0Þ ¼ 0 and rð0Þ ¼ 0, which shows that Eq. (1.1) admits a trivial solution. Firstly, we give the stability criterion of the exact solution (see [23]). Theorem 4.1. Let
l and r satisfy the local Lipschitz condition. If there exists a positive constant c such that
2hx; lðxÞi þ krðxÞk2 6 ckxk2 ;
8 x 2 Rd ;
ð4:2Þ
then the solution of Eq. (1.1) hold the property
Ekxt k2 6 Eknk2 ect :
ð4:3Þ
However, in order to examine stability of the numerical approximation Y n given by (4.1), we assume that there exist some with 2q > h2 and 2m > m such that nonnegative constants K; q; h; m; a > 1; m
hx y; f ðxÞ f ðyÞi 6 qkx yk2 ; hx y; gðxÞ gðyÞi 6 mkx ykaþ1 ;
kf ðxÞ f ðyÞk 6 Kkx yk; kxka kgðxÞk 6 m
ð4:4Þ ð4:5Þ
and
krðxÞ rðyÞk 6 hkx yk d
ð4:6Þ 2
for all x; y 2 R . Denote c ¼ 2q h > 0, then the three conditions above imply (4.3) by Theorem 4.1. Moreover, these conditions are also needed for the convergence.
248
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
Theorem 4.2. Let the conditions (4.4)–(4.6) hold. The there exist a stepsize h such that for any h < h , the discrete solution (4.1) holds the property
EkY n k2 6 Eknk2 enhch ;
ð4:7Þ
where ch > 0 satisfies limh!0 ch ¼ 2q h2 . Proof. The semi-tamed Euler approximation (4.1) gives 2
2
kY nþ1 k2 ¼ kY n k2 þ 2hY n ; f ðY n Þih þ kf ðY n Þk2 h þ krðY n ÞDW n k2 þ 2 þ 2 Y n þ f ðY n Þh þ
gðY n Þh ; rðY n ÞDW n : 1 þ kgðY n Þkh
hY n þ f ðY n Þh; gðY n Þhi kgðY n Þk2 h þ 2 1 þ kgðY n Þkh ð1 þ kgðY n ÞkhÞ ð4:8Þ
By the condition (4.5), we have 2
2
hkY n kaþ1 h hY n þ f ðY n Þh; gðY n Þhi mkY n kaþ1 h kf ðY n Þk kgðY n Þkh ½m K m 62 þ2 62 ; 1 þ kgðY n Þkh 1 þ kgðY n Þkh 1 þ kgðY n Þkh 1 þ kgðY n Þkh
ð4:9Þ
where K is the Lipschitz constant of f. Define Xn ¼ fx 2 X : kY n k > 1g, then we have 2
1Xn
kgðY n Þk2 h
ð1 þ kgðY n ÞkhÞ
2
6 1Xn
2
6 1Xcn
kgðY n Þkh mkY n kaþ1 h 6 1Xn 1 þ kgðY n Þkh 1 þ kgðY n Þkh
ð4:10Þ
and 2
1Xcn
kgðY n Þk2 h
ð1 þ kgðY n ÞkhÞ
m2 kY n k2a h2 1 þ kgðY n Þkh
6 1Xcn
m2 kY n kaþ1 h2 1 þ kgðY n Þkh
ð4:11Þ
:
mm ^ 2m . For any h < h , inserting (4.9)–(4.11) into (4.8) yields Define h1 ¼ 22K 1 Þm m ð2Kþm
2 kY nþ1 k2 6 kY n k2 þ 2hY n ; f ðY n Þih þ kf ðY n Þk2 h þ krðY n ÞDW n k2 þ 2 Y n þ f ðY n Þh þ
gðY n Þh ; rðY n ÞDW n : 1 þ kgðY n Þkh
ð4:12Þ
Taking expectations on the both sides of (4.12) gives that for any h < h1 , 2
EkY nþ1 k2 6 EkY n k2 þ 2E½hY n ; f ðY n Þih þ Ekf ðY n Þk2 h þ EkrðY n Þk2 h; where we used the independent increments of the Brownian motion W t . Hence, we obtain from (4.4) and (4.6) that 2
2
EkY nþ1 k2 6 EkY n k2 2qhEkY n k2 þ Kh EkY n k2 þ h2 hEkY n k2 ¼ ð1 2qh þ Kh þ h2 hÞEkY n k2 6 2
6 ð1 2qh þ Kh þ h2 hÞ
nþ1
Eknk2 6 Eknk2 ech ðnþ1Þh
2
for any h < h :¼ h1 ^ 2qKh . This completes the proof of Theorem 4.2. h
5. Numerical experiments This section presents some numerical experiments to illustrate our theoretical results. Let us firstly consider the scalar stochastic Ginzburg–Landau equation:
dxt ¼ ð2xt x3t Þdt þ xt dW t ;
t>0
ð5:1Þ
with the initial data x0 ¼ n ¼ 0:5. Then the drift lðxÞ ¼ 2x x3 do not satisfy the global Lipschitz condition, but can be decomposed into the Lipschitz continuous part f ðxÞ ¼ 2x and non-Lipschitz continuous part gðxÞ ¼ x3 . Fig. 1 depicts the root mean-square errors as a function of the stepsize h in log–log plot, where the expectation is approximated by the mean of 3 103 independent realizations. As expected, the semi-tamed Euler scheme gives an error that decreases proportional to 1=2 h . Then we investigate the numerical stability of the tamed scheme (1.3). We consider the scalar SDE:
dxt ¼ ð2xt x3t Þdt þ
pffiffiffi 2xt dW t ;
t>0
ð5:2Þ
It is easy to deduce from Theorem 4.1 that the solution (5.2) is exponentially mean square stable with the Lyapunov exponent not greater than 2, since
2xT lðxÞ þ jrðxÞj2 ¼ 4x2 2x4 þ 2x2 6 2x2 :
ð5:3Þ
249
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250 Convergence order
0
10
Semi−tamed Euler Slope of 0.5
Mean square error
−1
10
−2
10
−3
10
−3
−2
10
−1
10 Stepsize
10
Fig. 1. Mean square approximation error versus stepsize h to approximate (5.1).
Tamed Euler
Backward Euler
0
0
10
10
h = 1/4 h = 1/8 h = 1/16
2
E[Yn]
n
E[Y2]
h = 1/4 h = 1/8 h = 1/16
−10
−10
10
10
0
2
4
6
8
10
0
2
4
t
6
8
10
t
Split−step backward Euler
Semi−tamed Euler
0
0
10
10
h = 1/4 h = 1/8 h = 1/16
E[Yn]
2
n
E[Y2]
h = 1/4 h = 1/8 h = 1/16
−10
−10
10
10
0
2
4
6
t
8
10
0
2
4
6
8
10
t
Fig. 2. Computer simulation of Ejxt j2 , using different Euler schemes and different stepsizes.
We test the stability of the semi-tamed Euler, drift-tamed Euler, backward Euler and split-step backward Euler schemes with different stepsizes. We take h ¼ 1=4; 1=8 and 1=16 (h < h ¼ 0:4) and generate 7 103 sample paths for each numerical method. The mean square stability of the numerical solution is plotted in Fig. 2. Fig. 2 shows that the four numerical schemes are all mean square stable under the stepsizes h ¼ 1=4; 1=8 and 1=16, but the semi-tamed Euler scheme works better than the other three schemes. 6. Further remark The work proposes the semi-tamed Euler scheme to approximate the SDEs with non-Lipschitz continuous coefficients. We obtained its strong convergence order 0:5 and shown that the semi-tamed can reproduce the exponential mean square stability of the test equation. The numerical experiments also shown that the semi-tamed Euler scheme is superior to the tamed Euler and the backward Euler for reproducing mean square stability of the exact solution. We also meed to remark that (i) the conditions (2.3) and (4.5) on the non-Lipschitz continuous part g can be more complex, for example kg 0 ðxÞk 6 Kðkxka þ kxkb Þ; hx y; gðxÞ gðyÞi 6 t1 kx yk1þa t2 kx yk1þb and kgðxÞk 6 t1 kxka þ t2 kxkb ; (ii) the semi-tamed
250
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
idea can also be applied to the classical Milstein scheme and lead to the semi-tamed Milstein scheme, which can be considered as the extension of the tamed Milstein scheme developed in [13]. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2013.11.100. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
L. Arnold, Random dynamical systems, Springer Monographs in Mathematics, Springer-Verlag, Berlin, 1998. M. Hutzenthaler, A. Jentzen, Convergence of the stochastic Euler scheme for locally Lipschitz coefficients, Found. Comput. Math. 11 (6) (2011) 657–706. B. Schmalfuss, The random attractor of the stochastic Lorenz system, Zeitschrift fur Ange-wandte Mathematik und Physik 48 (1997) 951–975. P.E. Kloeden, E. Platen, Numerical solution of stochastic differential equations, Applications of Mathematics (New York), vol. 2, Springer-Verlag, Berlin, 1992. M. Arato, A famous nonlinear stochastic equation (Lotka–Volterra model with diffusion), Math. Comput. Model. 38 (7) (2003) 709–726. G. Maruyama, Continuous Markov processes and stochastic equations, Rend. Circ. Mat. Palermo 4 (2) (1955) 48–90. G.N. Milstein, Numerical integration of stochastic differential equations, Mathematics and Its Applications, vol. 313, Kluwer Academic, Dordrecht, 1995. M. Hutzenthaler, A. Jentzen, Non-globally Lipschitz counterexamples for the stochastic Euler scheme. arXiv:0905.0273v1. M. Hutzenthaler, A. Jentzen, P.E. Kloeden, Strong and weak divergence in finite time of Eulers method for stochastic differential equations with nonglobally Lipschitz continuous coefficients, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 467 (2011) 1563–1576. D.J. Higham, X. Mao, A.M. Stuart, Strong convergence of Euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40 (2) (2002) 1041–1063. X. Mao, L. Szpruch, Strong convergence and stability of implicit numerical methods for stochastic differential equations with non-globally Lipschitz continuous coefficients, J. Comput. Appl. Math. 238 (15) (2013) 14–28. M. Hutzenthaler, A. Jentzen, P.E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22 (4) (2012) 1611–1641. X. Wang, S. Gan, The tamed Milstein method for commutative stochastic differential equations with non-globally Lipschitz continuous coefficients, J. Differ. Equ. Appl. 19 (3) (2013) 466–490. M. Sauer, W. Stannat, Lattice approximation for stochastic reaction diffusion equations with one-sided Lipschitz condition. arXiv:1301.6350. D.J. Higham, X. Mao, A. Stuart, Exponential mean-square stability of numerical solutions to stochastic differential equations, LMS J. Comput. Math. 6 (2003) 297–313. D.J. Higham, Mean-square and asymptotic stability of the stochastic theta method, SIAM J. Numer. Anal. 38 (2000) 753–769. C. Huang, Exponential mean square stability of numerical methods for systems of stochastic differential equations, J. Comput. Appl. Math. 236 (2012) 4016–4026. Y. Saito, T. Mitsui, Stability analysis of numerical schemes for stochastic differential equations, SIAM J. Numer. Anal. 33 (1996) 2254–2267. Y. Saito, T. Mitsui, Mean-square stability of numerical schemes for stochastic differential systems, Vietnam J. Math. 30 (2002) 551–560. D.J. Higham, X. Mao, C. Yuan, Almost sure and moment exponential stability in the numerial simulation of stochastic differential equations, SIAM J. Numer. Anal. 45 (2007) 592–609. X. Zong, F. Wu, Choice of h and mean-square exponential stability in the stochastic theta method of stochastic differential equations, J. Comput. Appl. Math. 255 (2014) 837–847. H. Szpruch, X. Mao, Strong convergence of numerical methods for nonlinear stochastic differential equations under monotone conditions, Preprint, Technical Report. X. Mao, Stochastic differential equations and their applications, Horwood Publishing Series in Mathematics & Applications, Horwood Publishing Limited, Chichester, 1997. A. Klenke, Probability Theory: A Comprehensive Course, Springer, London, 2008.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Convergence and stability of the semi-tamed Euler scheme for stochastic differential equations with non-Lipschitz continuous coefficients Xiaofeng Zong 1, Fuke Wu ⇑,2, Chengming Huang 3 School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan, Hubei 430074, PR China
a r t i c l e
i n f o
Keywords: Stochastic differential equations Semi-tamed Euler scheme Strong convergence Mean square stability
a b s t r a c t Recently, explicit tamed schemes were proposed to approximate the SDEs with the non-Lipschitz continuous coefficients. This work proposes a semi-tamed Euler scheme, which is also explicit, to solve the SDEs with the drift coefficient equipped with the Lipschitz continuous part and non-Lipschitz continuous part. It is shown that the semitamed Euler converges strongly with the standard order one-half to the exact solution of the SDE. We also investigate the stability inheritance of the semi-tamed Euler schemes and reveal that this scheme does have advantage in reproducing the exponential mean square stability of the exact solution. Numerical experiments confirm the theoretical analysis. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction In this work, we study the numerical solution to the stochastic differential equations (SDEs) of the form
dxt ¼ lðxt Þdt þ rðxt ÞdW t ;
xð0Þ ¼ n;
ð1:1Þ d
dm
where W t is a m-dimensional Brownian motion, r : R # R is global Lipschitz continuous, lðxÞ ¼ f ðxÞ þ gðxÞ is globally one-side Lipschitz continuous, the function f : Rd # Rd is the Lipschitz continuous part of l, and g : Rd # Rd is non-Lipschitz continuous part. The SDE of this form does involve many stochastic models in the real world, such as stochastic Duffing–van der Pol oscillator [1,2], stochastic Lorenz equation [3,2], experimental psychology model [2], stochastic Ginzburg–Landau equation [4,2], stochastic Lotka–Volterra equations [5,4,2] and volatility processes [4,2], to name a few. Numerical analysis is an important tool in studying stochastic models, since most SDEs cannot be solved explicitly. In judging the quality of a numerical scheme, it is necessary to examine its convergence and stability. Since the drift coefficient is non-Lipschitz continuous, the classic explicit Euler–Maruyama (EM) method, investigated in Kloeden and Platen [4], Maruyama [6] and Milstein [7] for approximating the SDEs with globally Lipschitz continuous coefficients, may not converge in the strong mean square sense to the exact solution. Hutzenthaler et al. [8,9] shown that absolute moments of the explicit EM approximation for a SDE with a superlinearly growing and globally one-sided Lipschitz ⇑ Corresponding author. E-mail addresses: [email protected] (X. Zong), [email protected] (F. Wu), [email protected] (C. Huang). The research of this author was supported in part by the National Science Foundation of China (Grant No. 11001091). 2 The research of this author was supported in part by the National Science Foundation of China (Grant No. 11001091) and the Program for New Century Excellent Talents in University. 3 The research of this author was supported in part by the National Science Foundation of China (Grant Nos. 91130003 and 11371157) and the Fundamental Research Funds for the Central Universities (Grant No. 2013TS137). 1
0096-3003/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2013.11.100
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
241
continuous drift coefficient, diverge to infinity at a finite time point. Although the implicit method, such as the backward Euler, split-step backward Euler methods, holds the strong convergence [10,11], it requires additional computation effort to solve a implicit system. Recently, Hutzenthaler et al. [12] proposed a new explicit approximation, called drift-tamed Euler method
Y Nnþ1 ¼ Y Nn þ
lðY Nn ÞT=N þ rðY Nn ÞDW Nn ; 1 þ klðY Nn ÞkT=N
ð1:2Þ
where DW Nn ¼ W ðnþ1ÞT=N W nT=N . It is shown that if the drift coefficient function is globally one-side Lipschitz continuous and has an at most polynomially growing derivative, the drift-tamed Euler method (1.2) converge strongly to the exact solution with the standard convergence order 0:5. Following the ideas, Wang and Gan [13] proposed the tamed Milstein method and proved it converge strongly to the exact solution with the convergence order 1. For the SDEs with non-Lipschitz coefficients, Hutzenthaler et al. [2] also proposed a series of the numerical approximations, such as partially drift-implicit approximation, linear implicit scheme, drift-truncated Euler and increment-tamed Euler, etc. Sauer and Stannat [14] proved strong convergence of the finite differences approximation in space for stochastic reaction diffusion equations with multiplicative noise under a one-sided Lipschitz condition using the tamed idea. Numerical stability for SDEs is primarily concerned with ascertaining for what values of stepsize does a particular numerical method replicate the stability properties of the exact solution. Under the globally Lipschitz condition, the explicit Euler– Maruyama scheme can reproduce the exponential mean square stability of the exact solution (see [15–19]), but for the SDEs (1.1) with non-Lipschitz coefficients, Higham et al. [20] gave a counterexample to show that the explicit Euler scheme may be explosive even though the exact solution is stable. Although the implicit scheme can share the stability of the exact solution (see [17,15,21]), it also requires additional computation effort to solve a implicit system. Therefore, it is necessary to seek for a new numerical scheme, which is explicit and can replicate the stability properties of the exact solution for the SDEs with non-Lipschitz continuous coefficients. Inspired by the literature [22,12,2], we derive a non-Lipschitz term tamed approximation for SDE (1.1):
Y Nnþ1 ¼ Y Nn þ f ðY Nn ÞT=N þ
gðY Nn ÞT=N 1 þ kgðY Nn ÞkT=N
þ rðY Nn ÞDW Nn :
ð1:3Þ
We refer to this approximation as semi-tamed Euler scheme. The increment function defined in [2] is /ðx; h; yÞ ¼ gðxÞ f ðxÞh þ 1þkgðxÞkh h þ rðxÞy. We will firstly show that the semi-tamed Euler scheme converges strongly to the exact solution
with the standard convergence order 0:5. Then we prove that the semi-tamed Euler scheme (1.3) can not only reproduce the exponential mean square stability of the exact solution, but also work better than the tamed Euler (1.2), backward Euler and split-step backward Euler schemes at stability preservation under a stepsize restriction. The rest of the paper is organized as follows. Section 2 begins with notations and estimation of the moment. Section 3 aims to show the strong convergence of the semi-tamed Euler scheme. Section 4 devotes to investigating the exponential mean square stability of the semi-tamed scheme for SDEs. Section 6 provides some numerical experiments for demonstration. 2. Estimation of the pth moments Throughout the whole article, unless otherwise specified, we use the following notations. Let T 2 ð0; 1Þ be a fixed real number and N 2 N be the step number of the uniform mesh with the stepsize T=N. Let ðX; F; PÞ be a complete probability space with a filtration fFt gtP0 satisfying the usual conditions and W t be a m-dimensional Brownian motion defined on this 1=2
probability space. Moreover, we use the notation kxk ¼ ðjx1 j2 þ þ jxk j2 Þ ; hx; yi ¼ x1 y1 þ x2 y2 þ þ xk yk for all x; y 2 Rk ; k 2 N, and kAk :¼ supx2Rl ; kxk61 kAxk for all A 2 Rkl ; k; l 2 N. a _ b represents maxfa; bg and a ^ b denotes minfa; bg. In this work, we also make the following assumptions. Assumption 2.1. Let gðxÞ be a continuously differentiable function and there exist positive constants K > 1 and c > 0, such that for any x; y 2 Rd
kf ðxÞ f ðyÞk _ krðxÞ rðyÞk 6 Kkx yk;
ð2:1Þ
hx y; lðxÞ lðyÞi 6 Kkx yk2 ;
ð2:2Þ
kg 0 ðxÞk 6 Kkxkc :
ð2:3Þ
Remark 2.1. There many stochastic models hold Assumption 2.1, such as stochastic Ginzburg–Landau equation, stochastic Verhulst equation, volatility processes (see [2]). Note that conditions (2.1) and (2.2) imply that SDE (1.1) admits the pth moment bounded solution in any finite time T (see [23]).
242
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
In order to prove the uniform boundedness of pth moments of the numerical solution (1.3), we follow the ideas in [12] to introduce the dominating stochastic processes DNn and appropriate sub-events XNn
DNn
! n1 X N 2 N ¼ ðk þ knkÞ exp k þ sup ½kkDW k k þ ak
ð2:4Þ
06u6n k¼u
and
( N n
) j sup DNk ð 06k6n1
x2X
X ¼
xÞ 6 N
2c
; sup kDW Nk k 06k6n1
61 ;
ð2:5Þ
where
k ¼ ð1 þ 3KT þ Tkf ð0Þk þ krð0Þk þ KÞ
4
and
aNn ¼ 1fkY Nn kP1g
* Y Nn þ f ðY Nn ÞT=N kY Nn k
;
rðY Nn Þ kY Nn k
+
DW Nn :
ð2:6Þ
The following lemmas are needed for proving the uniform boundedness of the pth moment. N N Lemma 2.1. Let Y N n ; Dn and Xn be given by (1.3), (2.4) and (2.5). Then
1XNn Y Nn 6 DNn
ð2:7Þ
for all n 2 f0; 1; 2; . . . ; Ng. Proof. Note that kDW Nn k 6 1 on XNnþ1 for all n 2 f0; 1; . . . ; N 1g and all N 2 N. Then Assumption 2.1 implies
kY Nnþ1 k 6 kY Nn k þ kf ðY Nn ÞkT=N þ kgðY Nn ÞkT=N þ krðY Nn ÞkkDW Nn k 6 1 þ Tkf ðY Nn Þ f ð0Þk þ Tkf ð0Þk þ TkgðY Nn Þk þ krðY Nn Þ rð0Þk þ krð0Þk 6 1 þ 3KT þ Tkf ð0Þk þ krð0Þk þ K 6 k N nþ1
ð2:8Þ
kY Nn k
on X \ fx 2 X : 6 1g for all n 2 f0; 1; 2; . . . ; Ng. Using the Cauchy–Schwarz inequality and the estimate ab 6 for all a; b 2 R, we have
kY Nnþ1 k2 ¼ kY Nn þ f ðY Nn ÞT=N þ
a2 2
2
þ b2
gðY Nn ÞT=N
þ rðY Nn ÞDW Nn k2 1 þ kgðY Nn ÞkT=N gðY N ÞT=N 2 2hY Nn ; gðY Nn ÞT=Ni N 2 N 2 2 2 N N 2 N N n ¼ kY n k þ kf ðY n Þk T =N þ þ krðY n ÞDW n k þ 2hY n ; f ðY n ÞT=Ni þ 1 þ kgðY Nn ÞkT=N 1 þ kgðY Nn ÞkT=N þ
2hf ðY Nn ÞT=N; gðY Nn ÞT=Ni 1þ
kgðY Nn ÞkT=N
þ 2hY Nn þ f ðY Nn ÞT=N; rðY Nn ÞDW Nn i þ
hgðY Nn ÞT=N; rðY Nn ÞDW Nn i 1 þ kgðY Nn ÞkT=N
6 kY Nn k2 þ 2kf ðY Nn Þk2 T 2 =N2 þ 3kgðY Nn ÞT=Nk2 þ 2krðY Nn ÞDW Nn k2 þ 2hY Nn ; f ðY Nn ÞT=Ni þ 2hY Nn þ f ðY Nn ÞT=N; rðY Nn ÞDW Nn i þ
2hY Nn ; lðY Nn ÞT=Ni 1 þ kgðY Nn ÞkT=N
2hY Nn ; f ðY Nn ÞT=Ni 1 þ kgðY Nn ÞkT=N
on X for all n 2 f0; 1; . . . ; N 1g and all N 2 N. The globally Lipschitz condition of f and condition of l give that 2
2
ð2:9Þ
:
r as well as the one side Lipschitz
2
kf ðxÞk2 6 ðkf ðxÞ f ð0Þk þ kf ð0ÞkÞ 6 ðKkxk þ kf ð0ÞkÞ 6 ðK þ kf ð0ÞkÞ kxk2 ;
ð2:10Þ
pffiffiffi jhx; f ðxÞij 6 kxkkf ðxÞ f ð0Þk þ kxkkf ð0Þk 6 Kkxk2 þ kf ð0Þkkxk 6 ðK þ kf ð0ÞkÞkxk2 6 kkxk2 ;
ð2:11Þ
krðxÞk2 6 ðkrðxÞ rð0Þk þ krð0ÞkÞ2 6 ðKkxk þ krð0ÞkÞ2 6 ðK þ krð0ÞkÞ2 kxk2 ; 6 kkxk2 ;
ð2:12Þ
pffiffiffi hx; lðxÞi 6 Kkxk2 þ kxkklð0Þk 6 kkxk2
ð2:13Þ
d
for all x 2 R with kxk P 1 and the polynomial bound on g 0 yields that
pffiffiffi kgðxÞk2 6 2K 2 kxk2ðcþ1Þ 6 N kkxk2
ð2:14Þ
243
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
for all x 2 Rd with 1 6 kxk 6 N 1=ð2cÞ and all N 2 N. Firstly, resorting to (2.11) and (2.13) as well as the definition of aNn , we have
2hY Nn þ f ðY Nn ÞT=N; rðY Nn ÞDW Nn i ¼ 2aNn kY Nn k2 on fx 2 X : 1 6
2hY Nn ; 1þ
kY Nn k
6N
ðY Nn ÞT=Ni kgðY Nn ÞkT=N
l
kgðY Nn ÞkT=N
ð2:15Þ
g for all N 2 N, and
6
2jhY Nn ; f ðY Nn ÞT=Nij 1þ
1=2c
6
pffiffiffi T 2 k kY Nn k2 ; N
ð2:16Þ
pffiffiffi T 2 k kY Nn k2 : N
2
Note that 2ðK þ kf ð0ÞkÞ T 2 < k and 3T 2 þ 6T 6 we obtain
ð2:17Þ pffiffiffi k. Then substituting the inequalities (2.10)–(2.12), (2.14)–(2.16) into (2.9),
pffiffiffi pffiffiffi 2 kY Nnþ1 k2 6 kY Nn k2 ð1 þ 2ðK þ kf ð0ÞkÞ T 2 =N þ 3 kT 2 =N þ 2kkDW Nn k2 þ 6 kT=N þ 2an Þ ! pffiffiffi 2 2ðK þ kf ð0ÞkÞ T 2 þ ð3T 2 þ 6TÞ k ¼ kY Nn k2 1 þ þ 2kkDW Nn k2 þ 2aNn N 2k þ 2kkDW Nn k2 þ 2aNn 6 kY Nn k2 exp N
ð2:18Þ
on fx 2 X : 1 6 kY Nn k 6 N 1=2c g for all N 2 N. Then we can establish (2.7) by the mathematical induction. Repeating the proof Lemma 2.1 in [12] completes this proof. h Lemma 2.2. For all p P 1,
sup N2N; NP4kpT
" !# N 1 X < 1: E exp pk kDW Nk k2
ð2:19Þ
k¼0
Proof. This proof is similar to Lemma 3.3 in [12] with only different k. h Lemma 2.3. Let an be given by (2.6). Then for all p P 1
" !# n1 X < 1: sup supE sup exp pz aNk
z2f1;1g N2N
06n6N
ð2:20Þ
k¼0
P N Proof. Note that the discrete stochastic process z n1 k¼0 ak is an ðFnT=N Þn2f0;1;...;Ng -martingale for every z 2 f1; 1g and every P N N 2 N. Then one can easily verify that the discrete stochastic process expðz n1 k¼0 ak Þ is a positive ðFnT=N Þn2f0;1;...;Ng submartingale for every z 2 f1; 1g and every N 2 N. Using Doob’s maximal inequality (see [24]), we have
! n1 X N sup exp z ak n2f0;1;...;Ng k¼0
Lp ðX;RÞ
! n1 X p N 6 exp z ak p1 k¼0
ð2:21Þ
Lp ðX;RÞ
for all N 2 N; p 2 ð1; 1Þ and all z 2 f1; 1g. Moreover, we have that
" 2 # x þ f ðxÞT=N rðxÞ p2 T krðxÞ ðx þ f ðxÞT=NÞk2 p2 T krðxÞk2 kx þ f ðxÞT=Nk2 E pz1fkxkP1g ; DW Nk ¼ 1fkxkP1g 6 1fkxkP1g 4 kxk kxk N N kxk kxk4 2
6
p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞ2 N
for all x 2 Rd ; k 2 f0; 1; . . . ; N 1g; N 2 N; p 2 ð1; 1Þ and all z 2 f1; 1g. The proof of Lemma 3.4 in [12] gives
!
2 x þ f ðxÞT=N rðxÞ p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞ2 ; 6 exp E exp pz1fkxkP1g ; DW Nk kxk kxk N which also shows that 2
E½expðpzaNk ÞjFkT=N 6 exp
p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞ2 N
!
244
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
for all x 2 Rd ; k 2 f0; 1; . . . ; N 1g; N 2 N; p 2 ½1; 1Þ and all z 2 f1; 1g. Hence, we have
" !# " ! # n1 n2 X X ¼ E exp pz aNk E½expðpzaNN1 ÞjFðN1ÞT=N E exp pz aNk k¼0
k¼0
"
n2 X 6 E exp pz aNk
!#
k¼0
2
p2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞÞ2 exp N
!
2
6 6 expðp2 Tð1 þ K þ kf ð0ÞkTÞ ðK þ krð0ÞkÞÞ2 Þ; which together with (2.21) implies
" !# n1 X <1 sup 2 f1; 1gsup E sup exp pz aNk z
N2N
06n6N
k¼0
d
for all x 2 R ; k 2 f0; 1; . . . ; N 1g; N 2 N; p 2 ð1; 1Þ and all z 2 f1; 1g. h The following two lemmas are based on Lemmas 2.2 and 2.3 and can be proved by following the proof Lemma 3.5 and Lemma 3.6 in [12]. Lemma 2.4. Let DN n be given by (2.4). Then for all p 2 ½1; 1Þ
sup N2N; NP8kpT
E sup jDNn jp < 1:
ð2:22Þ
06n6N
Lemma 2.5. Let XNN be given by (2.5). Then for all p 2 ½1; 1Þ c
supðNp P½ðXNN Þ Þ < 1:
ð2:23Þ
N2N
Now, we establish the moment boundedness theorem of the approximation (1.3). Theorem 2.6. Let Y N n be given by (1.3), then for all p 2 ð1; 1Þ
sup sup EkY Nn kp < 1:
ð2:24Þ
N2N 06n6N
Proof. Firstly, we rewrite (1.3) as the form
Y Nn ¼ n þ
n1 n1 X X T=N f ðY Nk Þ þ k¼0
k¼0
T=NgðY Nn Þ 1 þ kgðY Nn ÞkT=N
¼ n þ f ð0ÞnT=N þ rð0ÞW NnT=N þ
þ
n1 X
rðY Nk ÞDW Nk
k¼0
n1 n1 X X T=Nðf ðY Nk Þ f ð0ÞÞ þ k¼0
k¼0
T=NgðY Nn Þ 1þ
kgðY Nn ÞkT=N
þ
n1 X ðrðY Nk Þ rð0ÞÞDW Nk k¼0
for all n 2 f0; 1; . . . ; N 1g and all N 2 N. By the Burkholder–Davis–Gundy type inequality in [12, Lemma 3.8], we have
kY Nn kLp ðX;Rd Þ
X n1 N 6 knkLp ðX;Rd Þ þ kf ð0ÞknT=N þ kr þ T=N½f ðY k Þ f ð0Þ p d k¼0 L ðX;R Þ X X n1 n1 gðY Nk Þ þ T=N þ ½rðY Nk Þ rð0ÞDW Nk p d k¼0 1 þ kgðY Nk ÞkT=NLp ðX;Rd Þ k¼0 L ðX;R Þ !1=2 m n1 X X 2 N 6 knkLp ðX;Rd Þ þ kf ð0ÞknT=N þ p nT=N kri ð0Þk þ T=N½f ðY k Þ f ð0Þ p d i¼1 k¼0 L ðX;R Þ !1=2 n1 m XX þNþp T=Nkri ðY Nk Þ ri ð0Þk2Lp ðX;Rd Þ ; ð0ÞW NnT=N kLp ðX;Rd Þ
k¼0 i¼1
which together with the global Lipschitz condition of f and
r gives
n1 X pffiffiffiffiffiffiffi kY Nn kLp ðX;Rd Þ 6 2ðknkLp ðX;Rd Þ þ kf ð0ÞkT þ p Tmkrð0Þk þ NÞ þ ðK þ 2p2 mK 2 Þ T=NkY Nn kLp ðX;Rd Þ
!
k¼0
for all n 2 f0; 1; . . . ; N 1g and all N 2 N. Resorting to Gronwall’s lemma, we have that for all N 2 N and all p 2 ½2; 1Þ
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
sup n2f0;1;...;Ng
kY Nn kLp ðX;Rd Þ 6 2ðknkLp ðX;Rd Þ þ C 0 þ NÞ exp ðK þ 2p2 mK 2 ÞT ;
245
ð2:25Þ
pffiffiffiffiffiffiffi where C 0 ¼ kf ð0ÞkT þ p Tmkrð0Þk. Note that the inequality above does not prove this theorem due to N 2 N on the righthand side of the inequality (2.25). However, exploiting (2.25) in an appropriate bootstrap argument will enable us to establish Theorem 2.6. Firstly,
sup sup k1ðXN Þc Y Nn kLp ðX;Rd Þ 6 sup sup ðk1ðXN Þc kL2p ðX;Rd Þ kY Nn kL2p ðX;Rd Þ Þ n
N2N 06n6N
n
N2N 06n6N
6 sup Nk1ðXN Þc kL2p ðX;Rd Þ sup sup N 1 kY Nn kL2p ðX;Rd Þ n
N2N
N2N 06n6N
1=2p c 6 sup N2p P½ðXNn Þ sup sup N1 kY Nn kL2p ðX;Rd Þ < 1; N2N
ð2:26Þ
N2N 06n6N
where we used Hölder’s inequality, estimate (2.25) and Lemma 2.5. In addition, By Lemmas 2.1 and 2.4, we obtain
sup sup k1XNn Y Nn kLp ðX;Rd Þ 6 sup sup kDNn kLp ðX;Rd Þ < 1; N2N 06n6N
N2N 06n6N
which together with (2.26) implies the desired assertion.
h
3. Strong convergence of the semi-tamed Euler scheme In order to establish our convergence theorem for the semi-tamed Euler method (1.3), we now introduce appropriate time continuous interpolations of the time discrete numerical approximations (1.3). More accurately, we define the time continuous approximation yN : ½0; T X ! Rd ; N 2 N given by
t nT gðY Nn Þ nT N f ðY Nn Þ þ yNt ¼ Y Nn þ t þ rðY Nn ÞðW t W nT=N Þ N 1 þ kgðY Nn ÞkT=N
ð3:1Þ
for all t 2 ½nT ; ðnþ1ÞT Þ; n 2 f0; 1; . . . ; N 1g and all N 2 N. It is easy to see that yN is Ft -adapted stochastic process. N N Theorem 3.1. Let Assumption 2.1 hold. Then there exists a family C p 2 ð0; 1Þ; p 2 ½1; 1Þ, of real numbers such that
"
#!1=p
E sup kxt yNt kp
6 C p N1=2
ð3:2Þ
t2½0;T
for all N 2 N and all p 2 ½1; 1Þ. Before proving this theorem, we firstly present a lemma, which will be used in the proof of Theorem 3.1. Lemma 3.2. Let Y N n be defined by (1.3). Then we have
sup sup Ekf ðY Nn Þkp < 1; N2N n2f0;...;Ng
sup sup EkgðY Nn Þkp < 1 N2N n2f0;...;Ng
and
sup sup EkrðY Nn Þkp < 1 N2N n2f0;...;Ng
for all p 2 ½1; 1Þ. Proof. This proof can be completed by Assumption 2.1 and Theorem 2.6. h
kgðxÞkgðxÞ ðxÞ ¼ lðxÞ T=N 1þkgðxÞkT=N Proof of Theorem 3.1 Define yNs ¼ Y Nn for s 2 ½nT ; ðnþ1ÞT Þ for n 2 f0; 1; . . . ; N 1g; N 2 N; l and rewrite N N the continuous approximation (3.1) as
yNs ¼ n þ
Z
s
0
l ðyNu Þdu þ
Z
s 0
rðyNu ÞdW u :
Then we obtain from (1.1) and (3.3)
xs yNs ¼
Z 0
s
ðyNu Þdu þ ½ðlðxu Þ l
Z 0
s
½rðxu Þ rðyNu ÞdW u
ð3:3Þ
246
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
for all s 2 ½0; T and N 2 N. Itô’s formula produces
+ Z s* kgðyNu ÞkgðyNu Þ du xu yNu ; 1 þ kgðyNu ÞkT=N 0 0 m Z s m Z s X X i þ2 hxu yNu ; ri ðxu Þ ri ðyNu ÞidW u þ kri ðxu Þ ri ðyNu Þk2 du
kxs yNs k2 ¼ 2
Z
s
hxu yNu ; lðxu Þ lðyNu Þidu þ 2T=N
0
i¼1
¼2
Z
hxu yNu ; lðxu Þ lðyNu Þidu þ 2
0
þ
i¼1
s
m Z s X 0
i¼1
kri ðxu Þ ri ðyNu Þk2 du þ 2
Z 0
hxu yNu ; lðyNu Þ lðyNu Þidu þ 2T=N
m Z s X 0
i¼1
2
0
s
+ Z s* kgðyNu ÞkgðyNu Þ xu yNu ; du 1 þ kgðyNu ÞkT=N 0
i
hxu yNu ; ri ðxu Þ ri ðyNu ÞidW u :
2
By the inequality ða þ bÞ 6 2a2 þ 2b for all a; b 2 R, the Cauchy–Schwarz inequality and Assumption 2.1, we have
kxs yNs k2 6 ð2K þ 2K 2 mÞ þ 2K 2 m
Z 0
s
Z
s 0
Z
kxu yNu k2 du þ 2
kyNu yNu k2 du þ 2
m Z X i¼1
0
s
0
kxu yNu k klðyNu Þ lðyNu Þkdu þ 2T=N
s
Z 0
s
kxu yNu k kgðyNu Þk2 du
i
hxu yNu ; ri ðxu Þ ri ðyNu ÞidW u 2
for all s 2 ½0; T and all N 2 N. The inequality ab 6 a2 =2 þ b =2 for all a; b 2 R hence yields that
Z t t kxu yNu k2 du þ klðyNu Þ lðyNu Þk2 du þ T 2 =N2 kgðyNu Þk4 du 0 0 0 Z t X m Z t i 2 N N 2 N N þ 2K m kyu yu k du þ 2 sup hxu yu ; ri ðxu Þ ri ðyu ÞidW u ; 06s6t i¼1 0 0
sup kxs yNs k2 6 2ðK þ K 2 m þ 1Þ
06s6t
Z
Z
t
which together with the Burkholder–Davis–Gundy type inequality in [12, Lemma 3.7] implies
Z t Z t N 2 sup xs yN k2 k p=2 d 6 2ðK þ K 2 m þ 1Þ kx y k du þ klðyNu Þ lðyNu Þk2Lp ðX;Rd Þ du p d u s u L ðX;R Þ L ðX;R Þ 06s6t 0 0 Z t Z t þ T 2 =N2 kgðyNu Þk4L2p ðX;Rd Þ du þ 2K 2 m kyNu yNu k2Lp ðX;Rd Þ du 0
þp
m Z X i¼1
0
t 0
!1=2
khxu yNu ; ri ðxu Þ ri ðyNu Þik2Lp=2 ðX;RÞ du
ð3:4Þ
for all t 2 ½0; T; N 2 N and p 2 ½4; 1Þ. The last term in (3.4) has been estimated in the proof Theorem 1.1 of [12] as
p
m Z X i¼1
0
t
!1=2
khxu yNu ; ri ðxu Þ ri ðyNu Þik2Lp=2 ðX;RÞ du 2
6
2 Z 1 p2 K 2 m t sup kxs yN k þ kxs yNs k2Lp ðX;Rd Þ ds: s 2 06s6t 2 p d 0 L ðX;R Þ
ð3:5Þ
2
Using the inequality ða þ bÞ 6 2a2 þ 2b for all a; b 2 R and inserting (3.5) into (3.4) produce
2 sup kxs yN k s p 06s6t
L ðX;Rd Þ
!Z Z t t p2 K 2 m 62 K þK mþ1þ kxu yNu k2Lp ðX;Rd Þ du þ T 2 =N2 kgðyNu Þk4L2p ðX;Rd Þ du þ ð2K 2 m 2 0 0 2 Z t Z t 1 N sup kyNu yNu k2Lp ðX;Rd Þ du þ klðyNu Þ lðyNu Þk2Lp ðX;Rd Þ du þ kx y k ; þ p2 K 2 mÞ s s 2 06s6t 0 0 Lp ðX;Rd Þ 2
which implies
2 Z t Z t 1 2 2 2 2 N 2 sup kxs yN k 6 2ðK þ 1 þ p K mÞ kx y k du þ T =N kgðyNu Þk4L2p ðX;Rd Þ du p d u s u L ðX;R Þ 2 06s6t 0 0 Lp ðX;Rd Þ Z t Z t kyNu yNu k2Lp ðX;Rd Þ du þ klðyNu Þ lðyNu Þk2p þ ð2K 2 m þ p2 K 2 mÞ du L ðX;Rd Þ 0 0 for all t 2 ½0; T; N 2 N and p 2 ½4; 1Þ. The Gronwall lemma and the inequality a; b; c 2 ½0; 1Þ give
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi pffiffiffi pffiffiffi a þ b þ c 6 a þ b þ b for all
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
sup kxs yN k s 06s6T
6
pffiffiffiffiffiffi pffiffiffiffiffiffiffi 2T expð2K þ 2 þ 2p2 K 2 mÞ 2p2 K 2 2m sup kyNu yNu kLp ðX;Rd Þ 06u6T
Lp ðX;Rd Þ
247
þT=N sup kgðY Nn Þk2L2p ðX;Rd Þ n2f0;1;...;Ng
þ sup kl 06u6T
ðyNu Þ
ðyNu ÞkLp ðX;Rd Þ
l
! :
ð3:6Þ
Additionally, (3.3) and the Burkholder–Davis–Gundy type inequality show that for all N 2 N and p 2 ½2; 1Þ
sup 06u6T
kyNu
yNu kLp ðX;Rd Þ
Z u T T T N N N 6 sup kf ðyu ÞkLp ðX;Rd Þ þ sup kgðyu ÞkLp ðX;Rd Þ þ sup rðyu ÞdW u p d N 06u6T N 06u6T N 06u6T u L ðX;R Þ pffiffiffiffiffiffiffi T T p Tm N N sup sup kri ðY Nn ÞkLp ðX;Rd Þ ; 6 sup kf ðY n ÞkLp ðX;Rd Þ þ sup kgðY n ÞkLp ðX;Rd Þ þ pffiffiffiffi N 06n6N N 06n6N N i2f1;2;...;mg06n6N
which together with Lemma 3.2 gives
sup N2N
pffiffiffiffi N sup kyNu yNu kLp ðX;Rd Þ < 1
ð3:7Þ
06u6T
for all p 2 ½1; 1Þ. Particularly, we obtain
sup sup kyNu kLp ðX;Rd Þ < 1
ð3:8Þ
N2N 06u6T
for all p 2 ½1; 1Þ. Assumption 2.1 gives that
sup kyNu yNu kLp ðX;Rd Þ ; sup klðyNu Þ lðyNu ÞkLp ðX;Rd Þ 6 K 1 þ sup kyNu kcLp ðX;Rd Þ
06u6T
06u6T
06u6T
which together with (3.8) produces
sup N2N
pffiffiffiffi N sup klðyNu Þ lðyNu ÞkLp ðX;Rd Þ < 1
ð3:9Þ
06u6T
for all N 2 N and p 2 ½1; 1Þ. Combining (3.6), (3.7) and (3.9) and Lemma 3.2 gives (3.2). h 4. Exponential mean square stability of the semi-tamed Euler scheme In this section, we examine whether the semi-tamed scheme can share the exponential mean square stability of the exact solution. Since our aim involves the asymptotic behavior at the larger time, we rewrite the semi-tamed scheme (1.3) as follows:
Y nþ1 ¼ Y n þ f ðY n Þh þ
gðY n Þh þ rðY n ÞDW n ; 1 þ kgðY n Þkh
ð4:1Þ
where h > 0 is the stepsize and DW n ¼ W ðnþ1Þh W nh is the Brownian increment. For the stability purpose, we assume f ð0Þ ¼ gð0Þ ¼ 0 and rð0Þ ¼ 0, which shows that Eq. (1.1) admits a trivial solution. Firstly, we give the stability criterion of the exact solution (see [23]). Theorem 4.1. Let
l and r satisfy the local Lipschitz condition. If there exists a positive constant c such that
2hx; lðxÞi þ krðxÞk2 6 ckxk2 ;
8 x 2 Rd ;
ð4:2Þ
then the solution of Eq. (1.1) hold the property
Ekxt k2 6 Eknk2 ect :
ð4:3Þ
However, in order to examine stability of the numerical approximation Y n given by (4.1), we assume that there exist some with 2q > h2 and 2m > m such that nonnegative constants K; q; h; m; a > 1; m
hx y; f ðxÞ f ðyÞi 6 qkx yk2 ; hx y; gðxÞ gðyÞi 6 mkx ykaþ1 ;
kf ðxÞ f ðyÞk 6 Kkx yk; kxka kgðxÞk 6 m
ð4:4Þ ð4:5Þ
and
krðxÞ rðyÞk 6 hkx yk d
ð4:6Þ 2
for all x; y 2 R . Denote c ¼ 2q h > 0, then the three conditions above imply (4.3) by Theorem 4.1. Moreover, these conditions are also needed for the convergence.
248
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
Theorem 4.2. Let the conditions (4.4)–(4.6) hold. The there exist a stepsize h such that for any h < h , the discrete solution (4.1) holds the property
EkY n k2 6 Eknk2 enhch ;
ð4:7Þ
where ch > 0 satisfies limh!0 ch ¼ 2q h2 . Proof. The semi-tamed Euler approximation (4.1) gives 2
2
kY nþ1 k2 ¼ kY n k2 þ 2hY n ; f ðY n Þih þ kf ðY n Þk2 h þ krðY n ÞDW n k2 þ 2 þ 2 Y n þ f ðY n Þh þ
gðY n Þh ; rðY n ÞDW n : 1 þ kgðY n Þkh
hY n þ f ðY n Þh; gðY n Þhi kgðY n Þk2 h þ 2 1 þ kgðY n Þkh ð1 þ kgðY n ÞkhÞ ð4:8Þ
By the condition (4.5), we have 2
2
hkY n kaþ1 h hY n þ f ðY n Þh; gðY n Þhi mkY n kaþ1 h kf ðY n Þk kgðY n Þkh ½m K m 62 þ2 62 ; 1 þ kgðY n Þkh 1 þ kgðY n Þkh 1 þ kgðY n Þkh 1 þ kgðY n Þkh
ð4:9Þ
where K is the Lipschitz constant of f. Define Xn ¼ fx 2 X : kY n k > 1g, then we have 2
1Xn
kgðY n Þk2 h
ð1 þ kgðY n ÞkhÞ
2
6 1Xn
2
6 1Xcn
kgðY n Þkh mkY n kaþ1 h 6 1Xn 1 þ kgðY n Þkh 1 þ kgðY n Þkh
ð4:10Þ
and 2
1Xcn
kgðY n Þk2 h
ð1 þ kgðY n ÞkhÞ
m2 kY n k2a h2 1 þ kgðY n Þkh
6 1Xcn
m2 kY n kaþ1 h2 1 þ kgðY n Þkh
ð4:11Þ
:
mm ^ 2m . For any h < h , inserting (4.9)–(4.11) into (4.8) yields Define h1 ¼ 22K 1 Þm m ð2Kþm
2 kY nþ1 k2 6 kY n k2 þ 2hY n ; f ðY n Þih þ kf ðY n Þk2 h þ krðY n ÞDW n k2 þ 2 Y n þ f ðY n Þh þ
gðY n Þh ; rðY n ÞDW n : 1 þ kgðY n Þkh
ð4:12Þ
Taking expectations on the both sides of (4.12) gives that for any h < h1 , 2
EkY nþ1 k2 6 EkY n k2 þ 2E½hY n ; f ðY n Þih þ Ekf ðY n Þk2 h þ EkrðY n Þk2 h; where we used the independent increments of the Brownian motion W t . Hence, we obtain from (4.4) and (4.6) that 2
2
EkY nþ1 k2 6 EkY n k2 2qhEkY n k2 þ Kh EkY n k2 þ h2 hEkY n k2 ¼ ð1 2qh þ Kh þ h2 hÞEkY n k2 6 2
6 ð1 2qh þ Kh þ h2 hÞ
nþ1
Eknk2 6 Eknk2 ech ðnþ1Þh
2
for any h < h :¼ h1 ^ 2qKh . This completes the proof of Theorem 4.2. h
5. Numerical experiments This section presents some numerical experiments to illustrate our theoretical results. Let us firstly consider the scalar stochastic Ginzburg–Landau equation:
dxt ¼ ð2xt x3t Þdt þ xt dW t ;
t>0
ð5:1Þ
with the initial data x0 ¼ n ¼ 0:5. Then the drift lðxÞ ¼ 2x x3 do not satisfy the global Lipschitz condition, but can be decomposed into the Lipschitz continuous part f ðxÞ ¼ 2x and non-Lipschitz continuous part gðxÞ ¼ x3 . Fig. 1 depicts the root mean-square errors as a function of the stepsize h in log–log plot, where the expectation is approximated by the mean of 3 103 independent realizations. As expected, the semi-tamed Euler scheme gives an error that decreases proportional to 1=2 h . Then we investigate the numerical stability of the tamed scheme (1.3). We consider the scalar SDE:
dxt ¼ ð2xt x3t Þdt þ
pffiffiffi 2xt dW t ;
t>0
ð5:2Þ
It is easy to deduce from Theorem 4.1 that the solution (5.2) is exponentially mean square stable with the Lyapunov exponent not greater than 2, since
2xT lðxÞ þ jrðxÞj2 ¼ 4x2 2x4 þ 2x2 6 2x2 :
ð5:3Þ
249
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250 Convergence order
0
10
Semi−tamed Euler Slope of 0.5
Mean square error
−1
10
−2
10
−3
10
−3
−2
10
−1
10 Stepsize
10
Fig. 1. Mean square approximation error versus stepsize h to approximate (5.1).
Tamed Euler
Backward Euler
0
0
10
10
h = 1/4 h = 1/8 h = 1/16
2
E[Yn]
n
E[Y2]
h = 1/4 h = 1/8 h = 1/16
−10
−10
10
10
0
2
4
6
8
10
0
2
4
t
6
8
10
t
Split−step backward Euler
Semi−tamed Euler
0
0
10
10
h = 1/4 h = 1/8 h = 1/16
E[Yn]
2
n
E[Y2]
h = 1/4 h = 1/8 h = 1/16
−10
−10
10
10
0
2
4
6
t
8
10
0
2
4
6
8
10
t
Fig. 2. Computer simulation of Ejxt j2 , using different Euler schemes and different stepsizes.
We test the stability of the semi-tamed Euler, drift-tamed Euler, backward Euler and split-step backward Euler schemes with different stepsizes. We take h ¼ 1=4; 1=8 and 1=16 (h < h ¼ 0:4) and generate 7 103 sample paths for each numerical method. The mean square stability of the numerical solution is plotted in Fig. 2. Fig. 2 shows that the four numerical schemes are all mean square stable under the stepsizes h ¼ 1=4; 1=8 and 1=16, but the semi-tamed Euler scheme works better than the other three schemes. 6. Further remark The work proposes the semi-tamed Euler scheme to approximate the SDEs with non-Lipschitz continuous coefficients. We obtained its strong convergence order 0:5 and shown that the semi-tamed can reproduce the exponential mean square stability of the test equation. The numerical experiments also shown that the semi-tamed Euler scheme is superior to the tamed Euler and the backward Euler for reproducing mean square stability of the exact solution. We also meed to remark that (i) the conditions (2.3) and (4.5) on the non-Lipschitz continuous part g can be more complex, for example kg 0 ðxÞk 6 Kðkxka þ kxkb Þ; hx y; gðxÞ gðyÞi 6 t1 kx yk1þa t2 kx yk1þb and kgðxÞk 6 t1 kxka þ t2 kxkb ; (ii) the semi-tamed
250
X. Zong et al. / Applied Mathematics and Computation 228 (2014) 240–250
idea can also be applied to the classical Milstein scheme and lead to the semi-tamed Milstein scheme, which can be considered as the extension of the tamed Milstein scheme developed in [13]. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2013.11.100. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
L. Arnold, Random dynamical systems, Springer Monographs in Mathematics, Springer-Verlag, Berlin, 1998. M. Hutzenthaler, A. Jentzen, Convergence of the stochastic Euler scheme for locally Lipschitz coefficients, Found. Comput. Math. 11 (6) (2011) 657–706. B. Schmalfuss, The random attractor of the stochastic Lorenz system, Zeitschrift fur Ange-wandte Mathematik und Physik 48 (1997) 951–975. P.E. Kloeden, E. Platen, Numerical solution of stochastic differential equations, Applications of Mathematics (New York), vol. 2, Springer-Verlag, Berlin, 1992. M. Arato, A famous nonlinear stochastic equation (Lotka–Volterra model with diffusion), Math. Comput. Model. 38 (7) (2003) 709–726. G. Maruyama, Continuous Markov processes and stochastic equations, Rend. Circ. Mat. Palermo 4 (2) (1955) 48–90. G.N. Milstein, Numerical integration of stochastic differential equations, Mathematics and Its Applications, vol. 313, Kluwer Academic, Dordrecht, 1995. M. Hutzenthaler, A. Jentzen, Non-globally Lipschitz counterexamples for the stochastic Euler scheme. arXiv:0905.0273v1. M. Hutzenthaler, A. Jentzen, P.E. Kloeden, Strong and weak divergence in finite time of Eulers method for stochastic differential equations with nonglobally Lipschitz continuous coefficients, Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 467 (2011) 1563–1576. D.J. Higham, X. Mao, A.M. Stuart, Strong convergence of Euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40 (2) (2002) 1041–1063. X. Mao, L. Szpruch, Strong convergence and stability of implicit numerical methods for stochastic differential equations with non-globally Lipschitz continuous coefficients, J. Comput. Appl. Math. 238 (15) (2013) 14–28. M. Hutzenthaler, A. Jentzen, P.E. Kloeden, Strong convergence of an explicit numerical method for SDEs with nonglobally Lipschitz continuous coefficients, Ann. Appl. Probab. 22 (4) (2012) 1611–1641. X. Wang, S. Gan, The tamed Milstein method for commutative stochastic differential equations with non-globally Lipschitz continuous coefficients, J. Differ. Equ. Appl. 19 (3) (2013) 466–490. M. Sauer, W. Stannat, Lattice approximation for stochastic reaction diffusion equations with one-sided Lipschitz condition. arXiv:1301.6350. D.J. Higham, X. Mao, A. Stuart, Exponential mean-square stability of numerical solutions to stochastic differential equations, LMS J. Comput. Math. 6 (2003) 297–313. D.J. Higham, Mean-square and asymptotic stability of the stochastic theta method, SIAM J. Numer. Anal. 38 (2000) 753–769. C. Huang, Exponential mean square stability of numerical methods for systems of stochastic differential equations, J. Comput. Appl. Math. 236 (2012) 4016–4026. Y. Saito, T. Mitsui, Stability analysis of numerical schemes for stochastic differential equations, SIAM J. Numer. Anal. 33 (1996) 2254–2267. Y. Saito, T. Mitsui, Mean-square stability of numerical schemes for stochastic differential systems, Vietnam J. Math. 30 (2002) 551–560. D.J. Higham, X. Mao, C. Yuan, Almost sure and moment exponential stability in the numerial simulation of stochastic differential equations, SIAM J. Numer. Anal. 45 (2007) 592–609. X. Zong, F. Wu, Choice of h and mean-square exponential stability in the stochastic theta method of stochastic differential equations, J. Comput. Appl. Math. 255 (2014) 837–847. H. Szpruch, X. Mao, Strong convergence of numerical methods for nonlinear stochastic differential equations under monotone conditions, Preprint, Technical Report. X. Mao, Stochastic differential equations and their applications, Horwood Publishing Series in Mathematics & Applications, Horwood Publishing Limited, Chichester, 1997. A. Klenke, Probability Theory: A Comprehensive Course, Springer, London, 2008.