A computational method for solving stochastic Itô–Volterra integral equations based on stochastic operational matrix for generalized hat basis functions

Accepted Manuscript A computational method for solving stochastic Itô-Volterra integral equations based on stochastic operational matrix for generalized hat basis functions M.H. Heydari, M.R. Hooshmandasl, F.M. Maalek Ghaini, C. Cattani

PII: DOI: Reference:

S0021-9991(14)00251-4 10.1016/j.jcp.2014.03.064 YJCPH 5186

To appear in:

Journal of Computational Physics

Received date: 18 September 2013 Revised date: 9 March 2014 Accepted date: 31 March 2014

Please cite this article in press as: M.H. Heydari et al., A computational method for solving stochastic Itô-Volterra integral equations based on stochastic operational matrix for generalized hat basis functions, Journal of Computational Physics (2014), http://dx.doi.org/10.1016/j.jcp.2014.03.064

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A computational method for solving stochastic Itˆo-Volterra integral equations based on stochastic operational matrix for generalized hat basis functions M. H. Heydari1 , M. R. Hooshmandasl2 , F. M. Maalek Ghaini3 , C. Cattani4 1,2,3 1,2,3

Faculty of Mathematics, Yazd University, Yazd, Iran.

The Laboratory of Quantum Information Processing, Yazd University, Yazd, Iran. 4

Department of Mathematics, University of Salerno, Fisciano, Italy e-mail:1 [email protected], 3

[email protected],

4

2

[email protected],

[email protected]

Abstract In this paper, a new computational method based on the generalized hat basis functions is proposed for solving stochastic Itˆ o-Volterra integral equations. In this way, a new stochastic operational matrix for generalized hat functions on the finite interval [0, T ] is obtained. By using these basis functions and their stochastic operational matrix, such problems can be transformed into linear lower triangular systems of algebraic equations which can be directly solved by forward substitution. Also, the rate of convergence of the proposed method is considered and it has been shown that it is O( n12 ). Further in order to show the accuracy and reliability of the proposed method, the new approach is compared with the block pulse functions method by some examples. The obtained results reveal that the proposed method is more accurate and efficient in comparison with the block pule functions method.

Keywords: Generalized hat basis functions; Stochastic operational matrix; Stochastic Itˆo-Volterra integral equations; Itˆo integral; Brownian motion process.

1

Introduction

Stochastic and deterministic functional equations are fundamental for modeling science and engineering phenomena. As the computational power increases, it becomes feasible to use more accurate functional equation models and solve more demanding problems. Moreover the study of stochastic or random functional equations will be very useful in application, because of the fact that they arise in many situations. For example, stochastic integral equations arise in the stochastic formulation of problems in reactor dynamics [1–3], in the study of the growth of biological populations [4], in the theory of automatic systems resulting in delay-differential equations [5], and in many other problems occurring in the general areas of biology, physics and engineering. Also, nowadays, there is an increasing demand to investigate the behavior of even more sophisticated dynamical systems in 1

physical, medical, engineering applications and finance [6–13]. These systems are often dependent on a noise source, like e.g. a Gaussian white noise, governed by certain probability laws, so that modeling such phenomena naturally requires the use of various stochastic differential equations [4, 14–20], or in more complicated cases, stochastic Volterra integral equations and stochastic integro-differential equations [21–28]. In most cases it is difficult to solve such problems explicitly. Therefore, it is necessary to obtain their approximate solutions by using some numerical methods [1, 2, 4–8, 14, 23–25]. Babolian and Mordad [29] have used hat basis function for solving systems of linear and nonlinear integral equations of the second kind by the hat functions. In [30], Tripathi et al. have used generalized hat basis functions to obtain approximate solutions of linear and nonlinear fractional differential equations. In this paper, generalized hat basis functions will be used to solve the following linear stochastic Volterra integral equation:  t  t X(t) = f (t) + K1 (τ, t)X(τ )dτ + K2 (τ, t)X(τ )dB(τ ), t ∈ [0, T ], (1) 0

0

where X(t), f (t), K1 (τ, t) and K2 (τ, t), for t, τ ∈ [0, T ], are some stochastic processes defined on the same probability space (Ω, F, P ), X(t) is the unknown function to be found, B(t) is a Brownian motion process and the second integral in (1) is the Itˆo integral [31, 32]. In order to compute the approximate solution on this equation, we first describe some properties of the generalized hat basis functions. Then, the new operational matrix of stochastic integration for the generalized hat functions is derived and applied to obtain approximate solution for the under study problem. Convergence and error analysis of the proposed method are also investigated and the efficiency of our method is shown on some concrete examples. It is worth noting that, the main advantage of the proposed method is that it reduces the problem under consideration into solving a linear lower triangular system of algebraic equations by expanding the solution in generalized hat functions with unknown coefficients and using the operational matrices of integration and stochastic integration. This paper is organized as follows: In section 2, a brief review of the generalized hat functions and their properties is described. In section 3, the stochastic integration operational matrix of generalized hat functions is obtained. In section 4, the proposed method is described. In section 5, convergence and error analysis of the proposed method are investigated. In section 6, some numerical examples are presented. Finally a conclusion is drawn in section 7, where a discussion on future applications is also given.

2

The generalized hat basis functions and their properties

The traditional hat basis functions are continuous functions, also called triangle, tent or triangular functions are defined on the interval [0, 1]. They are fundamental functions in signal analysis, since their Fourier transform coincide with the square of the sincfunction. Sinc-function is the Fourier transform of the rectangular (or box) function. The translated and dilated instances of the sinc-function can be used to define the so-called 2

Shannon wavelets, which have been shown to be an expedient tool for the solution of integro-differential equations [33–35]. This connection of triangle functions with Shannon wavelets, which are a basis for the L2 (R) functions, it guarantees that the triangle functions can be assumed as basis for the bounded energy functions. The generalized hat functions are extension of traditional hat functions on the finite interval [0, T ]. A set of these basis functions are usually defined on [0, T ] as [30]:  h−t h , 0 ≤ t ≤ h, φ0 (t) = (2) 0, otherwise, ⎧ ⎪ ⎪ ⎨ φi (t) =

⎪ ⎪ ⎩ 

φn (t) =

t−(i−1)h , h

(i − 1)h ≤ t ≤ ih,

(i+1)h−t , h

ih ≤ t ≤ (i + 1)h,

0,

otherwise,

t−(T −h) , h

T − h ≤ t ≤ T,

0,

otherwise,

i = 1, 2, . . . , n − 1,

(3)

(4)

where h = Tn and n is an arbitrary positive integer. Indeed, the interval [0, T ] is divided into n equidistant subintervals [ih, (i + 1)h], i = 0, 1, . . . , n − 1. From the definition of generalized hat functions, we have:  1, i = j, (5) φi (jh) = 0, i = j, and φi (t)φj (t) = 0, |i − j| ≥ 2.

(6)

In addition, these functions form a partition of unity so that

n 

φi (t) = 1.

i=0

An arbitrary function f ∈ L2 [0, T ] can be expanded by generalized hat functions as: f (t) 

n 

fi φi (t) = F T Φ(t) = Φ(t)T F,

(7)

i=0

where F  [f0 , f1 , . . . , fn ]T ,

(8)

Φ(t)  [φ0 (t), φ1 (t), . . . , φn (t)]T .

(9)

The important aspect of using the generalized hat basis functions in approximating a function f (t), lies in the fact that the coefficients fi in (7) are given by: fi = f (ih), i = 0, 1, . . . , n.

(10) 3

Similarly, an arbitrary function of two variables k(x, y) defined on L2 ([0, T ] × [0, T ]) can be expanded by generalized hat basis functions as: k(τ, t) = Φ(τ )T ΛΨ(t),

(11)

where Φ(τ ) and Ψ(t) are (n+1) and (m + 1) dimensional generalized hat function vectors, respectively, and Λ is the (n + 1) × (m + 1) generalized hat functions coefficients matrix with entries aij , i = 0, 1, . . . , n, j = 0, 1, . . . , m, as follows: aij = k(ih, js), T where s = m . In this paper, for convenience, we put n = m, so we From relation (6), we have: ⎛ φ0 (t)φ1 (t) φ20 (t) ⎜ ⎜ φ0 (t)φ1 (t) φ21 (t) φ1 (t)φ2 (t) ⎜ ⎜ .. .. ⎜ . . ⎜ T Φ(t)Φ(t) = ⎜ . .. ⎜ ⎜ ⎜ ⎜ ⎝

have aij = k(ih, jh). ⎞ ⎟ ⎟ ⎟ ⎟ .. ⎟ . ⎟ ⎟. .. .. ⎟ . . ⎟ ⎟ .. .. . . φn−1 (t)φn (t) ⎟ ⎠ 2 φn (t) φn−1 (t)φn (t) (12)

By considering (5) and expanding have: ⎛ φ0 (t) 0 ⎜ ⎜ 0 φ1 (t) Φ(t)Φ(t)T  ⎜ ⎜ .. .. ⎝ . . 0 0

entries of Φ(t)Φ(t)T by generalized hat functions, we ...

⎞

0

... 0 .. .. . . . . . φn (t)

⎟ ⎟ ⎟. ⎟ ⎠

Integration of vector Φ(t) defined in (9) can be expressed as [30]:  t Φ(τ )dτ  P Φ(t), t ∈ [0, T ],

(13)

(14)

0

where P is an (n + 1) × (n + 1) generalized operational matrix for integration and is given by: ⎞ ⎛ 0 1 1 ... 1 1 ⎜ 0 1 2 ... 2 2 ⎟ ⎟ ⎜ ⎟ h⎜ ⎜ 0 0 1 ... 2 2 ⎟ (15) P = ⎜ . . . . ⎟. . . 2 ⎜ .. .. .. . . .. .. ⎟ ⎟ ⎜ ⎝ 0 0 0 0 1 2 ⎠ 0 0 0 0 0 1 4

3

The stochastic integration operational matrix

Theorem 3.1. Let Φ(t) be the vector defined in (9). The Itˆ o integral of Φ(t) can be expressed as:  t Φ(τ )dB(τ )  Ps Φ(t), (16) 0

where the (n + 1) × (n + 1) stochastic operational matrix of integration is given by: ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ Ps = ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

0

α0 (h)

0 B(h) + α1 (h)

α0 (h)

0

...

α0 (h)

α0 (h)

β1 (h)

0

...

β1 (h)

β1 (h)

B(2h) + α2 (h) β2 (h) . . .

β2 (h)

β2 (h)

.. .

.. .

0

0

.. .

.. .

.. .

.. .

..

0

0

0

0

. . . B ((n − 1)h) + αn−1 (h)

0

0

0

0

...

.

0

βn−1 (h) B(T ) + αn (h) (17)

and ⎧  1 h ⎪ ⎪ α0 (h) = B(τ )dτ, ⎪ ⎪ ⎪ h 0 ⎪ ⎪ ⎪  ⎪ ⎨ 1 ih B(τ )dτ, i = 1, 2, . . . , n, αi (h) = − h (i−1)h ⎪ ⎪ ⎪   ⎪  (i+1)h ⎪ ih ⎪ 1 ⎪ ⎪ ⎪ B(τ )dτ − B(τ )dτ , i = 1, 2, . . . , n − 1. ⎩ βi (h) = − h (i−1)h ih

(18)

Proof. By considering definitions of φi (t) (i = 0, 1, . . . , n), and integration by parts, we have:  t  t φ0 (τ )dB(τ ) = φ0 (t)B(t) − φ0 (0)B(0) − φ0 (τ )B(τ )dτ, (19) 0

0

Also from (2), we obtain: 1 φ0 (t) = − χ[0,h] (t), h

(20)

where χ[0,h] (t) is the characteristic function of [0, h]. Expanding (19) in terms of generalized hat basis functions yields: 

t 0

φ0 (τ )dB(τ ) 

n 

a0j φj (t),

(21)

j=0

5

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

where the coefficient a0j is the value of Then from (19), we have:  a0j =

t 0

φ0 (τ )dB(τ ) at jth node point (jh).

jh 0

φ0 (τ )dB(τ ) = φ0 (jh)B(jh) +

1 h



jh 0

χ[0,h] (τ )B(τ )dτ, j = 0, 1, . . . , n. (22)

Now, using (22) and definitions of φ0 (t) and χ[0,h] (t), we obtain: ⎧ ⎪ ⎨ a00 = 0,  1 h ⎪ ⎩ a0j = B(τ )dτ, j = 1, 2, . . . , n. h 0

(23)

The Itˆo integral of φn (t) yields: 



t 0

φn (τ )dB(τ ) =

t

T −h

 φn (τ )dB(τ ) = φn (t)B(t)−φn (T −h)B(T −h)−

t T −h

φn (τ )B(τ )dτ. (24)

By considering (4) and expanding (24) in terms of generalized hat basis functions, we also have:  t  t n  φn (τ )dB(τ ) = φn (τ )dB(τ )  anj φj (t), (25) 0

T −h

j=0

where  anj =

jh T −h

φn (τ )dB(τ ) = φn (jh)B(jh) −

1 h



jh T −h

χ[T −h,T ] (τ )B(τ )dτ.

(26)

Now, using definitions of φn (t), χ[T −h,T ] (t) and relation (26), we obtain: ⎧ ⎪ ⎨ anj = 0, ⎪ ⎩ ann

j = 0, 1, . . . , n − 1,  1 T = B(T ) − B(τ )dτ. h T −h

(27)

Also, by considering definition of φi (t) for i = 1, 2, . . . , n − 1, we have: 



t 0

φi (τ )dB(τ ) =

t (i−1)h

 φi (τ )dB(τ ) = φi (t)B(t)−φi ((i − 1)h) B ((i − 1)h)−

t (i−1)h

φi (τ )B(τ )dτ. (28)

By considering (3) and expanding (28) in terms of generalized hat basis functions, we have:  t  t n  φi (τ )dB(τ ) = φi (τ )dB(τ )  aij φj (t), i = 1, 2, . . . , n − 1, (29) 0

(i−1)h

j=0

6

t where the coefficient aij is the value of (i−1)h φi (τ )dB(τ ) at jth node point (jh). Then from (28), we have:  1 jh δi (τ )B(τ )dτ, aij = φi (jh)B(jh) − h (i−1)h

(30)

where δi (t) = χ[(i−1)h,ih] (t) − χ(ih,(i+1)h] (t), i = 1, 2, . . . , n − 1. Now, using definitions of φi (t), δi (t) (i = 1, 2, . . . , n − 1), ⎧ aij = 0, ⎪ ⎪ ⎪ ⎪  ih ⎪ ⎪ ⎪ ⎨ aii = B(ih) − 1 B(τ )dτ, h (i−1)h ⎪   ⎪  (i+1)h ⎪ ih ⎪ 1 ⎪ ⎪ ⎪ B(τ )dτ − B(τ )dτ , ⎩ aij = − h (i−1)h ih

(31)

and (30), we obtain: j < i, j = i,

(32)

j > i,

which completes the proof.

4

The proposed method using stochastic operational matrix

In this section, we apply the stochastic operational matrix for generalized hat basis functions to solve linear stochastic Volterra integral equation:  t  t K1 (τ, t)X(τ )dτ + K2 (τ, t)X(τ )dB(τ ), t ∈ [0, T ], (33) X(t) = f (t) + 0

0

where X(t), f (t), K1 (τ, t) and K2 (τ, t), for τ, t ∈ [0, T ], are the stochastic processes defined on the same probability space (Ω, F, P ), and X(t) is an unknown function to be found. Also B(t) is a Brownian motion process and the second integral in (33) is an Itˆo integral. For solving this problem, we approximate X(t), f (t), K1 (τ, t) and K2 (τ, t) as follows: X(t)  X T Φ(t) = Φ(t)T X,

(34)

f (t)  F T Φ(t) = Φ(t)T F,

(35)

K1 (τ, t)  Φ(τ )T K1 Φ(t) = Φ(t)T K1T Φ(τ ), T

K2 (τ, t)  Φ(τ ) K2 Φ(t) = Φ(t)

T

(36)

K2T Φ(τ ),

(37)

where X and F are stochastic generalized hat coefficient vectors, and K1 and K2 are stochastic generalized hat coefficient matrices. Substituting (34)-(37) into (33) yields: T

T

X Φ(t)  F Φ(t)+X

T



t 0

T

Φ(τ )Φ(τ ) dτ

 K1 Φ(t)+X

T



t 0

T



Φ(τ )Φ(τ ) dB(τ ) K2 Φ(t). (38)

7

Now let K1i , K2i , Ri and Rsi respectively be the ith rows of the matrices K1 , K2 , P and Ps , DK i be a diagonal matrix with K1i as its diagonal entries and DK i be a diagonal matrix 1 2 with K2i as its diagonal entries. By the previous assumptions, we have: ⎛ 

t 0

Φ(τ )Φ(τ )T dτ



⎜ ⎜ ⎜ K1 Φ(t)  ⎜ ⎜ ⎜ ⎝

R1 Φ(t)K11 Φ(t) R2 Φ(t)K12 Φ(t) .. . Rn+1 Φ(t)K1n+1 Φ(t)

⎞

⎛

⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎠ ⎝

⎞

R1 Dk11 R2 Dk12 .. . Rn+1 Dkn+1

⎟ ⎟ ⎟ ⎟ Φ(t) = B1 Φ(t), ⎟ ⎟ ⎠

1

(39) where

⎛

⎜ ⎜ ⎜ h⎜ ⎜ B1 = ⎜ 2⎜ ⎜ ⎜ ⎝

1 0 k01

1 k02

...

1 k0n

0

...

1 knn

⎞

⎟ 1 1 1 ⎟ 0 k11 2k12 . . . 2k1n ⎟ ⎟ 1 1 0 0 k22 . . . 2k2n ⎟ ⎟. ⎟ .. .. .. .. ⎟ .. . . . . . ⎟ ⎠ 0

0

Also, we have: ⎛ 

t 0

⎜ ⎜ ⎜ T Φ(τ )Φ(τ ) dB(τ ) K2 Φ(t)  ⎜ ⎜ ⎜ ⎝ 

Rs1 Φ(t)K21 Φ(t) Rs2 Φ(t)K22 Φ(t) .. . Rsn+1 Φ(t)K2n+1 Φ(t)

⎞

⎛

⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎠ ⎝

Rs1 Dk21 Rs2 Dk22 .. . Rsn+1 Dkn+1

⎞ ⎟ ⎟ ⎟ ⎟ Φ(t) = B2 Φ(t), ⎟ ⎟ ⎠

2

(40) where

⎛

⎜ ⎜ ⎜ ⎜ ⎜ B2 = ⎜ ⎜ ⎜ ⎜ ⎝

0

2 α0 (h)k01

2 0 (B(h) + α1 (h)) k11

2 α0 (h)k02

...

2 α0 (h)k0n

2 β1 (h)k12

...

2 β1 (h)k1n

2 (B(2h) + α2 (h)) k22 ...

2 β2 (h)k2n

0

0

.. .

.. .

.. .

..

0

0

0

2 . . . (B(T ) + αn ) knn

.

.. .

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎠

Now, by substituting (39) and (40) into (33), we obtain: X T Φ(t)  F T Φ(t) + X T B1 Φ(t) + X T B2 Φ(t).

8

(41)

This equation yields: X T (I − B1 − B2 )  F T .

(42)

So, by putting M = (I − B1 − B2 )T and replacing  by =, we obtain the following system of algebraic equations: M X = F.

(43)

Finally by solving this system and determining X, we obtain the approximate solution of the problem.

5

Error analysis

In this section, we investigate the error analysis of the proposed method. From now on,  denotes the sup-norm which for any continuous function f (t) is defined on the interval [0, T ] by: f (t) = sup |f (t)|.

(44)

t∈[0,T ]

Theorem 5.1. Suppose f (t) ∈ C 2 ([0, T ]) and en (t) = f (t) − fn (t), t ∈ I = [0, T ], where n  fn (t) = f (ih)φi (t) is the generalized hat functions expansion of f (t). Then we have: i=0

en (t) ≤

T 2  f (t), 2n2

(45)

and so the convergence is of order two, that is   1 . en (t) = O n2

(46)

Proof. Let  enj (t) =

f (t) − fn (t), t ∈ Dj ,

(47)

t ∈ I − Dj ,

0,

T where Dj = {t|jh ≤ t < (j + 1)h, h = }, j = 0, 1, . . . , n − 1. n Then we have: enj (t) = f (t) −

n 

f (ih)φi (t) = f (t) − (f (jh)φj (t) + f ((j + 1)h)φj+1 (t))

i=0

 = f (t) − f (jh)



(j + 1)h − t h



 + f (jh + h) 9

t − jh h







= f (t) − f (jh) + (t − jh)

f (jh + h) − f (jh) h

 .

(48)

As h → 0, we get:   enj (t)  f (t) − f (jh) + (t − jh)f  (jh) .

(49)

Expanding f (t) in the Taylor’s series, in terms of (t − jh) yields: f (t) = f (jh) + (t − jh)f  (jh) +

(t − jh)2  f (ηj ), 2!

(50)

where ηj ∈ (jh, (j + 1)h). Now, from (49) and (50), we have: enj (t) =

(t − jh)2  f (ηj ). 2!

(51)

Since |t − jh| < h and nh = T , then we have: |enj (t)| ≤

T 2  |f (ηj )|. 2n2

(52)

Therefore en (t) = sup |enj (t)| ≤ t∈Dj

T2 T2 sup |f  (t)| = 2 f  (t). 2 2n t∈I 2n

Also from (53), we have:   1 en (t) = O . n2

(53)

(54)

This completes the proof. Theorem 5.2. Suppose f (τ, t) ∈ C 2 ([0, T ] × [0, T ]) and en (τ, t) = f (τ, t)−fn (τ, t), (τ, t) ∈ n n   D = [0, T ] × [0, T ], where fn (τ, t) = f (lh, rh)φl (τ )φr (t) is the generalized hat funcl=0 r=0

tions expansion of f (τ, t). Then we have:   2     2  ∂ f (τ, t)   ∂ 2 f (τ, t)  ∂ f (τ, t)  T2      ,   + 2 + e(τ, t) ≤ 2  2n ∂τ 2  ∂τ ∂t   ∂t2  and so the convergence is of order two, that is   1 en (τ, t) = O . n2

10

(55)

(56)

Proof. Let  enij (τ, t) =

f (τ, t) − fn (τ, t), (τ, t) ∈ Dij ,

(57)

(τ, t) ∈ D − Dij ,

0,

T where Dij = {(τ, t)|ih ≤ s < (i + 1)h, jh ≤ t < (j + 1)h, h = }, i, j = 0, 1, . . . , n − 1. n Then for i, j = 0, 1, . . . , n − 1, we have: enij (τ, t) = f (τ, t) −

n n  

f (lh, rh)φl (τ )φr (t)

l=0 r=0

= f (τ, t) − (f (ih, jh)φi (τ )φj (t) + f (ih, (j + 1)h)φi (τ )φj+1 (t) +f ((i + 1)h, jh)φi+1 (τ )φj (t) + f ((i + 1)h, (j + 1)h)φi+1 (τ )φj+1 (t)) .

(58)

Now, by considering two dimensional Taylor’s series of f (τ, t) in terms of (τ − ih) and (t − jh), as h → 0, we have the following error estimation: 1 enij (τ, t) = f (τ, t) − fn (τ, t) = 2!



∂ ∂ (τ − ih) + (t − jh) ∂τ ∂t

2 f (τ, t)|(ηi ,ξj ) ,

where ηi ∈ (ih, (i + 1)h), ξj ∈ (jh, (i + 1)h). Now, since |τ − ih|, |t − jh| < h and T = nh, we have:   2      ∂ f (ηi , ξj )   ∂ 2 f (ηi , ξj )  T 2  ∂ 2 f (ηi , ξj )    .   |enij (τ, t)| ≤ 2   + 2  ∂τ ∂t  +   2n ∂τ 2 ∂t2

(59)

(60)

Therefore   2     2  ∂ f (τ, t)   ∂ f (τ, t)   ∂ 2 f (τ, t)  T2     ,   en (τ, t) = sup |enij (τ, t)| ≤ 2 sup  + 2 + 2n (τ,t)∈D ∂τ 2  ∂τ ∂t   ∂t2  (τ,t)∈Dij (61) and consequently e(τ, t) ≤

T2 2n2

  2    2  2  ∂ f (τ, t)        + 2  ∂ f (τ, t)  +  ∂ f (τ, t)  .  ∂τ 2   ∂τ ∂t   ∂t2 

So from (62), we have:   1 . en (τ, t) = O n2

(62)

(63)

This completes the proof. Theorem 5.3. Suppose X(t) be the exact solution of (33) and Xn (t) be the generalized hat series approximate solution of (33) and also assume that: (H1) X(t) ≤ ρ, t ∈ I = [0, T ], 11

(H2) K1 (τ, t) ≤ M1 , (τ, t) ∈ D = I × I, (H3) K2 (τ, t) ≤ M2 , (τ, t) ∈ D = I × I, (H4) T (M1 + λ(h)) + ργ(h)B(t) < 1, then Γ(h) + T ρλ(h) + ργ(h)B(t) X(t) − Xn (t) ≤ , 1 − (T (M1 + λ(h)) + ργ(h)B(t))

(64)

where

h2  f (t), 2    2    2 2      h2   ∂ K1 (τ, t)  + 2  ∂ K1 (τ, t)  +  ∂ K1 (τ, t)  , λ(h) =       2 2 2  ∂τ   2 ∂τ ∂t   2 ∂t  2 2      h   ∂ K2 (τ, t)  + 2  ∂ K2 (τ, t)  +  ∂ K2 (τ, t)  . γ(h) =       2 2 2 ∂τ ∂τ ∂t ∂t Γ(h) =

Proof. From (33), we have: X(t) − Xn (t) = f (t) − fn (t) + 

t

+ 0



t 0

(K1 (τ, t)X(τ ) − K1n (τ, t)Xn (τ )) dτ

(K2 (τ, t)X(τ ) − K2n (τ, t)Xn (τ )) dB(τ ).

(65)

Now, we can write X(t) − Xn (t) ≤ f (t) − fn (t) + tK1 (τ, t)X(τ ) − K1n (τ, t)Xn (τ ) +B(t)K2 (τ, t)X(τ ) − K2n (τ, t)Xn (τ ).

(66)

Then by using theorems 5.1 and 5.2, and assumptions (H1) and (H2), we have: K1 (τ, t)X(τ ) − K1n (τ, t)Xn (τ ) ≤ K1 (τ, t)X(t) − Xn (t) +K1 (τ, t) − K1n (τ, t) (X(t) − Xn (t) + X(t)) ≤ (M1 + λ(h))X(t) − Xn (t) + ρλ(h),

(67)

and K2 (τ, t)X(τ ) − K2n (τ, t)Xn (τ ) ≤ K2 (τ, t)X(t) − Xn (t) +K2 (τ, t) − K2n (τ, t) (X(t) − Xn (t) + X(t)) ≤ (M2 + γ(h))X(t) − Xn (t) + ργ(h).

(68)

So X(t) − Xn (t) ≤ Γ(h) + t ((M1 + λ(h))X(t) − Xn (t) + ρλ(h)) +B(t)(M2 + γ(h))X(t) − Xn (t) + ργ(h)

(69)

Γ(h) + T ρλ(h) + ργ(h)B(t) . 1 − (T (M1 + λ(h)) + ργ(h)B(t))

(70)

or X(t) − Xn (t) ≤

Also from (70), we have X(t) − Xn (t) = O( n12 ). This completes the proof. 12

6

Numerical examples

In this section, we consider some numerical examples which their exact solutions are available to illustrate the efficiency and reliability of the proposed method. Note from now on, n is the number of basis functions and m is the number of iterations. Example 1. Let us consider the following linear stochastic Volterra integral equation [6]:  X(t) = 1 +

t 0

τ 2 X(τ )dτ +



t 0

(τ, t) ∈ [0, 0.5],

τ X(τ )dB(τ ),

(71)

where X(t) is a known stochastic process defined on the probability space (Ω, F, P ), and B(t) is a Brownian motion process with the exact solution   3  t t τ dB(τ ) . + X(t) = exp 6 0 A comparison between the unmetrical solutions given by the proposed method (GHM) for different values of n and m via the block pulse method (BPM) are shown in Figs. 1 and 2 and Table 1. From Figs. 1 and 2 it can be seen that the results obtained by

Exact and approximate solutions

1.2

Exact solution Proposed method Block pulse method

1.15

1.1

1.05

1

0.95 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

t

Fig. 1: The graphs of the approximate and exact solutions of Example 1 for n=32, m=25.

the proposed method are more accurate than the block pulse method. It is worth noting that the needed computational work in constructing matrices X, F, K1 and K2 for the proposed method is less than the block pulse method.

13

Exact and approximate solutions

1.05

1

0.95

0.9 Exact solution Proposed method Block pulse method 0.85 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

t

Fig. 2: The graphs of the approximate and exact solutions of Example 1 for n=64, m=75.

t 0.1 0.2 0.3 0.4 0.5

Table 1: Numerical results n = 32, m = 25 BPM GHM 0.9664232 0.9684280 1.0071571 1.0066520 1.0462828 1.0460315 1.0882984 1.0823571 1.0539179 1.0118890

of Example 1 for different values of n and n = 64, m = 75 Exact BPM GHM 0.9669176 1.0074152 1.0075817 1.0105371 0.9714866 0.9739018 1.0378587 0.9403179 0.9311638 1.0760953 0.9648926 0.9611705 0.9907228 1.0261553 1.0148483

Example 2. Consider the linear stochastic Volterra integral equation [6]:  t  t 1 cos(τ )X(τ )dτ + sin(τ )X(τ )dB(τ ), (τ, t) ∈ [0, 0.5], X(t) = + 12 0 0

m. Exact 1.0064729 0.9726610 0.9403846 0.9634115 0.9941971

(72)

where X(t) is a known stochastic process defined on the probability space (Ω, F, P ), and B(t) is a Brownian motion process with the exact solution    t 1 t sin(2t) sin(τ )dB(τ ) . (73) X(t) = exp − + sin(t) + + 12 4 8 0 A comparison between the unmetrical solutions given by the proposed method for some different values of n and m via the block pulse method are shown in Table 2 and Figs. 3 and 4. From Figs. 3 and 4, it can be seen that the proposed method provides an acceptable solution for this problem. 14

0.16

Exact and approximate solutions

0.15

Exact solution Proposed method Block pulse method

0.14 0.13 0.12 0.11 0.1 0.09 0.08 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

t

Fig. 3: The graphs of the approximate and exact solutions of Example 2 for n=32, m=25.

0.15

Exact and approximate solutions

0.14

Exact solution Proposed method Block pulse method

0.13 0.12 0.11 0.1 0.09 0.08 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

t

Fig. 4: The graphs of the approximate and exact solutions of Example 2 for n=64, m=75.

7

Applications in the mathematical finance

The market consists of a riskless cash bond, {B(t)}t≥0 , and a single risky asset with price process {S(t)}t≥0 governed by [6]: 15

t 0.1 0.2 0.3 0.4 0.5

Table 2: Numerical results n = 32, m = 25 BPM GHM 0.0890981 0.0891427 0.1016352 0.1020402 0.1165245 0.1160338 0.1307434 0.1301881 0.1355362 0.1309717



dB(t) = r(t)B(t)dt,

of Example 2 for different values of n and n = 64, m = 75 Exact BPM GHM 0.0890088 0.0925190 0.0927489 0.1024286 0.0984359 0.0987536 0.1151393 0.1044737 0.1033985 0.1294754 0.1164332 0.1169554 0.1283848 0.1320381 0.1309579

B0 = 1,

dS(t) = μ(t)S(t)dt + σ(t)S(t)dW (t), or

⎧ ⎪ ⎨ dB(t) = r(t)B(t)dt, B0 = 1,  t  t ⎪ ⎩ S(t) = S0 + μ(τ )S(τ )dτ + σ(τ )S(τ )dW (τ ), 0

m. Exact 0.0926465 0.0986228 0.1044119 0.1160552 0.1284231

(74)

(75)

0

where {W (t)}t≥0 is a P -Brownian motion generating the filtration {F(t)}t≥0 and r(t), μ(t) and σ(t) are {F(t)}t≥0 -predictable processes [6, 36]. A solution to these equations should take the form [6]:  t  B(t) = exp r(τ )dτ , (76) 0

 t  S(t) = S0 exp

0

   t 1 2 μ(τ ) − σs dτ + σ(τ )dW (τ ) . 2 0

For example, consider the following general stock model [6]: ⎧ ⎪ ⎨ dB(t) = sin(t)B(t)dt, B0 = 1, t ∈ [0, 0.8]  t  t 1 ⎪ ⎩ S(t) = ln(1 + τ )S(τ )dτ + τ S(τ )dW (τ ), τ, t ∈ [0, 0.8], + 10 0 0

(77)

(78)

with B(t) = exp (1 − cos(t)) and the exact solution    t t3 1 τ dW (τ ) . exp (1 + t) ln(1 + t) − t − + S(t) = 10 6 0 A comparison between the unmetrical solutions given by the proposed method for some different values of n and m via the block pulse method are shown in Table 3 and Figs. 5 and 6. From Figs. 5 and 6 it can be seen that the results obtained by the proposed method are more accurate than the block pulse method, because of the fact that the rate of convergence of the proposed method is more than the block pulse method. 16

0.17

Exact and approximate solutions

0.16

Exact solution Proposed method Block pulse method

0.15 0.14 0.13 0.12 0.11 0.1 0.09 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

t

Fig. 5: The graphs of the approximate and exact solutions of the above example for n=32, m=25.

0.13

Exact and approximate solutions

0.125

Exact solution Proposed method Block pulse method

0.12 0.115 0.11 0.105 0.1 0.095 0.09 0.085 0.08 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

t

Fig. 6: The graphs of the approximate and exact solutions of the above example for n=64, m=75.

8

Conclusion

Some stochastic differential equations can be written as stochastic Volterra integral equations given in (1). It may be impossible to find the exact solutions of such problems. So, 17

Table 3: Numerical results of the n = 32, m = 25 t BPM GHM 0.1 0.1947044 0.0970802 0.2 0.1893196 0.0941046 0.3 0.1055400 0.1036652 0.4 0.1107672 0.1099689 0.5 0.1229852 0.1211484 0.6 0.1349969 0.1315393 0.7 0.1477287 0.1406459 0.8 0.1320414 0.1218842

above stock model for different values of n and m. n = 64, m = 75 Exact BPM GHM Exact 0.0970619 0.1008782 0.1003858 0.1003672 0.0939736 0.1047496 0.1043171 0.1041724 0.1031824 0.0946055 0.0954179 0.0949809 0.1087882 0.0835175 0.0857424 0.0848358 0.1186166 0.0849324 0.0874110 0.0856135 0.1268358 0.1222197 0.1208435 0.1166484 0.1328019 0.1138578 0.1152089 0.1089121 0.1118127 0.1253731 0.1225266 0.1126076

it would be convenient to determine their numerical solutions using a stochastic numerical analysis. Employing generalized hat functions as basis functions to solve linear stochastic Volterra integral equations is very simple and effective. In this paper, the stochastic operational matrix of Itˆ o-integration for generalized hat functions was derived and applied for solving linear stochastic Volterra integral equations. The convergence and error analysis of the proposed method were investigated. Applicability and accuracy of the proposed method were checked on some examples. Moreover, the results of the proposed method were in good agreement with those obtained by using the block pulse functions method.

References [1] J. J. Levin and J. A. Nohel, “On a systemo f integro-differenetqiaula tionso ccurringin reactor dynamics,” J. Math. Mech., vol. 9, pp. 347–368, 1960. [2] R. K. Miller, “On a system of integro-differential equations occurring in reactor dynamics,” SIAM J. Appl. Math., vol. 14, pp. 446–452, 1966. [3] P. A. Cioica, S. Dahlke, “Spatial besov regularity for semilinear stochastic partial differential equations on bounded lipschitz domains,” International Journal of Computer Mathematics, vol. 89(18), pp. 2443–2459, 2012. [4] M. Khodabin, K. Maleknejad, M. Rostami, and M. Nouri, “Interpolation solution in generalized stochastic exponential population growth model,” Applied Mathematical Modelling, vol. 36, pp. 1023–1033, 2012. [5] M. N. Oˇguzt¨oreli, Time-Lag control systems, New York: Academic Press, 1966. [6] K. Maleknejad, M. Khodabin, and M. Rostami, “Numerical solutions of stochastic volterra integral equations by a stochastic operational matrix based on bloch pulse functions,” Mathematical and Computer Modelling, vol. 55, pp. 791–800, 2012. [7] M. Khodabin, K. Maleknejad, M. Rostami, and M. Nouri, “Numerical approach for solving stochastic volterra-fredholm integral equations by stochastic operational matrix,” Computre and Mathematics with Applications, vol. 64, pp. 1903–1913, 2012. 18

[8] K. Maleknejad, M. Khodabin, and M. Rostami, “A numerical method for solving mdimensional stochastic itˆ o-volterra integral equations by stochastic operational matrix,” Computre and Mathematics with Applications, vol. 63, pp. 133–143, 2012. [9] Y. Cao, D. Gillespie, L. Petzold, ”Adaptive explicit-implicit tau-leaping method with automatic tau selection”, J. Chem. Phys., vol. 126, pp. 224101, 2007. [10] P. Ru´e, J. Vill` a-Freixa, K. Burrage, ”Simulation methods with extended stability for stiff biochemical Kinetics”, BMC Syst. Biol., vol. 4(1), pp. 110, 2010. [11] K.D. Elworthy, A. Truman, H.Z. Zhao, J.G. Gaines, ”Approximate traveling waves for generalized KPP equations and classical mechanics”, Proc. Roy. Soc. London Ser. A, vol. 446, pp. 529–554, 1994. [12] E. Platen, N. Bruti-Liberati, Numerical solution of stochastic differential equations with jumps in Finance, Springer-Verlag, 2010. [13] M. Ehler, “Shrinkage rules for variational minimization problems and applications to analytical ultracentrifugation,” Journal of Inverse and Ill-posed Problems, vol. 19(45), pp. 593–614, 2011. [14] M. Khodabin, K. Maleknejad, M. Rostami, and M. Nouri, “Numerical solution of stochastic differential equations by second order rungekutta methods,” Applied Mathematical Modelling, vol. 53, pp. 1910–1920, 2011. [15] P. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations. Berlin: Springer-Verlag, 1999. [16] J. C. Cortes, L. Jodar, and L. Villafuerte, “Numerical solution of random differential equations: a mean square approach,” Mathematical and Computer Modelling, vol. 45, pp. 757–765, 2007. [17] J. C. Cortes, L. Jodar, and L. Villafuerte, “Mean square numerical solution of random differential equations: facts and possibilities,” Computers and Mathematics with Applications, vol. 53, pp. 1098–1106, 2007. [18] B. Oksendal, Stochastic Differential Equations, fifth ed., in: An Introduction with Applications. New York: Springer-Verlag, 1998. [19] H. Holden, B. Oksendal, J. Uboe, and T. Zhang, tochastic Partial Differential Equations, second ed. New York: Springer-Verlag, 2009. [20] A. Abdulle, A. Blumenthal, ”Stabilized multilevel Monte Carlo method for stiff stochastic differential equations”, Journal of Computational Physics, vol. 251, (2013), 445–460. [21] M. Berger and V. Mizel, “Volterra equations with itˆo integrals i,” Journal of Integral Equations, vol. 2, pp. 187–245, 1980.

19

[22] M. Murge and B. Pachpatte, “On second order ito type stochastic integrod-ifferential equations, analele stiintifice ale universitatii,” I. Cuzadin Iasi, Mathematica, pp. 25– 34, 1984. [23] M. Murge and B. Pachpatte, “Succesive approximations for solutions of second order stochastic integro-differential equations of ito type,” Indian Journal of Pure and Applied Mathematics, vol. 21 (3), pp. 260–274, 1990. [24] X. Zhang, “Euler schemes and large deviations for stochastic volterra equations with singular kernels,” Journal of Differential Equations, vol. 244, pp. 2226–2250, 2008. [25] S. Jankovic and D. Ilic, “One linear analytic approximation for stochastic integrodifferential eauations,” Acta Mathematica Scientia, vol. 308 (4), pp. 1073–1085, 2010. [26] X. Zhang, “Stochastic volterra equations in banach spaces and stochastic partial differential equation,” Acta Journal of Functional Analysis, vol. 258, pp. 1361–1425, 2010. [27] J. Yong, “Backward stochastic volterra integral equations and some related problems,” Stochastic Processes and their Applications, vol. 116, pp. 779–795, 2006. [28] A. Djordjevic, S. Jankovic, “On a class of backward stochastic Volterra integral equations”, Applied Mathematics Letters, doi.org/10.1016/j.aml.2013.07.006. [29] E. Babolian and M. Mordad, “A numerical method for solving systems of linear and nonlinear integral equations of the second kind by hat basis function,” Computers and mathematics with applications, vol. 62, pp. 187–198, 2011. [30] M. P. Tripathi, V. K. Baranwal, R. K. Pandey, and O. P. Singh, “A new numerical algorithm to solve fractional differential equations based on operational matrix of generalized hat functions,” Commun Nonlinear Sci Numer Simulat, vol. 18, pp. 1327– 1340, 2013. [31] B. K. Oksendal, Stochastic Differential Equations: An Introduction with Applications, 4th ed., Springer, 1995. [32] L. Arnold, Stochastic Differential Equations: Theory and Applications, Wiley, 1974. [33] C.Cattani, “Shannon Wavelets Theory”, Mathematical Problems in Engineering, vol. 2008, pp. 1–24, 2008. [34] C.Cattani, “Shannon Wavelets for the solution of Integro-Differential Equations”, Mathematical Problems in Engineering, vol. 2010, pp. 1–22, 2010. [35] C. Cattani, “Fractional calculus and Shannon wavelets”, Mathematical Problems in Engineering, vol. 2012, pp. 1–26, 2012. [36] A. Etheridge, A course in Financial Calculus, Cambridge University Press, 2002.

20