Chebyshev cardinal wavelets for nonlinear stochastic differential equations driven with variable-order fractional Brownian motion

Chaos, Solitons and Fractals 124 (2019) 105–124

Contents lists available at ScienceDirect

Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos

Frontiers

Chebyshev cardinal wavelets for nonlinear stochastic differential equations driven with variable-order fractional Brownian motion M.H. Heydari a, Z. Avazzadeh b,∗, M.R. Mahmoudi c a

Department of Mathematics, Shiraz University of Technology, Shiraz, Iran School of Mathematical Sciences, Nanjing Normal University, Nanjing 210023, China c Department of Statistics, Faculty of Science, Fasa University, Fasa, Iran b

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 4 November 2018 Revised 13 March 2019 Accepted 29 April 2019 Available online 13 May 2019

This paper is concerned with a computational approach based on the Chebyshev cardinal wavelets for a novel class of nonlinear stochastic differential equations characterized by the presence of variableorder fractional Brownian motion. More precisely, in the proposed approach, the solution of a nonlinear stochastic differential equation is approximated by the Chebyshev cardinal wavelets and subsequently the intended problem is transformed to a system of nonlinear algebraic equations. In this way, the nonlinear terms are significantly reduced, due to the cardinal property of the basis functions used. The convergence analysis of the expressed method is theoretically investigated. Moreover, the reliability and applicability of the approach are experimentally examined through the numerical examples. In addition, the presented method is implemented for some famous stochastic models, such as stochastic logistic problem, stochastic population growth model, stochastic Lotka–Volterra problem, stochastic Brusselator problem, stochastic Duffing-Van der Pol oscillator problem and stochastic pendulum model. As another new finding, a procedure is established for constructing the variable-order fractional Brownian motion. Indeed, the standard Brownian motion together with the block pulse functions and the hat functions are utilized for generating the variable-order fractional Brownian motion.

Keywords: Stochastic differential equations (SDEs) Chebyshev cardinal wavelets (CCWs) Variable-order fractional Brownian motion (V-Ofbm) Stochastic Lotka–Volterra system Stochastic Brusselator problem Stochastic Duffing-Van der Pol oscillator problem Stochastic pendulum model

© 2019 Published by Elsevier Ltd.

1. Introduction In this study, we present a newly established method based on the Chebyshev cardinal wavelets (CCWs) for the solution of the following class of nonlinear stochastic differential equations (SDEs):



dX (t ) = f (t, X (t ) )dt + g(t, X (t ) )dBH (t ) (t ),

t ∈ [0, 1],

X (0 ) = X0 , (1.1) where X(t) is a stochastic process in a certain probability space (, F, P ) and X0 is a deterministic initial value. Also, f and g are given smooth functions. Here, BH(t) (t) indicates a variable-order fractional Brownian motion (V-OfBm) of the Riemann–Liouville type with H(t) ∈ (0, 1) defined by Zeng et al. [1]

BH (t ) (t ) =



1

 H (t ) +

 1 2



t 0

(t − s )H (t )− 2 dB(s ),

∗

1

(1.2)

Corresponding author. E-mail addresses: [email protected] (M.H. Heydari), [email protected] (Z. Avazzadeh), [email protected] (M.R. Mahmoudi). https://doi.org/10.1016/j.chaos.2019.04.040 0960-0779/© 2019 Published by Elsevier Ltd.

where B(t) is the standard Brownian motion (Bm). Note that the classical fractional Brownian motion (fBm) is a particular case of the V-OfBm where H (t ) = H and also it is a standard Bm if H (t ) = 1 2. We remind that fractional calculus (derivatives and integrals of arbitrary orders) has been widely used to model the behavior of various phenomena in physics and engineering. For instance, this concept has been used in fluid mechanics [2], electromagnetism [3], dynamics of viscoelastic materials [4,5] and propagation of spherical flames [6]. In recent years, many researches have been developed on numerical solution of the problems involving fractional derivative and integral operators, e.g. [7–12]. Also, some modifications have been performed on the structure of these operators, for instance see [13–19]. As known, SDEs driven with the classical fBm have appeared in modeling of many phenomena in economic data [20], biology [21], turbulence [22] and medicine [23]. Moreover, the fBm is getting popular among scholars since the definition of classical fBm by using the standard Bm [24–30]. In the past few years, approximate methods for solution of the stochastic problems driven with the standard Bm and the classical fBm have drawn attention of many researches [31–46].

106

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

It is worth noting that the classical fBm can only model monofractal phenomena with the same global irregularity and constant ¨ memory as characterized by the constant Holder exponent H. In practice, global self-similarity seldom exists and fixed scaling is only held for a certain finite range of intervals. Furthermore, due to empirical data, the scaling exponent and the order of selfsimilarity usually have more than one value. Hence, in complex heterogeneous systems there commonly exist phenomena which represent multi-fractal properties with variable space and time dependent memory [47,48]. The V-OfBm is an alternative way to tackle these limitations as it has appeared in modelling of many phenomena containing variable irregularities or variable memory such as modelling of network traffic and signal processing [49,50], geophysics for terrain modelling [51], financial time series for stochastic volatility modelling [52] and anomalous diffusion with variable memory [53,54]. We emphasize that obtaining the exact solutions of stochastic problems driven with V-OfBm is mostly impossible. Hence, researchers mainly focus on computational approaches to obtain approximate solutions. For solving the SDEs expressed in Eq. (1.1), we need to construct the V-OfBm for obtaining numerical results. To do this, two families of the basis functions, namely the block pulse functions (BPFs), and the hat functions (HFs) are exploited. This procedure uses the standard Bm, the operational matrix (OM) of derivative of the BPFs and the OM of variable-order (V-O) fractional integration of the HFs. In the sequel, we convert the SDE expressed in Eq. (1.1) to a corresponding stochastic integral equation. Then, we expand the solution of the integral equation in terms of the CCWs. Precisely, by using both of the integration and differentiation operational matrices of the CCWs and applying the technique of fast computing nonlinear terms, we achieve a system of algebraic equations containing the coefficients of CCWs. Technically, this system is far easier than Eq. (1.1) and can even be solved by a few commands in mathematical softwares such as Maple, Mathematica or MATLAB. We remind that the CCWs are privileged due to their cardinality, orthogonality and exponential accuracy. Moreover, these wavelets inherit the properties of wavelets and the Chebyshev cardinal polynomials simultaneously. Hence, this set of basis functions is a very powerful tool in approximation theory. The convergence analysis of the approach is theoretically investigated. In addition, the accuracy of the method is experimentally examined by solving some test problems. The contents of the paper arefollowed as: in Section 2, some requisite definitions are reviewed e.g. V-O fractional calculus, the BPFs and the HFs. Constructing the V-OfBm is described in Section 3. In Section 4, the CCWs and their properties are analyzed. In Section 5, the description of method is presented. The convergence analysis is investigated in Section 6. The results of several numerical examples are demonstrated in Section 7. The approach is applied for solving some well-known models in Section 8. Finally, some conclusions are highlighted in Section 9. 2. Basic definitions Here, we present some basic definitions which are used in the following sections. 2.1. Variable-order (V-O) fractional calculus Many definitions are proposed since the theory of V-O fractional calculus is appeared. Here, we only hint at a definition of V-O fractional integral which is used in the sequel. Definition 2.1 [55]. The V-O fractional integral operator of order α (t) ≥ 0 of a function u(t) in the Riemann–Liouville type is given

by



 α (t ) 

u (t ) =

I



1

(α (t )) f (t ),

t

0

(t − s )α (t )−1 u(s )ds,

α (t ) > 0, (2.1) α (t ) = 0,

where  is the Gamma function and α (t) is a known continuous function. Remark 1. The below property is an immediate result of the above definition

(β + 1 ) t α (t )+β . (α (t ) + β + 1 )

Iα (t ) t β =

(2.2)

2.2. The block pulse functions (BPFs) The BPFs are defined on [0, T] by Kilicman et al. [56]



bi (t ) =

1,

i−1 T ˆ m

0,

otherwise,

≤t ≤

i T, ˆ m

(2.3)

ˆ . Any absolutely integrable function u(t) defined for i = 1, 2, . . . , m ˆ terms of the BPFs as over [0, T] can be approximated by m

u(t )  C T Bmˆ , where

ci =

CT

(2.4)

= [c1 c2 . . . cmˆ ] and Bmˆ = [b1 b2 . . . bmˆ

]T ,

and in which

   iT  iT ˆ ˆ m m ˆ ˆ m m ( 2i − 1 )T u(t )bi (t )dt = u(t )dt  u . T ( i−1 T ( i−1 2m T T ˆ ) ˆ ) m m (2.5)

Integration of the vector Bmˆ for q-times can be represented as follows



t



···

t

 0

0

Bmˆ (s )(ds )q  Pˆ(q ) Bmˆ (t ).

(2.6)

q−times

ˆ ×m ˆ matrix namely the OM of q-times integrawhere Pˆ(q ) is an m tion of the BPFs given by Kilicman et al. [56]

⎛

Pˆ(q ) =

T q ˆ m

1 ⎜0 ⎜0 ⎜

1 ( q + 1 )! ⎜

⎝0 0

λ1 1 0 0 0

λ2 λ1 1 0 0

... ... ... .. . 0

⎞ λmˆ −1 λmˆ −2 ⎟ λmˆ −3 ⎟ ⎟, .. ⎟ . ⎠

(2.7)

1

and λ j = ( j + 1 )q+1 − 2 jq+1 + ( j − 1 )q+1 . Moreover, it is shown in [56] that

dq Bmˆ (t )  Dˆ (q ) Bmˆ (t ), dt q

(2.8)

ˆ ×m ˆ matrix namely the OM of q-times differenwhere Dˆ (q ) is an m tiation of the BPFs and has the following form

⎛

1  q ⎜0 ⎜0 ˆ m ⎜ Dˆ (q ) = (q + 1 )! ⎜ T ⎝0 0 with ϑ j = −

j i=1

ϑ1 1 0 0 0

ϑ2 ϑ1 1 0 0

... ... ... .. . 0

⎞ ϑmˆ −1 ϑmˆ −2 ⎟ ϑmˆ −3 ⎟ ⎟, .. ⎟ . ⎠

(2.9)

1

λi ϑ j−i for all j = 1, 2, . . . , mˆ − 1 and ϑ0 = 1.

2.3. The hat functions (HFs) ˆ -set of the HFs is defined over [0, T] as follows [31,32] An m

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

 ϕ0 (t ) =

h−t , h 0,

0 ≤ t < h,

(2.10)

0,

bi (t ) 

(i − 1 )h ≤ t < ih, ih ≤ t < (i + 1 )h,

ϕmˆ −1 (t ) =

T − h ≤ t ≤ T,

(2.12)

otherwise,

ˆ − 1 ). Any function u(t ) ∈ C ([0, T ] ) can be approxwhere h = T /(m imated by the HFs as

u(t ) 

ˆ −1 m 

u(ih )ϕi (t )  F T mˆ (t ),

(2.13)

i=0

F = [ u ( 0 ) u ( h ) u ( 2 h ) . . . u ( T )] T ,

(2.14)

mˆ (t ) = [ϕ0 (t ) ϕ1 (t ) . . . ϕmˆ −1 (t )]T .

(2.15)

Theorem 2.2 [57]. Assume mˆ (t ) be the vector introduced in Eq. (2.15) and α (t ) : [0, T ] −→ R+ be a positive continuous function. The V-O fractional integration of order α (t) of mˆ (t ) in the Riemann– Liouville type can be represented as

(Iα (t ) mˆ )(t )  P¯ (α ) mˆ (t ),

(2.16)

where P¯ (α ) is the OM of V-O fractional integration of order α (t) for the HFs as

0 ⎜0 ⎜0 ⎜ P¯ (α ) = ⎜ . ⎜ .. ⎝ 0 0

ˆ m 

bi (( j − 1 )h )ϕ j−1 (t ) = QiT mˆ (t ),

j=1

where Qi is the ith row and Qij is the (i, j)th component of the matrix Q.  3. Constructiong the varaible-order fractional Brownian motion (V-OfBm) In this section, we introduce a new algorithm using the BPFs and HFs to construct the V-OfBm defined in Eq. (1.2) in the type of Riemann-Liouville. To do this, we present a two-step strategy where in the first step, the standard Bm is constructed by using the definition of this stochastic process and employing the spline interpolation. In the next step, the standard Bm is expanded in terms of the BPFs, and then the relation between these basis functions and the HFs is employed to achieve the V-OfBm. Step 1: Constructing the standard Brownian motion (Bm)

where

⎛

bi ( jh )ϕ j (t ) =

(2.19) ˆ − 2, i = 1, 2, . . . , m

otherwise,

t − (T − h ) , h 0,

ˆ −1 m  j=0

(2.11)



ˆ be the ith component of Proof. Assume bi (t) for i = 1, 2, . . . , m Bmˆ (t ). Approximating bi (t) by the HFs yields

otherwise,

⎧ t − ( i − 1 )h ⎪ ⎪ , ⎪ ⎨ h ϕi (t ) = (i + 1 )h − t , ⎪ ⎪ h ⎪ ⎩

107

ζ1 θ1 1 0 .. . 0 0

ζ2 θ1 2 θ2 2 .. . 0 0

... ... ... .. . 0 0

ζmˆ −2 θ1 mˆ −2 θ2 mˆ −2 .. .

θmˆ −2 mˆ −2 0

⎞ ζmˆ −1 θ1 mˆ −1 ⎟ θ2 mˆ −1 ⎟ ⎟ ⎟ .. ⎟ . ⎠

θmˆ −2 mˆ −1 θmˆ −1 mˆ −1

,

ˆ ˆ ×m m

where

  hα ( jh ) ( j − 1 )α ( jh)+1 + jα ( jh) (α ( jh ) − j + 1 ) , (α ( jh ) + 2 ) ˆ − 1, j = 1, 2, . . . , m

ζj = and

⎧ 0, ⎪ ⎪ ⎪ ⎪ ⎪ hα ( jh ) ⎪ ⎨ , θi j = (α (αjh( jh) )+ 2 ) ⎪ h ⎪ ⎪ ( j − i + 1 )α ( jh)+1 − 2( j − i )α ( jh)+1 ⎪ ⎪ (α ( jh ) + 2 ) ⎪  ⎩ + ( j − i − 1 )α ( jh )+1 ,

j < i, j = i,

j > i.

Lemma 2.3. Suppose that Bmˆ (t ) and mˆ (t ) are the BPFs and HFs vectors introduced in Eqs. (2.4) and (2.15), respectively. The vector Bmˆ (t ) can be approximated as

Bmˆ (t )  Q mˆ (t ),

Algorithm 1 The Maple commands for the standard Bm. > restart: > Digits := 15: > with(LinearAlgebra): > with(stats, random): > with(CurveFitting): > N := 50: > T := 2: > t := NT : > Times := convert([seq(i. t, i = 0.. N)], Array): > B := convert([seq(0, i = 0.. N)], Array): > for i from 2 to N+1 do √ B[i] := B[i-1]+ random[normald[0, t ]](1): end do: > B := unapply(Spline(Times, B, t, degree = 1), t): > plot(B(t), t = 0.. T, title = “The standard Bm”, labels = [“t”, “B(t)”]);

Step 2: Computing the V-O fractional integral of the standard Bm

(2.17) The standard Bm generated in the previous step can be expanded in terms of BPFs as

ˆ ×m ˆ matrix namely the BPFs matrix and where Q is an m

Qi j = bi (( j − 1 )h ),

We remind that the standard Bm B(t), t ∈ [0, T] provides the below properties [33]: (I) B(0 ) = 0, a.s. √ (II) B(t ) − B(s ) ∼ t − s N (0, 1 ), where 0 ≤ s < t ≤ T and N (0, 1 ) is normally distributed with mean 0 and variance 1. (III) The increments B(t ) − B(s ) and B(v ) − B(u ) for 0 ≤ s < t < u < v ≤ T are independent. Note that for numerical simulations, one can employ the discretized Bm. In fact, we can evaluate B(t) at some distinguished points ti and use an appropriate interpolation method for constructing B(t). In this case, we define t = T /N for N ∈ Z+ and Bi = B(ti ) where ti = i t and i = 0, 1, . . . , N. The first property means that B0 = 0 with the probability 1; and (II) and (III) imply that Bi = Bi−1 + dBi for i = 1, 2, . . . , N, such that each term of dBi is an √ independent random variable in the form of t − sN (0, 1 ). The procedure of forming the standard Bm on the interval [0,2] with N = 50 is presented in Algorithm 1.

ˆ. i, j = 1, 2, . . . , m

(2.18)

B(t )  C T Bmˆ (t ),

(3.1)

108

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

)T where C T = [c1 c2 . . . cmˆ ] and ci = B (2i2−1 ˆ m



ˆ, for i = 1, 2, . . . , m

and Bmˆ (t ) is introduced in Eq. (2.4). Replacing Eq. (3.1) into Eq. (1.2) yields

B

H (t )

(t ) 



1

   H (t ) + 12

t

0

(t − s )

H (t )− 12





d C Bmˆ (s ) . T

B

(t ) 



1

   H (t ) + 12

t

0

(t − s )

H (t )− 12





C Dˆ (1) Bmˆ (s ) ds, T

an by considering Lemma 2.3 and Theorem 2.2, we have

BH (t ) (t ) 

1

   H (t ) + 12 T

¯ (α )

A P

⎛



t 0

(3.3)

1⎜ ⎟ (t − s )H (t )− 2 ⎝C T Dˆ (1) Q mˆ (s )⎠ds 

mˆ (t )  Z mˆ (t ),



> b:= (i, t ) −→ piecewise

( 2i − 1 )T



(3.4)









,

,

2k

(4.1)

otherwise,

.

(4.2)



1

0

⎧ ⎨ ⎩

ψnm (ξ )ψn m (ξ )wn (ξ )dξ π

M2k+1 0,

1



(n, m ) = (n , m ), (n, m ) = (n , m ),

,

1 − 2k+1 ξ − 2n + 1

2

,

0,

ˆ ) , Matrix : > C := convert seq(c (i ), i = 1.m > λ := (i, j ) −→ piecewise i > j, 0, i = j, 1, i < j, ( j − i + 1 )

2k



2M

ψnm (ξ ), ψn m (ξ )wn =

wn ( ξ ) =



n − 1 n 

We emphasize that the set {ψnm (ξ )| n = 1, 2, . . . , 2k , m = 1, 2, . . . , M} provides an orthogonal basis for L2wn [0, 1]. Moreover,

where

iT ≤t ≤ ,1 . ˆ m

.

ξ∈

(4.3)

n − 1 n  2k

,

2k

, (4.4)

otherwise.

Thus, any function u(ξ ) ∈ L2wn [0, 1] can be approximated in terms of the CCWs as

2

− 2 ( j − i )2 + ( j − i − 1 )2 : ˆ , λ): > Pˆ(1 ) = T Matrix(m

k

ˆ 2m

> Dˆ (1 ) := MatrixInverse(Pˆ(1 ) ): > σ := (i, j ) −→ b(i, ( j − 1 )h ). ˆ , σ ): > Q := Matrix(m   > ϕ := (i, t ) −→ piecewise i = 0, piecewise 0 ≤ t ≤ ˆ 1≤i≤m − 2, piecewise (i − 1 )h ≤ t ≤ ih,

( 2m − 1 )π

=

ˆ m



ˆ 2m



ηm = − cos



( i − 1 )T

ξ∈

where ξnm = k+1 (ηm + 2n − 1 ), k ∈ Z+ ∪ {0}, n = 1, 2, . . . , 2k , and 2 ηm for m = 1, 2, . . . , M imply the roots of the Mth order Chebyshev polynomial [58] on [−1, 1], i.e.

AT T

H:= t −→ 0.5 + 0.2 sin(50t ). α := t −→H(t)+ 12 . ˆ := 400: m h := mˆ 1−1 :

> c:= i −→ B

⎧   M  ⎪ − ξn ξ ⎪ ⎪ , ⎨ ξnm − ξn ψnm (ξ ) =  = 1 ⎪ ⎪ = m ⎪ ⎩ 1

Algorithm 2 The Maple commands for the second step of constructing the V-OfBm.



The CCWs are defined on [0,1] by Heydari et al. [46]

0,

⎞

where α (t ) = H (t ) + 12 . The Maple commands for the second step of constructing the V-OfBm over the interval [0,2] when H (t ) = ˆ = 400 is presented in Algorithm 2. 0.5 + 0.2 sin(50t ) and m

> > > >

4. The Chebyshev cardinal wavelets (CCWs) and their properties

(3.2)

From Eq. (2.8), we achieve H (t )

terval [0,2]. The behavior of the obtained numerical solutions are illustrated in Fig. 1.

t−(i−1 )h , h

ˆ − 1, piecewise (T − 1 ) ≤ t ≤ T , i=m

u (ξ ) 

unm ψnm (ξ )  T (ξ ),

(4.5)

n=1 m=1

h−t h

where



,



ih ≤ t ≤ (i + 1 )h, (i+1h)h−t , t−(T −h ) h

2  M 

 = [u11 u12 . . . u1M |u21 u22 . . . u2M | . . . |u2k 1 u2k 2 . . . u2k M ]T , (ξ ) = [ψ11 (ξ ) ψ12 (ξ ) . . . ψ1M (ξ )|ψ21 (ξ ) ψ22 (ξ ) . . . ψ2M (ξ )|



:

> mˆ := t −→Transpose(convert([seq( ϕ (i, t ), i = 0.mˆ − 1)], Matrix)):  hα (( j−1 )h ) > ζ := (i, j, α ) −→ piecewise i < j, (α (( j − 1 )h ) + 2 )



 ( j − 2 )α (( j−1)h)+1 + ( j − 1 )α (( j−1)h) (α (( j − 1 )h ) − j + 2 ) :  hα (( j−1 )h ) > θ := (i, j, α ) −→ piecewise i = j, , (α (( j − 1 )h ) + 2 )  α (( j−1 ) h ) i < j, (αh(( j−1 )h )+2 ) ( j − i + 1 )α (( j−1 )h )+1  −2( j − i )α (( j−1 )h )+1 + ( j − i − 1 )α (( j−1 )h )+1 : > aa := (i, j ) −→ piecewise(i = 1, ζ (i, j, α ), θ (i, j, α ) ): ˆ , aa ): > P¯ (α ) := Matrix(m > A := Multiply(C, Multiply(Dˆ (1 ) , Q )): > Z := Multiply(A, P¯ (α ) ): > BH := t −→ Multiply(Z, mˆ (t ) )[1, 1]:

> plot(BH(t), t = 0. T, title = “The V-OfBm”, labels = [“t”, “BH(t)”]);

ˆ = 400 We have used the suggested scheme with N = 50 and m to construct the V-OfBm with some choices of H(t) over the in-

|ψ2k 1 (ξ ) ψ2k 2 (ξ ) . . . ψ2k M (ξ )]T ,

(4.6)

n = 1, 2, . . . , 2k , m = 1, 2, . . . , M.

(4.7)

... and

unm = u(ξnm ),

Lemma 4.1 [46]. Let G : [0, 1] × R −→ R be a continuous function, ! T (ξ ) are the CCWs expansions of u(ξ ) and the and T  (ξ ) and  identity function, respectively. Then, we get





! T , T (ξ ), G (ξ , u (ξ ) )  G  where





! T , T = [G (u˜11 , u11 ) G (u˜12 , u12 ) . . . G (u˜1M , u1M ) G 

|G (u˜21 , u21 ) G (u˜22 , u22 ) . . . G (u˜2M , u2M )| . . . |G (u˜2k 1 , u2k 1 ) G (u˜2k 2 , u2k 2 ) . . . G (u˜2k M , u2k M )]T .

(4.8)

Theorem 4.2 [46]. The differentiation of the vector  (ξ ) introduced in Eq. (4.6) can be represented as

d(ξ ) = D(1) (ξ ), dξ

(4.9)

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

109

ˆ = 400 for some choices of H(t). Fig. 1. The behavior of the V-OfBm where N = 50 and m

where D(1) is a 2k M × 2k M matrix namely the OM of differentiation of the CCWs as

⎛

Q 0 ⎜ D (1 ) = ⎜ . ⎝ .. 0

⎞

... ... .. . 0

0 Q .. . 0

0 0⎟ .. ⎟ ⎠, . Q

(4.10)

qi j =

⎪ 2k+1 ⎪ ⎪ ⎪ ⎪ η i − ηj ⎪ ⎪ =1 ⎪ ⎩ M 

i = j,

η − η   j , ηi − η

where the CCWs is as follows

⎛

L ⎜0

E L

P ( 1 ) = ⎜0 ⎜. ⎝ .. 0

0 .. . 0

⎜ ⎜

matrix

E E .. . .. . ...

... ... .. . L 0

P(1)

⎛ ei j =

1 2k+1

1

 = i

0

g(s, X (s ) )dBH (s ) (s ),

(5.1)

(5.2)

 = [x11 x12 . . . x1M |x21 x22 . . . x2M | . . . |x2k 1 x2k 2 . . . x2k M ]T . (5.3) Substituting Eq. (5.2) into Eq. (5.1) gives

T (t )  X0 eT (t ) + T1



⎞

t 0

(s )ds + T2



t 0

(s )dBH (s) (s ), (5.4)

eT

2 k M,

where denotes a unit row vector of size and the components of the vectors 1 and 2 are computed using Lemma 4.1 as follows

E E⎟

⎟ ⎟ ⎟ ⎠ E

E⎟,

(4.12)

L

⎞

⎞

⎜ M ⎟ ⎜  z − η  ⎟ ⎜ ⎟dz. ⎜ ηi − η ⎟ −1 ⎝ ⎠ =1



t

X (t )  T (t ),

namely the OM of integration of the

 ηj ⎜

z − η ⎟ M ⎜  ⎟ 1  ⎜ ⎟dz, li j = k+1 ⎜ ηi − η ⎟ 2 −1 ⎝ ⎠ =1  = i



! = [ξ11 ξ12 . . . ξ1M |ξ21 ξ22 . . . ξ2M | . . . |ξ k ξ k . . . ξ k ]T ,  2 1 2 2 2 M

(4.11)

in which the M × M matrices L and E are evaluated by

⎛

0

f (s, X (s ) )ds +

! is a 2k M known column vecwhere  (t) is defined in Eq. (4.6),  tor and  is a column vector of 2k M unknowns as follows

i = j.

(s )ds  P(1) (ξ ), 2k M × 2k M

t

! T (t ), t

Theorem 4.3 [46]. The integration of the vector  (ξ ) introduced in Eq. (4.6) can be represented as

0



where the second integral in the above relation is the stochastic integral including the V-OfBm. Now, we expand the identity function and unknown solution X(t) in terms of the CCWs as follows

 = i, j

 ξ

Now we go into detail about the method for solving SDE defined in Eq. (1.1). To do this, we first convert the problem into the following integral equation

X (t ) = X0 +

in which the entries of the M × M matrix Q are obtained by

⎧ M  2k+1 ⎪ ⎪ , ⎪ ⎪ ηi − ηr ⎪ ⎪ ⎪ r = 1 ⎪ ⎨r = i

5. Description of the method

(1 )i = f (ξi , xi ),

(2 )i = g(ξi , xi ),

i = 1, 2, . . . , 2k M, (5.5)

where i = (n − 1 )2k + m for n = 1, 2, . . . , 2k and m = 1, 2, . . . , M. By applying the integration by parts for the second integral in Eq. (5.4), we have



t 0

(s )dBH (s) (s ) = (t )BH (t ) (t ) − D(1)



0

t

(s )BH (s) (s )ds, (5.6)

where Theorem 4.2 is used. Also, the cardinal property of the Chebyshev cardinal wavelets gives

⎛

B H ( ξ1 ) ( ξ 1 ) 0 ⎜ ⎜ .. (t )BH (t ) (t )  ⎜ . ⎜ ⎝ 0

0

B H ( ξ2 ) ( ξ 2 ) .. . 0

(t )  (t ).

... ... .. . 0

0 0 .. .

⎞

⎟ ⎟ ⎟ ⎟ ⎠ H ( ξ2 k M ) B ( ξ2 k M ) (5.7)

110

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

From Eqs. (5.4)–(5.7) and employing Theorem 4.3, we achieve a residual function for the problem as follows





R(t ) = T − X0 eT − T1 P(1) − T2  − D(1)  P(1)

 (t ). (5.8)

Moreover, collocating the points ξ i for i = 1, 2, . . . , 2k M into Eq. (5.8) leads to the system of nonlinear algebraic equations as follows

  T − T1 P(1) − T2  − D(1)  P(1) = X0 eT .

(5.9)

Finally, we solve the obtained system and determine the elements of the vector  to insert in Eq. (5.2). In this case, the “fsolve” command of Maple 17 can be applied for solving the achieved system in Eq. (5.9). The method is described in the step-by-step Algorithm 3. Algorithm 3 The process of the proposed method. BH (t ) (t );

Z+ ,

Input: The numbers k ∈ M ∈ N; the V-OfBm X0 and the functions f, g : [0, 1] × R −→ R. T Output: The numerical solution: X (t ) 

 (t ). Step

1:

ηm = − cos

Generate

(2m−1 )π 2M

the value

% %2 |u|Hwm¯ ;N (a,b) = %u(m¯ ) %L2 (a,b) = |u|Hwm¯ (a,b) .

and

ξnm =

1, 2, . . . , 2 k .

n = 1, 2, . . . , 2k using Eq. (4.1). Step 3: Construct the vector (t ) as in Eq. (4.6). Step 4: Compute D(1 ) and P(1 ) using Eqs. (4.10) and (4.12), respectively. Step 5: Define the vector  with variable components. ! and e in Eqs. (5.3) and (5.4) like Eq. Step 6: Assign the vectors  (4.7). Step 7: Assign the vectors 1 and 2 using Eq. (5.5). Step 8: Assign the matrix  like  in Eq. (5.7).  Step 9: Put T − T1 P(1 ) − T2  − D(1 )  P(1 ) = X0 eT . Step 10: Solve the nonlinear system in Step 9 and compute the unknown vector .

(6.5)

w

2k M ¯ (0, 1 ) and I m Theorem 6.2 [46]. Let u ∈ Hw k,M u = n=1 m=1 n unm ψnm (t ), where unm = u(ξnm ). Then, the truncation error u − Ik,M u satisfies

u − Ik,M uL2wn (0,1) ≤ C¯m¯ (M − 1 )−m¯ $ k

1 2 j % %2 ¯ 2 m   %u( j ) % 2 n=1 j=min(m ¯ ,M )

where Ikn =

 n−1 2k

,

&1 / 2 ,

Lwn (Ikn )

2k+1

(6.6)



n 2k

. Moreover, in the maximum norm, we have

u − Ik,M uL∞ (0,1) ≤ Cˆm¯ (M − 1 )1/2−m¯ 2(k+1)/2 $ &1 / 2

1 2 j % %2 ¯ m  ( j ) %u % 2 max . k+1 L (I ) k n=1,2,...,2

(ηm + 2n − 1 ) for m = 1, 2, . . . , M and n = 2k+1 Step 2: Define the functions ψnm (t ) for m = 1, 2, . . . , M and 1

¯ ≤ N + 1, Remark 3. Note that |u|H m¯ ;N (a,b) ≤ uH m¯ (a,b) and for m w w we get

¯ ,M ) j=min(m

2

wn

(6.7)

kn

¯ (0, 1 ) be the analytical solution of Eq. m Theorem 6.3. Let X (t ) ∈ Hw n k (5.1) and Xnˆ (t ) where nˆ = 2 M is the CCWs solution of this equation. Furthermore, suppose that

(H1)

| f (t2 , X2 ) − f (t1 , X1 )| ≤ L1 |t2 − t1 |

+ L2 |X2 − X1 |,

(Lipschitz continuous),

(H2) |g(t2 , X2 ) − g(t1 , X1 )| ≤ L3 |t2 − t1 | + L4 |X2 − X1 |,

(Lipschitz continuous),

(6.8)

where t1 , t2 ∈ [0, 1], X1 , X2 ∈ R and Li (i = 1, 2, 3, 4 ) are positive constants. Then, we get

 % % X (t ) − Xnˆ (t )L∞ (0,1) ≤ L2 + L4 %BH (t ) (t )%L∞ (0,1) ρ (k, M, m¯ ),

(6.9) where

ρ (k, M, m¯ ) = Cˆm¯ (M − 1 )1/2−m¯ 2(k+1)/2 $

1 2 j % ¯ m %  %X ( j ) (t )%22 max

6. Convergence analysis To verify the validity of the CCWs method, we investigate the convergence analysis. To this end, we give some perquisite definitions and mathematical preliminaries. Definition 6.1 [58]. Assume (a, b) be a real bounded interval, w ¯ ≥ 0 be an integer. The Sobolev space be a weight function and m ¯ (a, b) is defined by m Hw

"

#

¯ , u( j ) (x ) ∈ L2w (a, b) . Hwm¯ (a, b) = u ∈ L2w (a, b) : for 0 ≤ j ≤ m (6.1) With the weighted inner product as

< u, v >m¯ ,w =

¯  m  j=0

b

u a

( j)

( j)

(x )v (x )w(x )dx,

uHwm¯ (a,b) =

¯ % m %  %u( j ) %22

(6.2)

&1 / 2

.

Lw (a,b)

j=0

(6.3)

|u|Hwm¯ ;N (a,b) =

¯ m  ¯ ,N+1 ) j=min(m

Lwn (Ikn )

2k+1

.

(6.10)

Proof. The approximate solution, Xnˆ (t ) satisfies

Xnˆ (t ) = X0 +



t 0

f (snˆ , Xnˆ (s ) )ds +



0

t

g(snˆ , Xnˆ (s ) )dBH (s ) (s ),

(6.11)

where snˆ is the wavelets expansion of the identity function. Subtracting Eq. (6.11) from Eq. (5.1), yields

X (t ) − Xnˆ (t ) =



0

t

( f (s, X (s ) ) − f (snˆ , Xnˆ (s ) ))ds

 0

t

(g(s, X (s ) ) − g(snˆ , Xnˆ (s ) ))dBH (s) (s ).

(6.12)

From Eqs. (6.8) and (6.12), and letting Y (t ) = t, we have

X (t ) − Xnˆ (t )L∞ (0,1) ≤ L1 Y (t ) − Ynˆ (t )L∞ (0,1) + L2 X (t ) − Xnˆ (t )L∞ (0,1)  + L3 Y (t ) − Ynˆ (t )L∞ (0,1) % % + L4 X (t ) − Xnˆ (t )L∞ (0,1) %BH (t ) (t )% ∞

L ( 0,1 )

Remark 2. To bound the approximation error, it is more appropriate to implement the following semi-norm

$

¯ ,M ) j=min(m

+

¯ (a, b) is a Hilbert space with the norm as m the Sobolev space Hw

$

n=1,2,...,2k

&1 / 2

% ( j ) %2 %u % 2

Lw (a,b)

&1 / 2

.

(6.4)

,

(6.13) where Ynˆ (t ) is the CCWs expansion of Y(t). Beside, choosing M > 1, from Eq. (6.7), we achieve

Y (t ) − Ynˆ (t )L∞ (0,1) = 0.

(6.14)

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Hence, Eqs. (6.7), (6.13) and (6.14) yield

CCWs method with the described method in [46] (where the VOfBm generated in the current paper is used instead of the classical fBm therein).

 % % X (t ) − Xnˆ (t )L∞ (0,1) ≤ L2 + L4 %BH (t )%L∞ (0,1) ρ (k, M, m¯ ), (6.15) 

and then the proof is completed. Corollary 6.4. Theorem 6.3 also yields

X (t ) − Xnˆ (t )L∞ (0,1) −→ 0,

as nˆ −→ ∞,

111

Example 1 [46]. Consider the SDE rely on the V-OfBm as

dX (t ) = σ¯ 2 cos(X (t )) sin (X (t )) dt − σ¯ sin (X (t )) dBH (t ) (t ), 3

X (0 ) = X0 , with the exact solution

(6.16)

which confirms the validity of the CCWs approximation method. 7. Numerical experiments To investigate the accuracy and applicability of the proposed method, some numerical examples are considered. Throughout the ˆ = 400 for constructing the numerical experiments, we assign m V-OfBm. Also, we compare the absolute error (AE) values of the

2





X (t ) = arccot σ¯ BH (t ) (t ) + cot (X0 ) . The method of this paper with M = 3 and three values k = 3, 4, 5 is used to solve this problem where σ¯ = X0 = 1/20. The AE values of the CCWs method and the one described in [46], for two selected H (t ) = 0.5 + 0.3 sin(π t ) and H (t ) = 0.6 − 0.2 exp(−2t ) are listed in Tables 1 and 2, respectively. In this case, for constructing the VOfBm, the assumed N is 120. The exact and approximate solutions are illustrated in Fig. 2, and also the log-plots of the AE functions for different k and H(t) are drawn in Fig. 3.

Fig. 2. The behavior of the results obtained by the CCWs method for two selected H(t) in Example 1.

Fig. 3. Logarithmic scale of the AE functions of the CCWs method for two selected H(t) in Example 1.

112

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 4. The behavior of the results obtained by the CCWs method for two selected H(t) in Example 2.

Fig. 5. Logarithmic scale of the AE functions of the CCWs method for two selected H(t) in Example 2. Table 1 The AE values of the CCWs method and the one described in [46] with M = 3 and H (t ) = 0.5 + 0.3 sin(π t ), where N = 120 at some points t ∈ [0, 1] in Example 1. The method of [46]

The proposed method

t

k=3

k=4

k=5

k=3

k=4

k=5

0.2 0.4 0.6 0.8 1.0

1.1093 × 10−5 1.6179 × 10−5 1.2966 × 10−5 1.9769 × 10−5 1.6575 × 10−5

2.6152 × 10−5 2.4825 × 10−5 2.8235 × 10−5 2.5299 × 10−5 2.3425 × 10−5

1.5015 × 10−5 1.9955 × 10−5 2.2597 × 10−5 2.6019 × 10−5 2.7041 × 10−5

2.7463 × 10−6 6.2365 × 10−6 1.3982 × 10−6 3.2472 × 10−6 2.8907 × 10−6

3.8793 × 10−6 2.4425 × 10−7 4.9411 × 10−7 3.1180 × 10−6 2.5445 × 10−6

1.4045 × 10−6 3.9707 × 10−7 3.4729 × 10−7 9.6406 × 10−7 3.9374 × 10−7

Regarding to Tables 1 and 2, and Figs. 2 and 3, it can be concluded that by increasing k the accuracy is improving. Precisely, the approximate solutions converge to the analytical solutions by increasing the number of the CCWs as it is already proven in Section 6. Moreover, the accuracy of results obtained by using the CCW method is higher than the accuracy of the results achieved by the method in [46].

Example 2 [46]. Consider the SDE rely on the V-OfBm as









dX (t ) = σ¯ 2 X (t ) X 2 (t ) − 1 dt + σ¯ 1 − X 2 (t ) dBH (t ) (t ), X (0 ) = X0 , with the exact solution

  (1 + X0 ) exp 2σ¯ BH (t ) (t ) + X0 − 1   X (t ) = . (1 + X0 ) exp σ¯ BH (t ) (t ) − X0 + 1

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

113

Table 2 The AE values of the CCWs method and the one described in [46] with M = 3 and H (t ) = 0.6 − 0.2 exp(−2t ), where N = 120 at some points t ∈ [0, 1] in Example 1. The method of [46] t 0.2 0.4 0.6 0.8 1.0

The proposed method

k=3

k=4 −6

4.2790 × 10 1.2906 × 10−5 1.7811 × 10−5 3.8843 × 10−5 2.8988 × 10−5

k=5 −5

4.7035 × 10 4.4857 × 10−5 5.9377 × 10−5 6.2745 × 10−5 6.2068 × 10−5

k=3 −5

3.2261 × 10 4.6245 × 10−5 6.1916 × 10−5 7.2320 × 10−5 7.0370 × 10−5

k=4 −6

5.0 0 08 × 10 2.8066 × 10−5 5.1922 × 10−6 5.0562 × 10−6 6.1019 × 10−6

k=5 −6

2.9191 × 10 2.6119 × 10−6 2.1989 × 10−6 4.9267 × 10−6 2.6561 × 10−6

1.8266 × 10−6 1.0112 × 10−6 1.5571 × 10−6 1.7644 × 10−7 1.7949 × 10−6

Table 3 The AE values of the CCWs method and the one described in [46] with M = 2 and H (t ) = 0.7 + 0.2 sin(π t ), where N = 90 at some points t ∈ [0, 1] in Example 2. The method of [46]

The proposed method

t

k=4

k=5

k=6

k=4

k=5

k=6

0.2 0.4 0.6 0.8 1.0

1.3372 × 10−3 2.1738 × 10−3 2.7887 × 10−3 3.4938 × 10−3 6.4426 × 10−3

1.6538 × 10−4 6.0265 × 10−4 3.3729 × 10−4 5.8469 × 10−4 1.6964 × 10−3

8.9095 × 10−4 9.8840 × 10−5 5.7934 × 10−4 7.9824 × 10−4 3.1934 × 10−4

3.2850 × 10−5 7.7123 × 10−5 1.1454 × 10−4 3.7300 × 10−6 1.0563 × 10−3

2.4033 × 10−4 1.2578 × 10−5 5.1047 × 10−5 3.6659 × 10−6 5.0603 × 10−4

1.4634 × 10−5 8.2132 × 10−6 4.6802 × 10−5 2.9566 × 10−6 3.0478 × 10−4

Table 4 The AE values of the CCWs method and the one described in [46] with M = 2 and H (t ) = 0.7 − 0.25 exp(−t ), where N = 90 at some points t ∈ [0, 1] in Example 2. The method of [46] t 0.2 0.4 0.6 0.8 1.0

The proposed method

k=4

k=5 −3

9.2683 × 10 9.4482 × 10−3 1.0324 × 10−2 1.2324 × 10−2 1.6127 × 10−2

k=6 −3

4.1978 × 10 3.5887 × 10−3 3.2583 × 10−3 2.0882 × 10−3 3.1585 × 10−4

k=4 −3

4.6954 × 10 2.2058 × 10−3 1.4168 × 10−3 5.9840 × 10−4 2.8317 × 10−4

k=5 −4

1.0811 × 10 2.4154 × 10−4 4.9532 × 10−4 1.9303 × 10−4 1.9402 × 10−3

k=6 −4

6.9681 × 10 1.4780 × 10−4 8.0529 × 10−5 1.6314 × 10−4 8.4869 × 10−4

1.2367 × 10−4 3.6764 × 10−4 5.7225 × 10−4 6.4132 × 10−5 7.1071 × 10−4

Table 5 The AE values of the CCWs method and the one described in [46] with M = 2 and H (t ) = 0.5 + 0.3 sin(2t ), where N = 125 at some points t ∈ [0, 1] in Example 3. The method of [46]

The proposed method

t

k=4

k=5

k=6

k=4

k=5

k=6

0.2 0.4 0.6 0.8 1.0

5.9187 × 10−3 8.2297 × 10−3 1.0692 × 10−2 1.0921 × 10−2 1.1450 × 10−2

4.5700 × 10−3 1.3903 × 10−3 1.4159 × 10−3 3.2230 × 10−3 1.2380 × 10−3

3.5653 × 10−3 2.7855 × 10−3 4.3962 × 10−3 7.2613 × 10−3 6.7031 × 10−3

1.0894 × 10−3 3.1306 × 10−4 2.0265 × 10−4 1.4900 × 10−5 1.0142 × 10−3

1.0753 × 10−3 1.7989 × 10−4 1.5228 × 10−4 2.9180 × 10−4 4.2322 × 10−4

3.8017 × 10−5 1.4971 × 10−4 1.8112 × 10−4 8.3251 × 10−5 1.0831 × 10−4

Table 6 The AE values of the CCWs method and the one described in [46] with M = 2 and H (t ) = 0.4 + 0.3t 2 , where N = 125 at some points t ∈ [0, 1] in Example 3. The method of [46] t 0.2 0.4 0.6 0.8 1.0

The proposed method

k=4

k=5 −3

7.0546 × 10 2.5215 × 10−4 9.8373 × 10−3 2.0408 × 10−2 8.5276 × 10−3

k=6 −3

4.3165 × 10 9.5217 × 10−3 1.1411 × 10−2 6.7045 × 10−3 1.0984 × 10−2

k=4 −3

2.4885 × 10 1.6036 × 10−2 1.0883 × 10−2 5.3272 × 10−3 1.0050 × 10−2

The proposed method is applied for solving this problem with M = 2 and different values of k for chosen σ¯ = 1/30 and X0 = 0.01. The AE values of the CCW method and the one described in [46] for two selected H (t ) = 0.7 + 0.2 sin(π t ) and H (t ) = 0.7 − 0.25 exp(−t ), are respectively shown in Tables 3 and 4. In this case, for constructing the V-OfBm, the assumed N is 90. The exact and approximate solutions are illustrated in Fig. 4, and also the log-plots of the AE functions for different k and H(t) are drawn in Fig. 5.

k=5 −3

3.3143 × 10 7.3848 × 10−3 1.1772 × 10−3 1.6250 × 10−4 1.6101 × 10−3

k=6 −3

4.2019 × 10 7.3124 × 10−4 3.6961 × 10−4 1.1570 × 10−3 5.5401 × 10−4

9.9881 × 10−5 5.3348 × 10−4 1.3488 × 10−4 1.7675 × 10−4 1.2658 × 10−4

Figs. 4 and 5, and Tables 3 and 4 conclude that by increasing the number of CCWs, the accuracy is exponentially improving. Moreover, Tables 3 and 4 show that this approach is more accurate than the one described in [46]. Example 3 [34]. Consider the SDE rely on the V-OfBm

dX (t ) =

' 1 2 σ¯ X (t ) dt + σ¯ 1 + X 2 (t ) dBH (t ) (t ), 2

X (0 ) = X0 ,

114

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 6. The behavior of the results obtained by the CCWs method for two selected H(t) in Example 3.

Fig. 7. Logarithmic scale of the AE functions using the CCWs method for two selected H(t) in Example 3.

with the exact solution





X (t ) = sinh σ¯ BH (t ) (t ) + arcsinh(X 0 ) . We solve the system using the CCWs method with M = 2 and different values of k for chosen σ¯ = 1/20 and X0 = 0. The AE values of the CCW method and the one described in [46] for two choices of H (t ) = 0.5 + 0.3 sin(2t ) and H (t ) = 0.4 + 0.3t 2 are shown in Tables 3. In this case, for constructing the V-OfBm, the assumed N is 125. The exact and approximate solutions are illustrated in Fig. 6, and also the log-plots of the AE functions for different k and H(t) are drawn in Fig. 7. The obtained results confirm the appropriateness of CCW method for solving this kind of problems. 8. Applications in applied problems To investigate the practical applications of the CCWs method, we discuss some famous stochastic models occurred in the real world.

8.1. Mathematical ecology Stochastic logistic equation Logistic equation is a famous model in the ecology science. The classic type of it is represented in [59] as follows





X (t ) dX (t ) = rX (t ) 1 − , dt κ

X (0 ) = X0 ,

where X(t) denotes the size of population, κ > 0 express the carrying capacity of environment and r is the growth rate of population. It is investigated in [59] that lim X (t ) = 0 for r < 0 and lim X (t ) = t→∞

t→∞

κ for r > 0. Let environmental noise affects on the growth rate r as

r → r + σ¯ ξ¯ (t ), where ξ¯ (t ) = dBdt(t ) and B(t) is a standard Bm, and σ¯ is a given constant. It enhances the classic logistic equation to the following

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

115

Fig. 8. The behavior of the results achieved by the CCWs method for two selected H(t) in the logistic equation.

Fig. 9. The behavior of the results achieved by the CCWs method for two selected H(t) in the population growth model.

stochastic model



dX (t ) = rX (t ) 1 −

X (t )





κ

dt + σ¯ X (t ) 1 −

X (t )

The proposed method is applied with (k = 5, M = 5 ) for solving Eq. (8.1) for two choices of H (t ) = 0.5 + 0.3 cos(10 0 0t ) and H (t ) = 0.3 + 0.3 exp(−t ) with N = 350. The behavior of the numerical results for X0 = 0.3, r (t ) = 0.2, κ = 1 and two choices of σ¯ (t ) are illustrated in Fig. 8.



κ

dB(t ),

X (0 ) = X0 . Since the natural growth of populations commonly vary with respect to t, studying the following non-autonomous form of stochastic logistic equations is preferred [59]



dX (t ) = r (t )X (t ) 1 −





X (t )

dt + σ¯ (t )X (t ) 1 −

κ

X (t )

κ



dB(t ),

X (0 ) = X0 . Furthermore, we can introduce the stochastic model of the logistic equation with the V-OfBm by upgrading B(t) to BH(t) (t) as follows



dX (t ) = r (t )X (t ) 1 − X (0 ) = X0 .

X (t )

κ





dt + σ¯ (t )X (t ) 1 −

X (t )

κ



dBH (t ) (t ),

(8.1)

Stochastic population growth model in a closed system One of the studied model for balancing the population size is the standard logistic growth affected by the build up of toxins. This deleterious effect on the size of population mathematically can be described as follows [32]

du = γ u − ζ u − u dt¯



t¯ 0

u(s )ds,

u ( 0 ) = u0 ,

(8.2)

where the parameter γ denotes the coefficient of the birth rate, ζ > 0 indicates the coefficient of the crowding and the term

116

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 10. The behavior of the results achieved by the CCWs method for three selected H(t) in the Lotka–Volterra system.

Fig. 11. The behavior of the results achieved by the CCWs method for some values of σ¯ in the Brusselator model.

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 12. The behavior of the results achieved by the CCWs method for some values of σ¯ in the Brusselator model.

Fig. 13. The behavior of the results obtained by the CCWs method for κ = 11 in the Duffing-Van der Pol oscillator problem.

117

118

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 14. The behavior of the results achieved by the CCWs method for κ = 12 in the Duffing-Van der Pol oscillator model.

Fig. 15. The behavior of the results obtained by the CCWs method for κ = 13 in the Duffing-Van der Pol oscillator problem.

including the integral denotes the “total amount” or a “total metabolism” of toxins produced since the initial time. By change  γζ of variable t = ζ t¯ and defining X = γζ u, and by putting κ =  , we get



dX (t ) = κ X (t ) 1 − X (t ) − dt



t 0



X (s )ds ,

dX (t ) = rX (t ) 1 − X (t ) −

 × 1 − X (t ) −



t 0



t 0

X (0 ) = X0 .



 × 1 − X (t ) −



t 0





t 0

X (0 ) = X0 . (8.3)



1i = f (ξi , xi , x¯i ),

2i = g(ξi , xi , x¯i ),

i = 1, 2, . . . , 2k M,

¯ and where x¯i is the ith component of the vector 













f t , X (t ), X¯ (t ) = rX (t ) 1 − X (t ) − X¯ (t ) ,





As a numerical example, we have solved Eq. (8.4) where r = 1 and X0 = 0.1 for two values of σ¯ by applying the CCWs method with (k = 6, M = 2 ). The numerical solutions for two selected H (t ) = 0.5 + 0.4 sin(20 0 0t ) and H (t ) = 0.6 − 0.5 exp(−t ) (where N = 350 is considered for constructing the V-OfBm) are drawn in Fig. 9. 8.2. Biological systems

X (s )ds dt + σ¯ X (t )



X (s )ds dBH (t ) (t ),

¯ T (t ), X (s )ds  T P (1) (t )  

g t , X (t ), X¯ (t ) = σ¯ X (t ) 1 − X (t ) − X¯ (t ) .

X (s )ds dt + σ¯ X (t )

The fractional stochastic model driven with the V-OfBm for Eq. (8.3) can be written as

dX (t ) = rX (t ) 1 − X (t ) −

t

and



X (s )ds dB(t ),



0

The stochastic model of the above equation can be obtained by putting κ = r + σ¯ ξ¯ (t ) in which ξ¯ (t ) = dBdt(t ) and σ¯ is a nonstochastic coefficient which indicates the intensity of noise. Hence, we have



The CCWs method can simply be applied for solving Eq. (8.4). Note that by expanding the solution in terms of the CCWs as X(t) ࣃ T  (t), we have

X (0 ) = X0 . (8.4)

Lotka–Volterra is another well-known nonlinear system appeared in biology [60]. This system proposed for the observed periodic variations in a predator-prey system. One of the simplest stochastic model of the Lotka–Volterra system is described

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 16. The behavior of the results achieved by the CCWs method for κ = 14 in the Duffing-Van der Pol oscillator model.

Fig. 17. The behavior of the results achieved by the CCWs method for κ = 11, 12, 13, 14 in the Duffing-Van der Pol oscillator model.

119

120

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 18. The behavior of the results achieved by the CCWs method for some values of σ¯ in the pendulum model.

in [34,60] as

We first transform Eq. (8.6) into the following equivalent stochastic system



dX (t ) = (b1 − a1Y (t ) )X (t )dt + σ¯ 1 X (t )dB1 (t ),

X (0 ) = X0 ,

dY (t ) = (b2 − a2 X (t ) )Y (t )dt + σ¯ 2Y (t )dB2 (t ),

Y (0 ) = Y0 ,

where X(t) indicates the preys number and Y(t) express the predators number, and B1 (t) and B2 (t) are two independent Brownian motions. Here, we study a fractional stochastic model driven with the V-OfBm of Lotka–Volterra system as



(t ) dX (t ) = (b1 − a1Y (t ) )X (t )dt + σ¯ 1 X (t )dBH (t ), 1

dY (t ) = (b2 − a2 X (t ) )Y (t )dt + σ¯ 2Y (t )

(t ) dBH 2

(t ),

X (0 ) = X0 , Y (0 ) = Y0 .

(8.5)

The CCWs method described in Section 5 can simply be expanded for solving Eq. (8.5). To this end and the next uses, we implement the CCWs method for the following general form



⎧  t  t (s ) ⎪ ⎨X (t ) = X0 + f1 (s, X (s ), Y (s ) )ds + g1 (s, X (s ), Y (s ) )dBH ( s ), 1 ⎪ ⎩Y (t ) = Y0 +

 0t 0

f2 (s, X (s ), Y (s ) )ds +

 0t 0

(s ) g2 (s, X (s ), Y (s ) )dBH ( s ). 2

(8.7) Now, we expand the identity function and X(t) and Y(t) in terms of the CCWs as follows

! T (t ), t

X (t )  T1 (t ),

Y (t )  T2 (t ),

(8.8)

! is a 2k M known column vecwhere  (t) is defined in Eq. (4.6),  tor (the CCWs coefficients vector of the unit function), and 1 and 2 are column vectors of 2k M unknowns as follows

! = [ξ11 ξ12 . . . ξ1M |ξ21 ξ22 . . . ξ2M | . . . |ξ k ξ k . . . ξ k ]T ,  2 1 2 2 2 M (t ) dBH 1

(t ),

X (0 ) = X0 ,

1 = [x11 x12 . . . x1M |x21 x22 . . . x2M | . . . |x2k 1 x2k 2 . . . x2k M ]T ,

(t ) dY (t ) = f2 (t, X (t ), Y (t ) )dt + g2 (t, X (t ), Y (t ) )dBH (t ), 2

Y (0 ) = Y0 .

2 = [y11 y12 . . . y1M |y21 y22 . . . y2M | . . . |y2k 1 y2k 2 . . . y2k M ]T .

dX (t ) = f1 (t, X (t ), Y (t ) )dt + g1 (t, X (t ), Y (t ) )

(8.6)

(8.9)

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Substituting Eq. (8.8) into Eq. (8.7) leads to

⎧  t  t ⎪ T T T T ⎪  ( t )  X e ( t ) +  ( s ) ds +  (s )dBH1 (s) (s ), 0 ⎨ 1 11 12 0

0

0

0

 t  t ⎪ ⎪ ⎩T2 (t )  Y0 eT (t ) + T21 (s )ds + T22 (s )dBH2 (s) (s ), (8.10) where the components of the vectors j1 and j2 for j = 1, 2 are computed as follows

( j1 )i = f j (ξi , xi , yi ),

( j2 )i = g j (ξi , xi , yi ),

j = 1, 2, i = 1, 2, . . . , 2k M,

(8.11)

in which i = (n − 1 )2k + m, n = 1, 2, . . . , 2k and m = 1, 2, . . . , M. Applying the integration by parts for the second integrals in Eq. (8.10) yields

⎧ t  t ⎪ ⎨ (s )dBH1 (s) (s ) = (t )BH1 (t ) (t ) − D(1) (s )BH1 (s) (s )ds, 0 t

⎪ ⎩

0

(s )

(s ) dBH 2

(s ) = (t )

(t ) BH 2

(t ) − D

(1 )

0 t 0

(s )BH2 (s) (s )ds, (8.12)

where Theorem 4.2 is used. Also, the cardinal property of the Chebyshev cardinal wavelets gives

⎛

H ξ B1 ( 1 ) (ξ1 ) ⎜ 0 ⎜ ⎜ . H (t ) . (t )B1 (t )  ⎜ . ⎜

⎝

0

. . .

... ... .. .

0

0

0

H ξ B1 ( 2 )

(ξ2 )

⎟ ⎟ ⎟ ⎟ ⎟ H ( ξ2 k M ) ⎠ B1 (ξ2k M )

(t )  1 (t ), and

(8.13)

⎛

H ξ B2 ( 1 ) (ξ1 ) ⎜ 0 ⎜ ⎜ . H (t ) . (t )B2 (t )  ⎜ . ⎜

⎝

⎞

0 0 . . .

0

. . .

... ... .. .

0

0

0

H ξ B2 ( 2 )

(ξ2 )

⎞

0 0 . . .

⎟ ⎟ ⎟ ⎟ ⎟ H ( ξ2 k M ) ⎠ B2 (ξ2k M )

(t )  2 (t ).

(8.14)

From Eqs. (8.10)–(8.14) and employing Theorem 4.3, we achieve the residual functions for the problem as follows







 (t ),  T   T T ( 1) T ( 1) ( 1) R2 (t ) = 2 − Y0 e − 21 P − 22 2 − D 2 P (t ). R1 (t ) = T1 − X0 eT − T11 P(1) − T12 1 − D(1) 1 P(1)

(8.15) Moreover, by taking the collocation points ξ i for i = 1, 2, . . . , 2k M into Eq. (8.15), we achieve the following system of algebraic equations



  T1 − T11 P(1) − T12 1 − D(1) 1 P(1) = X0 eT ,   T2 − T21 P(1) − T22 2 − D(1) 2 P(1) = Y0 eT .

(8.16)

Finally, by solving the above system and determining the elements of the vectors 1 and 2 , and inserting in Eq. (8.8), we get an approximate solution. The CCWs method is described concisely in Algorithm 4. As a numerical example, we have implemented the method explained in this section with (k = 6, M = 2 ) for solving Eq. (8.5), where a1 = 0.3, a2 = 0.1, b1 = 2.0, b2 = 1.5, σ¯ 1 = 0.2, σ¯ 2 = 0.4, X0 = 1 and Y0 = 1. The behaviors of the numerical solutions for three choices H (t ) = 0.5, H (t ) = 0.5 + 0.3 sin(50t ) and H (t ) = 0.5 + 0.3 sin(150t ) (where N = 150 is considered in constructing the V-OfBm) are plotted in Fig. 10.

121

Algorithm 4 The process of the proposed method for the stochastic model defined in Eq. (8.6). (t ) Input: The numbers k ∈ Z+ , M ∈ N; the V-OfBms BH (t ) and 1

(t ) BH (t ); the values X0 and Y0 , and the functions 2 f1 , f2 , g1 , g2 : [0, 1] × R × R −→ R. Output: The numerical solutions: X (t )  T1 (t ) Y (t )  T2 (t ).

Step

1:

Generate

ηm = − cos

(2m−1 )π 2M



and

and

ξnm =

(ηm + 2n − 1 ) for m = 1, 2, . . . , M and n = 2k+1 Step 2: Define the functions ψnm (t ) for m = 1, 2, . . . , M and n = 1, 2, . . . , 2k using Eq. (4.1). Step 3: Construct the vector (t ) as in Eq. (4.6). Step 4: Compute the matrices D(1 ) and P(1 ) using Eqs. (4.10) and (4.12), respectively. Step 5: Define the vectors 1 and 1 with variable components. ! and e in Eqs. (8.8) and (8.8) like Eq. (4.7). Step 6: Assign  Step 7: Assign the vectors  j1 and  j2 for j = 1, 2 using Eq. (8.11). Step 8: Assign the matrices 1 and 2 like in Eqs. (8.13) and (8.14).  T   1 − T11 P(1) − T12 1 − D(1) 1 P(1)  = X0 eT , Step 9: Put . T2 − T21 P(1) − T22 2 − D(1) 2 P(1) = Y0 eT . Step 10: Solve the algebraic system in the previous step and determine 1 and 2 . 1

1, 2, . . . , 2 k .

8.3. Brusselator problem The fractional stochastic Brusselator equation driven with the V-OfBm can be defined as follows

  ⎧ 2 dX (t ) = ( − 1 )X (t ) + (1 + X (t ) ) Y (t ) dt ⎪ ⎪ ⎪ ⎨ +σ¯ X (t )(1 + X (t ) )dBH (t ) (t ),   2 ⎪ dY (t ) = −  X (t ) + (1 + X (t ) ) Y (t ) dt ⎪ ⎪ ⎩ −σ¯ X (t )(1 + X (t ) )dBH (t ) (t ),

X (0 ) = X0 , Y (0 ) = Y0 , (8.17)

where σ¯ and ϖ are real constants. It is worth to mention that the deterministic Brusselator model (σ¯ = 0) was developed for a particular model for bifurcations in chemical reactions [61]. We remind that the standard stochastic Brusselator model is investigated in [34,61]. The method proposed in Section 8.2 can simply be used for the model expressed in Eq. (8.17). As a particular case, we examine Eq. (8.17) with  = 2 and some values of σ¯ , starting at (X0 , Y0 ) = (−0.1, 0.0 ). As a numerical simulation, we used the proposed method with (k = 5, M = 3 ) for this problem where H (t ) = 0.5 + 0.3 sin(40 0 0t ) and N = 290. The behavior of the numerical results are plotted in Figs. 11 and 12. 8.4. Physics problems The method proposed in Section 8.2 is successfully employed for solving some well-known stochastic equations in the physics. Duffing-Van der Pol Oscillator A common type of the stochastic Duffing-Van der Pol oscillator driven with white noise, ξ¯ (t ) = dBdt(t ) is as follows [62]

x¨ + x˙ −



  − x2 x = σ¯ xξ¯ ,

(8.18)

where ϖ is a real constant and σ¯ controls the strength of the induced multiplicative noise. Note that the corresponding twodimensional form includes the components X and Y relying on the

122

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

Fig. 19. The behavior of the results achieved by the CCWs method for some values of σ¯ in the pendulum model.

displacement x and speed x˙ , respectively which can be represented as follows



dX (t ) = Y (t )dt ,

X (0 ) = X0 ,

dY (t ) =

Y (0 ) = Y0 .



   − X 2 (t ) X (t ) − Y (t ) dt + σ¯ X (t )dB(t ),

(8.19) In addition, we generalize the above system to a fractional stochastic one using the V-OfBm as follows



dX (t ) = Y (t )dt , dY (t ) =







 − X 2 (t ) X (t ) − Y (t ) dt + σ¯ X (t )dBH (t ) (t ),

X (0 ) = X0 , Y (0 ) = Y0 .

(8.20) The stochastic system generated in Eq. (8.20) can be solved using the CCWs method discussed in Section 8.2. For a numerical simulation, we investigate the stochastic problem (8.20) with  = 0 and various values σ¯ = 0.0 (deterministic problem) and σ¯ = 0.5, 1.0 (stochastic problems) on the interval [0, 8], starting at (X0 , Y0 ) = (−κε , 0 ) for κ = 11, 12, 13, 14 and ε = 0.2. The obtained results using the CCWs method with (k = 6, M = 3 ) where H (t ) = 0.5 + 0.4 sin(20 0 0t ) and N = 350 are shown in Figs. 13–17.

Pendulum problem The nonlinear stochastic pendulum equation is given in [61] by

x¨ + σ¯ x˙ ξ¯ + sin(x ) = 0,

x ( 0 ) = x0 ,

x˙ (0 ) = x˙ 0 ,

(8.21)

where σ¯ is a real constant and ξ¯ (t ) = dBdt(t ) . The above classical SDE can be rewritten in a canonical stochastic system as follows



dX (t ) = Y (t )dt , dY (t ) = − sin (X (t ) )dt − σ¯ Y (t )dB(t ),

X ( 0 ) = x0 , Y (0 ) = x˙ 0 .

(8.22)

The stochastic system defined in Eq. (8.22) can be developed to a fractional stochastic one as follows



dX (t ) = Y (t )dt , dY (t ) = − sin (X (t ) )dt − σ¯ Y (t )dBH (t ) (t ),

X ( 0 ) = x0 , Y (0 ) = x˙ 0 .

(8.23)

The method proposed in Section 8.2 with (k = 6, M = 3 ) and H (t ) = 0.5 + 0.4 sin (30 0 0t ) where N = 360 is applied for solving this problem. The behaviors of the oscillator with (x0 , x˙ 0 ) = (π /3, 1 ) for some various values of σ¯ on [0, 3π ] are plotted in Figs. 18 and 19.

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

9. Conclusion In this paper, a new technique introduced for generating variable-order fractional Brownian motion (V-OfBm). The blockpulse functions (BPFs) and the hat functions (HFs) simultaneously are employed in the process of constructing this stochastic process. An innovative computational method expressed for solving a common type of nonlinear stochastic differential equations (SDEs) involving the V-OfBm. The scheme is based on the Chebyshev cardinal wavelets (CCWs) and their operational matrices of differentiation and integration. The proposed algorithm transforms the intended problem into a system of algebraic equations. We discussed the convergence of the CCWs method in the Sobolev space. Some numerical examples investigated to examine the accuracy and applicability of the suggested method. The obtained results confirm that the new scheme is an applicable and useful tool for solving such nonlinear stochastic problems. Moreover, to show the significance of the CCWs method, we illustrated the numerical procedure for solving the stochastic type of some well-known mathematical models, including the logistic model, population growth model in a closed system, and Lotka–Volterra, Brusselator, Duffing-Van der Pol oscillator and pendulum problems. In addition to the mentioned privileges, the approach can be extended for solving more various types of stochastic problems. Declaration of competing interest In this paper, firstly, a new computational technique is established to construct the variable-order fractional Brownian motion (V-OfBm). The operational matrix (OM) of derivative of the BPFs and the OM of variable-order (V-O) fractional integration of the HFs are employed for constructing the V-OfBm. Then, a new computational method based on the Chebyshev cardinal wavelets is proposed to solve a common class of nonlinear stochastic differential equations (SDEs) driven by (V-OfBm). These basis functions possess many valuable properties, such as spectral accuracy, orthogonality and cardinality. Using the interpolation property of these wavelets, a new algorithm is presented for computing nonlinear terms in such problems. The main advantage of the proposed method is the reduction of the problem to a simpler one, which consists of solving a system of nonlinear algebraic equations. The convergence analysis of the expressed method is proved. Moeover, the useful applications of the proposed method are illustrated by solving some well-known stochastic models, i.e. stochastic logistic problem, stochastic population growth model in a closed system, stochastic Lotka–Volterra problem, stochastic Brusselator problem, stochastic Duffing-Van der Pol oscillator problem and stochastic pendulum problem. References [1] Zeng C, Yang Q, Chen YQ. Synthesis of multi fractional gaussian noises based on variable-order fractional operators. Signal Process 2011;91:1645–50. [2] Kulish VV, Lage JL. Application of fractional calculus to fluid mechanics. J Fluids Eng 2002;124(3):803–6. [3] Engheta N. On fractional calculus and fractional multipoles in electromagnetism. Antennas Propag 1996;44:554–66. [4] Lederman C, Roquejoffre JM, Wolanski N. Mathematical justification of a nonlinear integro-differential equation for the propagation of spherical flames. Ann di Matem 2004;183:173–239. [5] Bagley RL, Torvik PJ. Fractional calculus in the transient analysis of viscoelastically damped structures. AIAA J 1985;23:918–25. [6] Meral FC, Royston TJ, Magin R. Fractional calculus in viscoelasticity: an experimental study. Commun Nonlinear Sci Numer Simul 2010;15:939–45. [7] Heydari MH, Hooshmandasl MR, Ghaini FMM. An efficient computational method for solving fractional biharmonic equation. Comput Math Appl 2014;68(9):269–87. [8] Heydari MH, Hooshmandasl MR, Maalek Ghaini FM, Cattani C. Wavelets method for solving fractional optimal control problems. Appl Math Comput 2016;286:139–54.

123

[9] Heydari MH, Hooshmandasl MR, Maalek Ghaini FM, Cattani C. Wavelets method for the time fractional diffusion-wave equation. Phys Lett A 2015;379:71–6. [10] Heydari MH, Hooshmandasl MR, Mohammadi F. Legendre wavelets method for solving fractional partial differential equations with dirichlet boundary conditions. Appl Math Comput 2014;234:267–76. [11] Heydari MH, Hooshmandasl MR, Mohammadi F, Cattani C. Wavelets method for solving systems of nonlinear singular fractional volterra integro-differential equations. Commun Nonlinear Sci Numer Simul 2014;19:37–48. [12] Heydari MH, Hooshmandasl MR, Maalek Ghaini FM, Li M. Chebyshev wavelets method for solution of nonlinear fractional integrodifferential equations in a large interval. Adv Math Phys 2013;2013:12. doi:10.1155/2013/482083. (Article ID 482083) [13] He JH. Fractal calculus and its geometrical explanation. Results Phys 2018;10:272–6. [14] Li XX, Tian D, He CH, He JH. A fractal modification of the surface coverage model for an electrochemical arsenic sensor. Electroch Acta 2019;296:491–3. [15] Wang Q, Shi X, He JH, Li Z. Fractal calculus and its application to explanation of biomechanism of polar bear hairs. Fractals 2018;26(6). doi:10.1142/ S0218348X1850086X. [16] He JH. A tutorial review on fractal space-time and fractional calculus. Int J Theor Phys 2014;53(11). [17] He JH, Elagan SK, Li ZB. Geometrical explanation of the fractional complex transform and derivative chain rule for fractional calculus. Phys Lett A 2012:376. [18] Wang Y, Zhang YF, Liu ZJ. An explanation of local fractional variational iteration method and its application to local fractional modified KDV equation. Therm Sci 2018;22(1A). [19] Wang KL, Liu SY. He’S fractional derivative and its application for fractional Fornberg-Whitham equation. Therm Sci 2017;21(5). [20] Mandelbrot BB, Thornton T, Mallon MW. The variation of certain speculative prices. J Bus 1963;36:394–419. [21] Gilden DL, Thornton T, Mallon MW. 1/F noise in human cognition. Science 1995;267:1837–9. [22] Perez DG, Zunino L, Garavaglia M. Modeling turbulent wavefront phase as a fractional Brownian motion: a new approach. J Opt Soc Am 2004;21:1962–9. [23] Osorio I, Frei M. Hurst parameter estimation for epileptic seizure detection. Commun Inf Syst 2007;7:167–76. [24] Mishura YS. Stochastic calculus for fractional Brownian motion and related processes. Berlin: Springer; 2008. [25] Duncan TE, Hu Y, Pasik-Duncan B. Stochastic calculus for fractional Brownian motion. i. theory. SIAM J Control Optim 20 0 0;38:582–612. [26] Alòs E, Mazet O, Nualart D. Stochastic calculus with respect to gaussian processes. Ann Probab 20 0 0;29:766–801. [27] Elliott RC, Hoek JV. A general fractional white noise theory and applications to finance. Math Finance 2003;13:301–30. [28] Carmona P, Coutin L, Montseny G. Stochastic integration with respect to fractional Brownian motion. Ann I H Poincar Probab Stat 2003;31:27–68. [29] Biagini F, Øksendal B, Sulem A, Wallner N. An introduction to white noise theory and Malliavin calculus for fractional Brownian motion. Proc R Soc 2004;460:347–72. [30] Jolis M. On the wiener integral with respect to the fractional Brownian motion on an interval. J Math Anal Appl 2007;330:1115–27. [31] Heydari MH, Hooshmandasl MR, Maalek Ghaini FM, Cattani C. A computational method for solving stochastic itô-volterra integral equations based on stochastic operational matrix for generalized hat basis functions. J Comput Phys 2014;270:402–15. [32] Heydari MH, Hooshmandasl MR, Cattani C, Ghaini FMM. An efficient computational method for solving nonlinear stochastic itôvolterra integral equations: application for stochastic problems in physics. J Comput Phys 2015;283:148–68. [33] Heydari MH, Hooshmandasl MR, Loghmani GB, Cattani C. Wavelets Galerkin method for solving stochastic heat equation. Int J Comput Math 2016;9. [34] Heydari MH, Hooshmandasl MR, Shakiba A, Cattani C. Legendre wavelets Galerkin method for solving nonlinear stochastic integral equations. Nonlinear Dyn 2016;85(2):1185–202. [35] Khodabin M, Maleknejad K, Rostami M, Nouri M. Interpolation solution in generalized stochastic exponential population growth model. Appl Math Model 2012;36:1023–33. [36] Maleknejad K, Khodabin M, Rostami M. Numerical solutions of stochastic volterra integral equations by a stochastic operational matrix based on block pulse functions. Appl Math Model 2012;55:791–800. [37] Mirzaee F, Hadadiyan E. A collocation technique for solving nonlinear stochastic itô-volterra integral equations. Appl Math Comput 2014;247:1011–20. [38] Mirzaee F, Hadadiyan E. Approximation solution of nonlinear stratonovich volterra integral equations by applying modification of hat functions. J Comput Appl Math 2016;302:272–84. [39] Mohammadi F. A wavelet-based computational method for solving stochastic itô-volterra integral equations. J Comput Phys 2015;298:254–65. [40] Mohammadi F. Efficient Galerkin solution of stochastic fractional differential equations using second kind chebyshev wavelets. Boletim da Sociedade Paranaense de Matemtica 2015;35(1):195–215. [41] Taheri Z, Javadi S, Babolian E. Numerical solution of stochastic fractional integro-differential equation by the spectral collocation method. J Comput Appl Math 2017. doi:10.1016/j.cam.2017.02.027.

124

M.H. Heydari, Z. Avazzadeh and M.R. Mahmoudi / Chaos, Solitons and Fractals 124 (2019) 105–124

[42] Lisei H, Soós A. Approximation of stochastic differential equations driven by fractional Brownian motion. Semin Stoch Anal Random Fields Appl 2008;59:227–41. [43] Mishura Y, Shevchenko G. The rate of convergence for euler approximations of solutions of stochastic differential equations driven by fractional Brownian motion. Stochastic 2008;80(5):489–511. [44] Mirzaee1 F, Hamzeh A. Stochastic operational matrix method for solving stochastic differential equation by a fractional Brownian motion. Int J Appl Comput Math 2017. doi:10.1007/s40819- 017- 0362- 0. [45] Hashemi B, Khodabin M, Maleknejad K, Mediterr K. Numerical solution based on hat functions for solving nonlinear stochastic Itô volterra integral equations driven by fractional Brownian motion. Mediterr J Math 2017;14:24. doi:10. 10 07/s0 0 0 09- 016- 0820- 7. [46] Heydari MH, Mahmoudi MR, Shakiba A, Avazzadeh Z. Chebyshev cardinal wavelets and their application in solving nonlinear stochastic differential equations with fractional Brownian motion. Commun Nonlinear Sci Numer Simulat 2018;64:98–121. [47] Kobelev Y, Kobelev L, Klimontovich Y. Statistical physics of dynamic systems with variable memory. Doklady Phys 2003;48(6):285–9. [48] Sun HG, Chen W, Chen YQ. Variable-order fractional differential operators in anomalous diffusion modeling. Phys A 2009;388(21):4586–92. [49] Krishna PM, Gadre VM, Desai UB. Multifractal based network traffic modeling. New York: Springer; 2003. [50] Sheng H, Chen YQ, Qiu TS. Fractional processes and fractional-order signal processing: techniques and applications. New York: Springer; 2012. [51] Echelard A, Barrière O, Lévy-Véhel J. Terrain modeling with multifractional Brownian motion and self-regulating processes, in computer vision and graphics, lecture notes in computer science. Berlin Heidelberg: Springer; 2010.

[52] Lebovits J, Corlay S, Lévy-Véhel J. Multifractional stochastic volatility models. Math Finance 2014;24(2):364–402. [53] Lim SC, Muniandy SV. Multifractional stochastic volatility models. Phys Rev E 2002;66:021114. [54] Marguez-Lago T, Leier A, Burrage K. Anomalous diffusion and multifractional Brownian motion: simulating molecular crowding and physical obstacles in systems biology. IET Systems Biol 2012;6(4):134–42. [55] Chen Y, Liu L, Sun Y. Numerical solution for the variable order linear cable equation with Bernstein polynomials. Appl Math Comput 2014;238:329–41. [56] Kilicman A, Zhour ZAAA. Kronecker operational matrices for fractional calculus and some applications. Appl Math Comput 2014;238:329–41. [57] Heydari MH, Avazzadeh Z. A new wavelet method for varaible-order fractional optimal control problems. Asian J Control 2018;20(5):1–14. [58] Canuto C, Hussaini M, Quarteroni A, Zang T. Spectral methods in fluid dynamics. Berlin: Springer; 1988. [59] Liu M, Wang K. A note on stability of stochastic logistic equation. Applied Mathematical Letters 2013;26:601–6. [60] Aarató M. nonlinear stochastic equation (Lotka–Volterra model with diffusion). Math. Comput Model 2003;38:709–26. [61] Henderson D, Plaschko P. Stochastic differential equations in science and engineering. Singapore: World Scientific; 2006. [62] Kloeden PE, Platen E. Numerical solution of stochastic differential equations. Berlin: Springer; 1995.