Theoretical and numerical analysis for Volterra integro-differential equations with Itô integral under polynomially growth conditions

Applied Mathematics and Computation 360 (2019) 70–82

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Theoretical and numerical analysis for Volterra integro-differential equations with Itô integral under polynomially growth conditions Huizi Yang a, Zhanwen Yang b,∗, Shufang Ma c a

Department of Statistic, School of Mathematics and Statistic Science, Ludong University, Yantai 264025, PR China Department of Mathematics, Harbin Institute of Technology, Harbin 150001, PR China c Department of Mathematics, North Forest University, Harbin 150040, PR China b

a r t i c l e

i n f o

MSC: 60H20 97N40 91B70 Keywords: Volterra integro-differential equations with Itô integral Semi-implicit Euler method Boundedness Hölder condition Strong convergence order

a b s t r a c t In this paper, we theoretically and numerically deal with nonlinear Volterra integrodifferential equations with Itô integral under a one-sided Lipschitz condition and polynomially growth conditions. It is proved that both the exact solutions and vector fields are bounded and satisfy a Hölder condition in the pth moment sense. Analogously, the boundedness and Hölder condition in the pth moment sense are preserved by the semi-implicit Euler method for sufficiently small step-size. Moreover, by the local truncated errors, we prove the strong convergence order 1. Finally, numerical simulations on stochastic control models and stochastic Ginzburg–Landau equation illustrate our results. © 2019 Elsevier Inc. All rights reserved.

1. Introduction In most problems, there are real phenomena depending on the effect of ‘white noise’ random forces in almost all areas of science, finance and engineering. These problems are intrinsically nonlinear and complex in nature and at least mathematically modeled and described by various generalized differential and integro-differential equations with the Itô type. Volterra integro-differential equations (VIDEs) with Itô integral play an key role in control theory and engineering. There are some results on theoretical analysis for VIDEs with Itô integral (see [22,23]) and many investigations on the controllability and stability (see e.g. [3,20,24]). Since exact solutions are rarely available, the numerical approximations become increasingly important in many applications. There has been a great deal of results on numerical approximations to ordinary differential equations (ODEs)

x (t ) = f (x(t )),

t ∈ I := [0, T ],

x ( 0 ) = x0 ,

where t ≥ 0 and the initial value x(0 ) = x0 (see [8,11,12,18,26]). For example, under a global Lipschitz condition of the nonlinear function f, there exists a unique exact solution in the interval I and the forward Euler method (FEM) and backward Euler method (BEM) are of order 1. Although the C1 -smooth condition of the nonlinear function f is not sufficient for the existence of the unique exact solution in the whole interval I, the FBM and BEM are still of order 1 whenever the exact

∗

Corresponding author. E-mail address: [email protected] (Z. Yang).

https://doi.org/10.1016/j.amc.2019.03.053 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

71

solution exists in the interval I, since the C1 -smooth condition ensures the boundedness and Hölder condition of the exact solution x(t) and the vector field f(x(t)) (see [26]). Some results are also extended to deterministic Volterra integro-differential equations (VIDEs)

dx(t ) = f (x(t )) + dt



t

0

σ (t, s )x(s )ds, t ∈ I, x(0 ) = x0 ,

(see [1,5,6,10]). Similarly to ODEs, the C1 -smooth conditions of the nonlinear function f and the kernel σ yield the boundedness and Hölder condition of the exact solution and the vector field existing in the whole interval I, which implies that the BEM is of order 1 (see the detailed discussion in Appendix A). Therefore, a basic task of this paper is to discuss the convergence order of numerical methods for nonlinear VIDEs with an Itô integral

dx(t ) = f (x(t )) + dt



0

t

σ (t, s )x(s )dB(s ), t ∈ I, x(0 ) = x0 .

(1)

under the non-global Lipschitz condition on f. With the development in the recent years, the strong convergence order of the Euler–Maruyama (EM) method, the semiimplicit Euler method, and so on, for stochastic ordinary differential equations (SODEs) are investigated in many literatures under a global Lipschitz condition (see [2,4,7,14,17,21,25,27]). However, there are few papers concerning on the strong convergence order for SODEs under a non-global Lipschitz condition. In [13], Higham et al. investigate the EM method and the semi-implicit Euler method under an assumption on the boundedness of both exact and numerical solutions in the pth moment sense and provide the strong convergence order of the semi-implicit Euler method with the condition that the nonlinear function f behaviors polynomially. Unfortunately, for SODEs with a superlinearly growing and one-sided Lipschitz continuous function f, the numerical solutions of the EM method are unbounded in the pth moment sense, though the pth moment boundedness of the exact solutions is ensured by a one-sided Lipschitz condition (see in [15]). Based on the boundedness of numerical solutions in the pth moment sense, they in [16] prove the strong convergence order of numerical methods if the nonlinear function f is one-sided Lipschitz continuous and has an at most polynomially growing derivative. In this paper, we consider the strong convergence order of the semi-implicit method for Eq. (1) under a one-sided Lipschitz condition and a polynomially growing condition. With some elementary inequalities in Section 2, the boundedness and Hölder condition of the exact solution and the vector field f are presented in Section 3. In Section 4, the semi-implicit Euler method for Eq. (1) is introduced and the numerical processes are provided by a one-step mapping. The one-sided Lipschitz continuous function f yields a global Lipschitz condition of the one-step mapping, which leads to the boundedness of numerical solutions in the pth moment sense. In Section 5, the strong convergence order is studied by the local truncated error based on a fundamental convergence analysis similarly to [2]. Compared to the results on SODEs, it is more interesting that the strong convergence order of the semi-implicit Euler method for VIDEs with an Itô integral is 1. In Section 6, some numerical experiments for stochastic control models and stochastic Ginzburg–Landau equations are presented to illustrate our theoretical results. 2. Preliminaries In this section, some notations and elementary inequalities on Itô intergal are introduced. In Eq. (1), the nonlinear function f : R → R is a C1 -smooth function, the kernel σ is a C 1 (D, R ) continuous function on D := {(t, s): 0 ≤ s ≤ t ≤ T}, the unknown solution x(t) is an Ft -adapted sampling-differentiable stochastic process in the complete probability space (, F, P ) with a filtration {Ft }t≥0 satisfying the usual conditions (i.e., it is right continuous and increasing while F0 contains all P-null sets) and B(t) is a standard Brownian motion with respect to the filtration {Ft }t≥0 . Let L p (, R ) be a family of Rvalued random variable X such that E|X|p < ∞ and M p ([0, T ]; R ) be the family of R-valued Ft -adapted processes such that T E 0 |g(t )| p dt < ∞. The following four properties are very useful (see [17,21]): (i) Young inequality

ab ≤

1 p 1 q a + b . p q

Here a and b are two positive constants, and (ii) Hölder inequality 1

1 p

+

1

E[(X T Y )] ≤ (E |X | p ) p (E |Y |q ) q , if p > 1, 1p + 1q = 1, X ∈ L p (, R ), Y ∈ Lq (, R ). (iii) Properties of the Paley–Wiener–Zygmund integral



b



1

(B(r ) − B(a ))dr 2 = (b − a )3 , 3 a  b  E (B(r ) − B(a ))dr = 0. E

a

1 q

= 1.

72

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

(iv) Burkholder–Davis–Gundy inequality

 t   t  p p 2 2 E | g ( s )d B ( s ) | ≤ C p E | | g( s ) | d s | , 0

where g ∈

0

M p ([0, T ]; R )

and Cp is a positive constant, which depends only on p.

3. Theoretical analysis of VIDEs with Itô integral From now on, we investigate the boundedness and Hölder condition of the exact solution x(t) and the vector field f(x(t)) of Eq. (1) in the pth moment sense, which are useful for the numerical analysis of the strong convergence order. Different from the deterministic problems, the boundedness of each sample path of exact solutions of Eq. (1) is variable with respect to the samples. Hence, the sample paths are not uniformly bounded. In Higham et al. [13], the strong convergence order is considered by the assumptions that the exact solutions of SODEs are momently bounded. Later on, it is shown that the moment boundedness is ensured by a one-sided Lipschitz condition (see in [16]). 3.1. Moment boundedness A one-sided Lipschitz condition on f is imposed in order to guarantee the moment boundedness of the exact solutions for Eq. (1). Assumption 3.1. Assume that the nonlinear function f satisfies a one-sided Lipschitz condition, i.e., there exists a constant ν such that

(x − y )( f (x ) − f (y )) ≤ ν|x − y|2 ,

(2)

for all x, y ∈ R. Remark 3.2. It is clear that the one-sided Lipschitz condition on the nonlinear function implies that there exists a positive constant Lp > 0 such that for all x ∈ R and even integer p,

x p−1 f (x ) ≤ L p (1 + |x| p ), which ensures the existence and uniqueness of the exact solution to Eq. (1) (see Theorem 3.6 of the Chapter 2 in [21]). Moreover, the one-sided Lipschitz condition also results to the pth-moment boundedness for any p ≥ 1. Theorem 3.3. Let the nonlinear function f satisfy a one-sided Lipschitz condition and E|x0 |p < ∞ for all p ≥ 2. Then there exists a constant H1 = H1 ( p) such that for all t ∈ I

E[|x(t )| p ] ≤ H1 ( p) < ∞. Proof. Without loss of generality, assume that p is an even integer. According to (1) and Assumption 3.1, we have

|x(t )| p + 1 = |x ( 0 )| p + 1 + ≤ |x ( 0 )| p + 1 +

 

t

0 t 0

px(s ) p−1 f (x(s ))ds + pL p (1 + |x(s )| p )ds +



t

0



0

t

px(s ) p−1 px(s ) p−1



s

0



0

s

σ (s, r )x(r )dB(r )ds σ (s, r )x(r )dB(r )ds.

(3)

Hence, by introducing the stopping times if necessary, and using the Fubini Theorem, the Young inequality and the Burkholder–Davis–Gundy inequality, we have for t ∈ I

E[|x(t )| p + 1] ≤ E[|x(0 )| p + 1] +



t

pL p E[(1 + |x(s )| p )]ds +

0

≤ E[|x(0 )| p + 1] + pL p ≤ E[|x(0 )| p + 1] + pL p  +

t 0



 

t 0 t 0

E[(1 + |x(s )| p )]ds +



t 0

t 0

 t 0



 

s 0

p

s

σ (s, r )x(r )dB(r )]ds  s p  p−1 1 p· E[|x(s )| p ] + p · E σ (s, r )x(r )dB(r ) ds

E[ px(s ) p−1

0

p

E[(1 + |x(s )| p )]ds

p−1 1 p· E[|x(s )| p ] + p · C p E p p

≤ E[|x(0 )| p + 1] +



σ 2 (s, r )x2 (r )dr

p

0

 2p

p [ pL p + ( p − 1 ) + T 2 C p σ ∞ ]E[|x(s )| p + 1]ds,

ds

(4)

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

73

where Cp is a given constant in Section 2 and σ ∞ = max{|σ (t, s )| : (t, s ) ∈ D}. Let G(t ) = sup0≤s≤t E[|x(s )| p + 1]. Then it follows from (4) and the Gronwall inequality that p G(t ) ≤ E[|x(0 )| p + 1] + [ pL p + ( p − 1 ) + T C p σ ∞ ]



t 0

G ( s )d s

p 2

p ≤ E[|x(0 )| p + 1]eT ( pL p +( p−1)+ σ ∞Cp T ) .

Therefore, the proof is completed.

(5)



Remark 3.4. For simplicity, we prove the moment boundedness of one-dimensional VIDEs with Itô integral. Whereas, the result also holds for multi-dimensional problems under the condition

< x − y, f (x ) − f (y ) >≤ ν x − y 2 for x, y ∈ Rd . Similarly to Theorem 4.1 in [21], the detailed proof is processed by replacing |x(t )| p + 1 into (1 + x(t ) 2 ) p/2 . For deterministic problems, the boundedness of the exact solutions is equivalent to that of the vector field f(x(t)), t ∈ I. For VIDEs with Itô integral, different from the global Lipschitz condition, the one-sided Lipschitz condition is not enough for the moment boundedness of the stochastic process f(x(t)). Hence, a polynomially growth condition on the nonlinear function f is introduced. Assumption 3.5. Assume that the nonlinear function f satisfies a polynomially growth condition, i.e., there exist a constant ν 0 > 0 and an integer k ≥ 1 such that

| f ( x ) | ≤ ν0 ( 1 + | x | k ) ,

(6)

for all x ∈ R. Corollary 3.6. Assume that Assumptions 3.1 and 3.5 hold and that E|x0 |p < ∞ for all p ≥ 1. Then there exists a constant H2 = H2 ( p) such that for all t ∈ I

E | f (x(t ))| p ≤ H2 ( p).

(7) 

Proof. The proof is straightforward and omitted. 3.2. Hölder continuous

Although the exact solution x(t) is sampling-differentiable, the Hölder condition is not trivial in the sense of the pth moment norm, while it is a corollary of the moment boundedness of the stochastic process f(x(t)). Theorem 3.7. Assume that the nonlinear function f satisfies Assumptions 3.1 and 3.5 and that E|x0 |p < ∞ for all p ≥ 1. Then the exact solutions satisfy the Hölder condition with exponent 1 with respect to any pth-moment norm. Namely, there exists a constant H3 = H3 ( p) such that for all 0 ≤ t1 ≤ t2 ≤ T

E |x(t1 ) − x(t2 )| p ≤ H3 (t2 − t1 ) p . Proof. According to Eq. (1), it is easily seen that

x(t2 ) − x(t1 ) =



t2

t1

f (x(s ))ds +



t2

t1



s 0

σ (s, r )x(r )dB(r )ds.

(8)

Following from the elementary inequality, Hölder inequality, Burkholder–Davis–Gundy inequality, Theorem 3.3 and Corollary 3.6, we get

E |x(t2 ) − x(t1 )| p  t  2 p−1 ≤2 E[| f (x(s ))ds| p ] + E[|

t2





s

σ (s, r )x(r )dB(r )ds| ]   t2  t2   s p−1 p−1 p p−1 p−1 p ≤ 2 (t2 − t1 ) E[| f (x(s ))| ]ds + 2 (t2 − t1 ) E | σ (s, r )x(r )dB(r )| ds. t1

t1

p

0

t1

t1

0

Therefore, with a similar progress to the proof of Theorem 3.3, ones obtain that

E |x(t2 ) − x(t1 )| p ≤ 2 p−1 H2 ( p)(t2 − t1 ) p + 2 p−1 (t2 − t1 ) p−1 T ≤ H3 (t2 − t1 ) p . This completes the proof.

p−2 2

p C p σ ∞



t2

t1



s 0

E[|x(r )| p ]drds (9)



74

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

In view of the derivation of the convergence for deterministic problems, the Hölder condition of f(x(t)) in the pth moment norm is more important than that of the exact solutions. Such a property is ensured by a stronger condition to Assumption 3.5 for SODEs in Hutzenthaler et al. [16]. Assumption 3.8. Assume that the derivative of the nonlinear function f satisfies a polynomially growth condition, there exist constants ν 2 > 0 and an integer k ≥ 1 such that

| f  ( x ) | ≤ ν2 ( 1 + | x | k ) ,

(10)

for all x ∈ R. Corollary 3.9. Let Assumptions 3.1 and 3.8 hold. Then there exists a constant H4 = H4 ( p) such that

E[| f (x(t )) − f (x(s ))| p ] ≤ H4 ( p)|t − s| p , for 0 ≤ s ≤ t ≤ T. Proof. Let 0 ≤ s ≤ t ≤ T. Then Assumption 3.8 ensures that there exists a random variable θ ∈ [0, 1] such that

| f (x(t )) − f (x(s ))| = | f  (x(s ) + θ (x(t ) − x(s )))||x(t ) − x(s )| ≤ ν2 |x(t ) − x(s )| + 2k−1 ν2 (|x(s )|k |x(t )) − x(s )| + |x(t ) − x(s )|k+1 ).

(11)

According to the Hölder inequality and Theorem 3.7, we have

E | f (x(t )) − f (x(s ))| p ≤ 3 p−1 [ν2p H3 ( p)|t − s| p + 2 p(k−1) ν2p (E[|x(t )|2kp ] ) 2 (E[|x(t ) − x(s )|2 p ] ) 2 + 2 p(k−1) ν2p H3 ( p(k + 1 ))|t − s| p(k+1) ] 1

1

≤ 3 p−1 [ν2p H3 ( p) + 2 p(k−1) ν2p (c (2 pk )) 2 (H (2 p)) 2 + 2 p(k−1) ν2p H3 ( p(k + 1 ))T pk ]|t − s| p . 1

This completes the proof.

1

(12)



4. Numerical methods for VIDEs with Itô integral In this section, the semi-implicit Euler method for VIDEs with Itô integral is introduced and the moment boundedness of the numerical solutions is also discussed. Under the same conditions, numerical solutions reproduce the moment boundedness. From now on, H stands for positive real constants dependent on p, k, T, but independent of h and its values may change between occurrences. 4.1. The semi-implicit Euler method The semi-implicit Euler method to Eq. (1) is given by

Yn+1 = Yn + h[ f (Yn+1 ) +

n

σ (tn+1 , tl )Yl Bl ],

(13)

l=0

where Y0 = x0 . Here h is a given step-size, tn = nh and Bn = B(tn+1 ) − B(tn ), n = 0, . . . , N − 1, is a random variable with Gaussian distribution N(0, h). Remark 4.1. The solvability of the semi-implicit Euler method is resulted from the contractive mapping principle when the nonlinear function f satisfies a global Lipschitz condition. While by the equivalent form

Yn+1 − h( f (Yn+1 ) − f (0 )) = Yn + h f (0 ) + h

n l=0

σ (tn+1 , tl )Yl Bl ,

(14)

the solvability is also available for hν < 1 under Assumption 3.1 (see [9]). Hence there always exists a unique numerical solution to (13) for all h > 0 when the one-sided Lipschitz constant ν is non-positive. Moreover, the numerical solution of (13) is represented by

Yn+1 = h (Yn + h f (0 ) + h

n

σ (tn+1 , tl )Yl Bl ),

l=0

where the one-step mapping h (x ) : R → R is implicitly defined by

h (x ) − h( f (h (x )) − f (0 )) = x. The moment boundedness of the numerical solutions is yielded from the global Lipschitz condition of the one-step mapping displayed the following lemma.

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

75

Lemma 4.2. Let p ≥ 1 be an integer and Assumption 3.1 hold. Then the one-step mapping satisfies a global Lipschitz condition, i.e., there exist two constants h∗ = h∗ ( p) > 0 and C = C ( p, ν ) such that for all x ∈ R



h (x )

2 p

≤ (1 + C ( p, ν )h )x2 p ,

whenever 0 < h < h∗ . Proof. In view of

x2 = (h (x ) − h( f (h (x )) − f (0 )))2 = h (x )2 − 2hh (x )( f (h (x )) − f (0 )) + h2 | f (h (x )) − f (0 )|2 ≥ h (x )2 − 2hh (x )( f (h (x )) − f (0 )) ≥ (1 − 2ν h )h (x )2 , the desired result is ensured by sufficiently small step-size.

(15)



4.2. Boundedness of numerical solutions It is ready to prove the moment boundedness of numerical solutions by induction. Theorem 4.3. Under the conditions of Theorem 3.3, there exists a constant H¯ 1 = H¯ 1 ( p) such that for 0 < h < h∗ (p)

E |Yn |2 p ≤ H¯ 1 ( p), where p ≥ 1 and n = 0, 1, . . . , N. Proof. By (14) and Lemma 4.2, we obtain that for 0 < h < h∗ and n = 1, . . . , N,





|Yn+1 |2 p = h Yn + h f (0 ) + h

n

2 p σ (tn+1 , tl )Yl Bl

l=0

≤ (1 + C ( p, ν )h )|Yn + h f (0 ) + h

n

σ (tn+1 , tl )Yl Bl |2 p .

(16)

l=0

Hence it follows from the definition of the convex function

|Yn+1 |2 p



≤ (1 + C ( p, ν )h ) Yn + h f (0 ) + h





= (1 + C ( p, ν )h ) (1 + h )



n−1

2 p σ (tn+1 , tl ) Bl Yl

l=0 n−1

Yn h + ( f (0 ) + σ (tn+1 , tl ) Bl Yl ) (1 + h ) (1 + h )

≤ (1 + C ( p, ν )h )(1 + h )2 p

l=0

n−1

1 h |Yn |2 p + ( f (0 ) + σ (tn+1 , tl ) Bl Yl )2 p (1 + h ) (1 + h ) l=0

 = (1 + C ( p, ν )h )(1 + h )

2 p−1

 2 p

|Yn | + h( f (0 ) + 2p

n−1

σ (tn+1 , tl ) Bl Yl )

l=0

≤ (1 + C ( p, ν )h )[(1 + (2 p − 2 )2 h ) + x(1 + (2 p − 2 )2 )h]



2p

 |Yn |2 p + h( f (0 ) +

n−1

l=0

≤ (1 + Hh )|Yn |2 p + Hh( f (0 ) +

n

l=0

which implies that

σ (tn+1 , tl ) Bl Yl )2 p ,



 σ (tn+1 , tl ) Bl Yl )2 p

76

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

E |Yn+1 |2 p

⎡ 2 p ⎤ n−1

≤ E |Yn |2 p (1 + Hh ) + HhE ⎣ f (0 ) + σ (tn+1 , tl ) Bl Yl ⎦ l=0

⎡ 2 p ⎤ n

≤ E |Yn |2 p (1 + Hh ) + Hh| f (0 )|2 p + HhE ⎣ σ (tn+1 , tl ) Bl Yl ⎦ l=0

≤ E |Yn |2 p (1 + Hh ) + Hh| f (0 )|2 p + Hh p+1 σ 2 p

n

E[(Yl )2 p ].

(17)

l=0

Similarly to the discussion on the exact solutions, we denote Gn = max{E |Ym |2 p : 0 ≤ m ≤ n}. Then

Gn+1 ≤ eHT E |Y0 |2 p + H | f (0 )|2 p (eHT − 1 ) + Hh p

n

Gl ,

l=1



The proof is complete by the discrete Gronwall inequality.

Remark 4.4. From the proof of Theorem 4.3, we obtain that the numerical solutions are bounded in the pth moment norm for all step-size h > 0 whenever ν ≤ 0. Moreover, the result in Theorem 4.3 also holds for multi-dimensional problems. In the end of this section, we present the moment boundedness of the numerical vector field in a similar way to the exact vector field. Corollary 4.5. Let p ≥ 1 be an integer and the conditions of Corollary 3.6 hold. Then there exist two constants h∗ = h∗ ( p) and H¯ 2 = H¯ 2 ( p) such that

E | f (Yn )| p ≤ H¯ 2 ( p), whenever 0 < h ≤ h∗ and n = 0, 1, . . . , N. Proof. The proof is straightforward and will be omitted.



5. Convergence analysis for VIDEs with Itô integral 5.1. Local truncated error In this section we investigate the convergence order by the approach of the fundamental convergence theorem similar to Baker and Buckwar [2]. Let x˜(tn+1 ) satisfy (13) when Yi = x(ti ) for i = 0, 1, . . . , n and define the local truncated error by

Rn+1 = x(tn+1 ) − x˜(tn+1 )  tn+1   = x(tn ) + f (x(s )) + tn

− x(tn ) − h f (x˜(tn+1 )) − h

s 0

 σ (s, r )x(r )dB(r ) ds

n

σ (tn+1 , tl )x(tl ) Bl .

(18)

l=0

Remark 5.1. It is easy to see from (18) that

 Rn+1 =

tn+1

tn

f (x(s ))ds − h f (x˜(tn+1 )) +



tn+1

tn



s 0

σ (s, r )x(r )dB(r )ds − h

n

σ (tn+1 , tl )x(tl ) Bl ,

(19)

l=0

of which an equivalent form is

Rn+1 − h( f (x(tn+1 )) − f (x(tn+1 ) − Rn+1 )) = Jn(1) + Jn(2) , where

Jn(1) = Jn(2) =



tn+1

tn



tn+1

tn

f (x(s ))ds − h f (x(tn+1 )), 

s 0

σ (s, r )x(r )dB(r )ds − h

n

σ (tn+1 , tl )x(tl ) Bl .

l=0

Remark 5.2. Similarly to the discussion in Section 4, the truncated error is represented by

Rn+1 = h (Jn(1) + Jn(2) , x(tn+1 )),

(20)

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

77

where the one-step mapping h (x, Z ) : R ×  → R is implicitly defined by

h (x, Z ) − h( f (h (x, Z ) + Z ) − f (Z )) = x, which satisfies a global Lipschitz condition with the constant 1 + C ( p, ν ) in Lemma 4.2. Lemma 5.3. If Assumptions 3.1 and 3.8 hold, then there exists a constant h∗ > 0 such that for h ∈ (0, h∗ ), (i) E |Jn(1 ) |2 p ≤ Hh4 p , (ii) E |Jn(2 ) |2 p ≤ Hh3 p . Proof. By Hölder inequality, we have

|Jn(1) |2 p ≤ h2 p−1



tn+1

| f (x(s )) − f (x(tn+1 ))|2 p ds.

tn

(21)

By Corollary 3.9, we obtain

E |Jn(1) |2 p ≤ Hh4 p .

(22) Jn2 ,

According to the definition of

Jn(2) =

 n−1 

l=0

+

tn+1

tn



tl+1

tl

it is rewritten as the form

(σ (s, r )x(r ) − σ (tn+1 , r )x(r )

 σ (tn+1 , r )x(r ) − σ (tn+1 , tl )x(r ) + σ (tn+1 , tl )x(r ) − σ (tn+1 , tl )x(tl ))dB(r ) ds.

(23)

Hence taking the expectation, we obtain from the elementary inequality, Hölder inequality and Assumption 3.8 that

E |Jn(2) |2 p ≤ n2 p−1

n−1

E|



tn+1



tn

l=0

 (2hL1 |x(r )| + σ ∞ |x(r ) − x(tl )| )dB(r ) ds|2 p

tl+1

tl

  n−1

≤ 22 p−1 n2 p−1 E|

tn+1

tn

l=0

|



tl+1

tl

2hL1 |x(r )|dB(r )ds|2 p + E |



tn+1

tn

 tl

tl+1

 σ ∞ |x(r ) − x(tl )| )dB(r )ds|2 p .

Then we obtain according to the similar progress to the proof of Theorem 3.3 that

E |Jn(2) |2 p ≤ 22 p−1 n2 p−1 h2 p−1  +

  tn+1 E |

tn

tl

 n−1 

l=0

tl+1

tn+1

E

tn

  |

tl

tl+1

 2hL1 |x(r )|dB(r )|2 p ds

σ ∞ |x(r ) − x(tl )| )dB(r )|

  2p

ds

≤ Hh3 p . This completes the proof.



The estimations of the local truncated errors and the numerical vector field are presented based on Lemma 5.3. Theorem 5.4. Under the conditions of Corollary 3.9, it holds for 0 ≤ h ≤ 4ν1 p ,

E |Rn+1 |2 p ≤ Hh3 p for n = 0, . . . , N − 1. Proof. This proof is similar to the proof of Theorem 4.3. Using the formula (20) and Lemma 4.2, we have

|h (Jn(1) + Jn(2) , x(tn+1 ))|2 p ≤ (1 + C ( p, ν )h )|Jn(1) + Jn(2) |2 p , Hence the proof is completed by Lemma 5.3.

(24)



Corollary 5.5. Under the conditions of Corollary 3.9, there exists a constant h∗ > 0 such that for h ∈ (0, h∗ ), p ≥ 1 and n = 0, 1, . . . , N − 1,

E[| f (x˜n+1 ) − f (x(tn+1 ))|2 p ] ≤ Hh3 p .

(25)

Finally, we consider the mean error of Rn+1 analogously to the fundamental convergence theorem of SODEs in [2]. Corollary 5.6. Under the conditions of Corollary 3.9, there exists a constant h∗ such that for 0 < h < h∗ and n = 0, 1, . . . , N − 1

78

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

E[Rn+1 ] ≤ Hh2 . Proof. According to the definition (18), Hölder inequality, Assumption 3.8, Theorem 5.4, and Corollaries 3.9 and 5.5, we have

E[Rn+1 ]  ≤E

tn+1

tn

 ( f (x(s )) − f (x˜(tn+1 )))ds



 n−1 

+E

 ≤ hE

l=0 tn+1

tn

tn+1



tn

tl+1

tl

 (2hL1 |x(r )| + |σ (tl , r )x(r ) − σ (tn+1 , tl )x(tl ))dB(r ) ds



| f (x(s )) − f (x˜(tn+1 ))|ds

≤ Hh2 .

(26)

This completes the proof.



5.2. Convergence order Let εn = x(tn ) − Yn denote the global error between the exact solution and the numerical solution. Then the 1.0 strong convergence order of the semi-implicit Euler method for Eq. (1) is presented in our main theorem. Theorem 5.7. Let Assumptions 3.1 and 3.8 hold. Then there exists a constant h∗ such that for 0 < h ≤ h∗ ,

max E |εn |2 ≤ Hh2 .

n=0,...,N

Proof. By the definitions, the global error and the local truncated error have the relationship

εn+1 = x(tn+1 ) − Yn+1 = x(tn+1 ) − x˜n+1 + x˜n+1 − Yn+1 = Rn+1 + εn + h( f (x˜n+1 ) − f (x(tn+1 ))) + h( f (x(tn+1 )) − f (Yn+1 )) + δn+1 , where

 δn+1 = h

n

(27)

 σ (tn+1 , tl )εl Bl .

(28)

l=0

By Assumption 3.1, it holds that



εn+1 − h( f (x(tn+1 )) − f (Yn+1 ))

2

≥ |εn+1 |2 − 2hεn+1 ( f (x(tn+1 )) − f (Yn+1 ))) ≥ |εn+1 |2 − 2ν h|εn+1 |2 = (1 − 2ν h )|εn+1 |2 .

(29)

Hence, we obtain

(1 − 2ν h )|εn+1 |2

≤ Rn+1 + εn + h( f (x˜(tn+1 )) − f (x(tn+1 ))) + δn+1 2 ≤ |εn |2 + 3|Rn+1 |2 + 3h2 | f (x˜(tn+1 )) − f (x(tn+1 ))|2 + 3|δn+1 |2 + 2εn Rn+1 + 2hεn ( f (x˜(tn+1 )) − f (x(tn+1 ))) + 2εn δn+1 .

(30)

According to Theorems 3.3 and 4.3, we have

E |δn+1 |2 ≤ h2 T σ 2∞ E

n−1 

l=0

tn+1

tn

|x(tl ) − Yl |2 dr.

(31)

Therefore, the strong convergence order is verified by Theorem 5.4 and Corollaries 5.5 and 5.6 in a similar to Theorem 3 in [2] and the detailed proof is omitted.  Remark 5.8. It is clear see that the at most 2k + 2th moment boundedness of Rn+1 is required when we estimate the meansquare moment of the global error ε n .

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

79

Table 1 Strong convergence order of the semi-implicit Euler method for Example 6.1. Stepsizes

2−5 2−6 2−7 2−8

σ (t, s ) = 1

σ (t, s ) = cos(t − s )

σ (t, s ) = ts

Error

Order

Error

Order

Error

Order

0.0061 0.0031 0.0015 0.0 0 07

– 0.9714 1.0 0 07 0.9657

0.0066 0.0033 0.0016 0.0 0 07

– 0.9887 1.0053 0.9776

0.0032 0.0015 0.0 0 07 0.0 0 03

– 1.0136 0.9983 1.0291

Table 2 Strong convergence order of the semi-implicit Euler method for Example 6.2. Stepsizes

2−5 2−6 2−7 2−8

σ (t, s ) = 1

σ (t, s ) = cos(t − s )

σ (t, s ) = ts

Error

Order

Error

Order

Error

Order

0.0055 0.0028 0.0014 0.0 0 07

– 0.9608 1.0264 0.9770

0.0049 0.0025 0.0012 0.0 0 06

– 0.9891 0.9944 0.9740

0.0042 0.0021 0.0011 0.0 0 05

– 0.9724 0.9895 1.0039

6. Numerical simulation In this section, several examples are presented for illustrating the foregoing convergence results for Eq. (1) with the three kinds of kernels : (1) σ (t, s ) = 1 is a simplest kernel, (2) σ (t, s ) = cos(t − s ) is a convolution kernel, (3) σ (t, s ) = ts is a general kernel. We choose 40 0 0 sample paths to simulate the global errors between the exact and numerical solutions and the strong convergence orders with four different stepsizes h = 2−5 , 2−6 , 2−7 , 2−8 by Matlab The strong convergence error is defined by

h = E |x(T ) − YN |2 and the strong convergence order is numerically computed by

Order = log( h / h )/ log(2 ). 2

Example 6.1. In the first example, we consider the following problem

dx(t ) = − sin x3 (t ) + dt



t 0

σ (t, s )x(s )dB(s ), t ∈ [0, 1], x(0 ) = 1,

(32)

which satisfies a global Lipschitz condition, which is used to model stochastic control systems in [20,24]. In Table 1, we list the strong convergence errors and orders. It illustrates that the strong convergence order 1 is independent of the kernels, which is different from the previous related results for stochastic Volterra integral equations (see in [19]). Example 6.2. In order to verify our strong convergence analysis for VIDEs with a non-global Lipschitz condition, we consider a dissipative system with an Itô integral

dx(t ) = −x(t ) − x3 (t ) + dt



t

0

σ (t, s )x(s )dB(s ), t ∈ [0, 1], x(0 ) = 1,

(33)

of which the nonlinear function satisfies a one-sided Lipschitz condition and a polynomially growth condition. The digits in Table 2 indicate that the strong convergence order of the semi-implicit Euler method is 1, which is agreement with our theoretical analysis. Example 6.3. We apply the semi-implicit Euler method to an integral version of stochastic Ginzburg–Landau equation in [17], i.e.,

dx(t ) = x(t ) − x3 (t ) + dt



0

t

σ (t, s )x(s )dB(s ), t ∈ [0, 1], x(0 ) = 1.

(34)

From the digits in Table 3, the strong convergence order 1 is available for VIDEs with a positive one-sided Lipschitz constant, which is agreement with our theoretical analysis. 7. Future works Due to a lot of real and practical problems modeled by Volterra integral equations with non-smooth and weakly singular kernels, we will investigate the strong convergence order of the semi-implicit Euler method for a kind of VIDEs with nonsmooth and weakly singular kernels. To avoid the errors from quadratures, we are more interested with the semi-implicit

80

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82 Table 3 Strong convergence order of the semi-implicit Euler method for Example 6.3. Stepsizes

σ (t, s ) = 1 Error

2−5 2−6 2−7 2−8

0.0056 0.0028 0.0014 0.0 0 07

σ (t, s ) = cos(t − s )

σ (t, s ) = ts

Order

Error

Order

Error

Order

– 0.9886 1.0041 0.9586

0.0063 0.0031 0.0016 0.0 0 07

– 1.0444 0.9623 1.0038

0.0036 0.0018 0.0 0 09 0.0 0 04

– 0.9790 1.0019 1.0084

Table 4 Strong convergence order of the semi-implicit Euler method with distribution. Stepsizes

2−5 2−6 2−7 2−8

σ (t, s ) = (t − s ) 2

σ (t, s ) = (t − s )− 4

Error

Order

Error

Order

0.0013 0.0 0 06 0.0 0 03 0.0 0 01

– 1.0 0 04 1.0055 1.0027

0.0188 0.0108 0.0063 0.0036

– 0.7995 0.7715 0.7997

1

1

Euler method with Gaussian distributions, i.e.,



Y¯n+1 = Y¯n + h f (Y¯n+1 ) +

n



ωn,l Y¯l ,

(35)

l=0

t where ωn,l is a random variable with Gaussian distribution N (0, t l+1 |σ (tn+1 , s )|2 ds ). In Table 4, we list the strong converl gence errors and orders of the method Eq. (35) to Eq. (34) with a non-smooth kernel

dx(t ) = x(t ) − x3 (t ) + dt



t 0

(t − s )α x(s )dB(s ), t ∈ [0, 1],

x ( 0 ) = 1,

where α = 12 , − 14 . Table 4 illustrates the strong convergence order of the semi-implicit Euler method with distribution is still 1 for non-smooth continuous kernels. Moreover, the strong convergence order for VIDEs with weakly singular kernels varies with α ∈ (−1, 0 ). In the sequel papers, we will analyze the strong convergence order theoretically. Acknowledgment This research was supported by the National Natural Science Foundation of China (NSFC 1187011327 and NSFC 11771128). The first author is also supported by the Natural Science Foundation of Shandong Province, China (No. ZR2019PA024) and an MOE (Ministry of Education in China) Youth Foundation Project of Humanities and Social Sciences (Project No. 19YJC630180). Appendix A. Volterra integro-differential equations We consider deterministic nonlinear Volterra integro-differential equations (VIDEs) of the form

x (t ) = f (x(t )) +



t 0

σ (t, s )x(s )ds, t ∈ I,

(A.1)

with initial value x(0 ) = x0 , where f : R → R is a nonlinear function, the kernel σ is continuous on D and x is the unknown solution. It is well-known that there exists a unique solution to Eq. (A.1) in I when the nonlinear function f satisfies a global Lipschitz condition and then the present numerical analysis illustrates that the convergence orders of the EM, the BEM and the semi-implicit Euler method are 1. While we are more interested with the convergence analysis under a C1 -smooth condition of the nonlinear function f. Since the C1 -smooth condition yields a local Lipschitz condition on f, there always exists a unique local solution to Eq. (A.1) and without loss of generality we assume that the local solution is continuable to the interval I. The following propositions are useful for the numerical analysis and their proofs are trivial. Proposition A.1. Let f ∈ C 1 (R, R ) and x(t) be the unique continuous solution to Eq. (A.1) in I. Then x(t) and f(x(t)) are bounded. Proposition A.2. Let f ∈ C 1 (R, R ) and x(t) be the unique continuous solution to Eq. (A.1) in I. Then there exists a constant Q0 > 0 such that it holds for all 0 ≤ s, t ≤ T,

|x(t ) − x(s )| ∨ | f (x(t )) − f (x(s ))| ≤ Q0 |t − s|.

(A.2)

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

81

Proposition A.3. Assume that σ ∈ C1 (D) and x(t) is the unique continuous solution to Eq. (A.1) in I. Then there exists a constant Q1 > 0 such that it holds for all 0 ≤ s, t ≤ T,

 t   s    σ (t, r )x(r )dr − σ (s, r )x(r )dr ≤ Q1 |t − s|.  0  0

(A.3)

Applying the semi-implicit Euler method to VIDEs, one obtains the following numerical scheme

Yn+1 = Yn + h[ f (Yn+1 ) + h

n

σ (tn+1 , tl )Yl ],

(A.4)

l=0

where h is the step-size, tn = nh and Yn is an approximation to x(tn ). For the convergence analysis, we only need to show that the numerical solutions are bounded for sufficiently small step-size, since the C1 -smooth condition of the nonlinear function f plays the same role as the global Lipschitz condition on bounded domains. Instead of the numerical solutions, we analyze the boundedness of the errors n = x(tn ) − Yn , which satisfies the following difference equation

n+1 = n + hFn ( n+1 ) + h2

n

σ (tn+1 , tl ) l + Mn (h ),

(A.5)

l=0

where Fn (e ) = f (x(tn+1 )) − f (x(tn+1 ) − e ) and

Mn ( h ) =



tn+1

tn

f (x(s )) − f (x(tn+1 ))ds +

n 

+h

l=0

tl+1

tl



tn+1

tn



s 0

σ (s, r )x(r )dr −



tn+1 0

 σ (tn+1 , r )x(r )dr ds

(σ (tn+1 , r )x(r ) − σ (tn+1 , tl )x(tl ) )dr.

For convenience, we introduce some notations:

M = max |x(t )| ∨ max | f (x(t ))|, t∈I

t∈I

L = max | f  (x )| ∨ max |x|≤2M

(t,s )∈D

|σ (t, s )|.

Remark A.4. Since Mn (h) is independent of the numerical solutions, it follows from Propositions A.2 and A.3 that for f ∈ C 1 (R, R ) and σ ∈ C1 (D), there exists a constant > 0 such that

|Mn ( h )| ≤ h2 , for all 0 ≤ n ≤ N − 1. Remark A.5. The C1 -smooth condition of the nonlinear function f yields that for all 0 ≤ n ≤ N − 1 and |e| ≤ M

|Fn (e )| ≤ L|e|. It is ready to show the main theorem in this section. Theorem A.6. Assume that f ∈ C 1 (R, R ), σ ∈ C1 (D) and x(t) be the unique continuous solution to Eq. (A.1) in I. Then there exists an h∗ > 0 such that for all 0 < h < h∗ , the errors of the semi-implicit Euler method satisfy

| n | ≤ M

e(2L+1)nh for 0 ≤ n ≤ N. e(2L+1)T

Proof. Let h∗ ≤ min

1

,



M

2L 2 e(2L+1 )T

(A.6)

be such that for all 0 < h < h∗ ,

(1 + h )e−(2L+1)h + Lh < 1   e(2L+1 )nh and denote Bn = x ∈ R : |x| ≤ M . (2L+1 )T

e In view of 0 = 0, (A.6) holds for n = 0. By induction, we suppose that (A.6) holds up to n ≤ N − 1 and define a mapping

by

T (e ) = n + hFn (e ) + h2

n

σ (tn+1 , tl ) l + Mn (h ).

l=0

By 0 < h < h∗ , we obtain from (A.6) and Remarks A.4 and A.5 that for e ∈ Bn+1 ,

82

H. Yang, Z. Yang and S. Ma / Applied Mathematics and Computation 360 (2019) 70–82

|T (e )| ≤ | n | + Lh|e| + h2 L

n

| l | + h2

l=0

e(2L+1)lh e(2L+1)nh e(2L+1)(n+1)h hM 2 ≤ M (2L+1)T + LhM + h L M (2L+1)T + (2L+1)T e e(2L+1)T e 2e n

M ≤ (2L+1)T e

l=0



e(2L+1)nh − 1 1 e(2L+1)nh + Lhe(2L+1)(n+1)h + h2 L (2L+1)h + h 2 e −1

M ≤ (2L+1)T (1 + h )e(2L+1)nh + Lhe(2L+1)(n+1)h e

e(2L+1)(n+1)h e(2L+1)(n+1)h ≤M (1 + h )e−(2L+1)h + Lh ≤ M (2L+1)T . e(2L+1)T e



(A.7)

Together with

|T (e ) − T (e¯ )|

≤ Lh|e − e¯| ≤

1 2

|e − e¯|

for e, e¯ ∈ Bn+1 , T : Bn+1 → Bn+1 is a contractive mapping. Thus the error n+1 is the fixed point n+1 = T ( n+1 ) ∈ Bn+1 satisfying (A.6). Hence the proof is completed.  Therefore, the convergence theorem is formulated in the following. Theorem A.7. Assume that f ∈ C 1 (R, R ), σ ∈ C1 (D) and x(t) be the unique continuous solution to Eq. (A.1) in I. Then the semiimplicit Euler method for Eq. (A.1) is of convergence order 1. References [1] K.E. Atkinson, The Numerical Solution of Integral Equations of the Second Kind, Cambridge University Press, Cambridge, 1997. [2] C.T.H.B. Baker, E. Buckwar, Exponential stability in p-th mean of solutions, and of convergent euler-type solutions to stochastic delay differential equations, J. Comput. Appl. Math. 184 (2005) 404–427. [3] K. Balachandran, S. Karthikeyan, Controllability of nonlinear Ito type stochastic integro-differential systems, J. Frankl. Inst. 345 (2008) 382–391. [4] A.F. Bastani, M. Tahmasebi, Strong convergence of split-step backward euler method for stochastic differential equations with non-smooth drift, J. Comput. Appl. Math. 236 (2012) 1903–1918. [5] H. Brunner, A survey of recent advances in the numerical treatment of Volterra integral and integro-differential equations, J. Comput. Appl. Math. 8 (1982) 213–229. [6] H. Brunner, Collocation Methods for Volterra Integral and Related Function Equations, Cambridge Unversity Press, Cambridge, 2004. [7] E. Buckwar, C. Kelly, Non-normal drift structures and linear stability analysis of numerical methods for systems of stochastic differential equations, Comput. Math. Appl. 64 (2012) 2282–2293. [8] J.C. Butcher, Numerical Methods for Ordinary Differential Equations, Wiley, Chichester, 2008. [9] K. Dekker, J.G. Verwer, Stability of Runge–Kutta Methods for Stiff Nonlinear Differential Equations, North-Holland, Amsterdam, 1984. [10] E. Hairer, C. Lubich, S.P. Norsett, Order of convergence of one-step methods for Volterra integral equation of the second kind, SIAM J. Numer. Anal. 20 (1983) 569–579. [11] E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, Springer-Verlag, Berlin, 1987. [12] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, Springer-Verlag, Berlin, 1991. [13] D.J. Higham, X.R. Mao, A.M. Stuart, Strong convergence of euler-type methods for nonlinear stochastic differential equations, SIAM J. Numer. Anal. 40 (2002) 1041–1063. [14] C.M. Huang, Exponential mean square stability of numerical methods for systems of stochastic differential equations, J. Comput. Appl. Math. 236 (2012) 4016–4026. [15] M. Hutzenthaler, A. Jendtzen, P.E. Kloeden, Strong and weak divergence in finite time of Euler’s method for stochastic differential equations with non-globally Lipschitz continuous coefficients, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 467 (2010) 1563–1576. [16] M. Hutzenthaler, A. Jendtzen, P.E. Kloeden, Strong convergence of an explicit numerical method for SDEs with non-globally Lipschitz continuous coefficients, Ann. Appl. Probab. 22 (2012) 1611–1641. [17] P.E. Kloden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer-Verlog, Berlin, 1992. [18] J.D. Lambert, Numerical Methods for Ordinary Differential Systems, Wiley, Chichester, 1991. [19] H. Liang, Z.W. Yang, J.F. Gao, Strong super convergence of the Euler–Maruyama method for linear stochastic Volterra integral equations, J. Comput. Appl. Math. 317 (2017) 447–457. [20] B. Liu, Stability of solutions for stochastic impulsive systems via comparison approach, IEEE Trans. Autom. Control 53 (2008) 2128–2133. [21] X.R. Mao, Stochastic Differential Equations and Applications, Horwood Pub Limited, 1997. [22] X.R. Mao, Stability of stochastic intergro-differential equations, Stoch. Anal. Appl. 18 (20 0 0) 10 05–1017. [23] X.R. Mao, M. Riedle, Mean square stability of stochastic intergro-differential equations, Syst. Control Lett. 55 (2006) 459–465. [24] L.J. Shen, J.P. Shi, J.T. Sun, Complete controllability of impulsive stochastic integro-differential systems, Automatica 46 (2010) 1068–1073. [25] D.E. Stewart, Rigid-body dynamics with friction and impact, SIAM Rev. 42 (20 0 0) 3–39. [26] A.M. Stuart, A.R. Humphries, Dynamical Systems and Numerical Analysis, Cambridge University Press, 1998. [27] F. Wu, X. Mao, P. Kloeden, Discrete Razumikhin-type technique and stability of the Euler–Maruyama method to stochastic functional differential equations, Discrete Contin. Dyn. Syst. 33 (2013) 885–903.