Numerical solution of stochastic integral equations by using Bernoulli operational matrix

Available online at www.sciencedirect.com

ScienceDirect Mathematics and Computers in Simulation 165 (2019) 238–254 www.elsevier.com/locate/matcom

Original Articles

Numerical solution of stochastic integral equations by using Bernoulli operational matrix Rebiha Zeghdane Department of Mathematics, Faculty of Mathematics and Informatics, University of Bordj-Bou-Arreridj, Algeria Received 28 March 2018; received in revised form 7 March 2019; accepted 13 March 2019 Available online 1 April 2019

Abstract In this paper, a new computational method based on stochastic operational matrix for integration of Bernoulli polynomials is proposed for solving nonlinear Volterra–Fredholm–Hammerstein stochastic integral equations. By using this new operational matrix of integration and the so-called collocation method, nonlinear Volterra–Fredholm–Hammerstein stochastic integral equation is reduced to nonlinear system of algebraic equations with unknown Bernoulli coefficients. This work is inspired by Bazm (2015), where the authors study the deterministic integral equations. In order to show the rate of convergence of the suggested approach, we present theorems on convergence analysis and error estimation. Some illustrative error estimations and examples are provided and included to demonstrate applicability and accuracy of the technique. c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights ⃝ reserved. Keywords: Bernoulli polynomials; Stochastic operational matrix; Itô integral; Collocation method; Numerical solution

1. Introduction Stochastic and deterministic integral equations are fundamental for modeling science and engineering phenomena. For Fredholm–Hammerstein integral equations, the classical method of successive approximations was introduced in [16]. In [11,12], the authors applied a collocation type method to solve Volterra–Hammerstein integral equations. The Fredholm–Volterra–Hammerstein type of equations appear in many problems of mathematical physics, theory of elasticity contacts problems [1,24]. 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. These equations 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 [5,13,14,21], or in more complicated cases, stochastic integral equations and stochastic integrodifferential equations. 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 [6,9–11,14,17–23,29]. The available set of orthogonal functions can be divided into three classes, piecewise constant basis functions, set of orthogonal E-mail address: [email protected]. https://doi.org/10.1016/j.matcom.2019.03.005 c 2019 International Association for Mathematics and Computers in Simulation (IMACS). Published by Elsevier B.V. All rights 0378-4754/⃝ reserved.

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

239

polynomials and the third one, the set of sine–cosine functions in Fourier series. There are some works that used Bernoulli polynomials as basis for numerically solving integral equations such as [30,31]. In this work, some operational matrix for integration and differentiation of Bernoulli polynomials is used to convert the considered problems to systems of algebraic equations. One of the advantages of Bernoulli basis functions, in comparison with the basis used in [7,20], is that they provide more accurate approximations of the problem solution with a fewer number of basis functions. On the other hand, the implementation of Bernoulli operational matrix method is very simple. In this paper, we consider the numerical solution of the second kind stochastic Volterra–Fredholm–Hammerstein integral equation. Bernoulli polynomials basis functions will be used to solve the following stochastic integral equation. ∫ x ∫ 1 ∫ x u(x) = f (x) + K 1 (x, t)N1 (t, u(t))dt + K 2 (x, t)N2 (t, u(t))dt + K 3 (x, t)N3 (t, u(t))d W (t), (1) 0

0

0

where f (x), K 1 (x, t), K 2 (x, t), K 3 (x, t), N1 , N2 , N3 are the stochastic processes defined on the same filtered probability space (Ω , F, Ft , P) with natural filtration {Ft } and u(t) is an ∫ x unknown stochastic process which ( ) should be computed. W (t), t ≥ 0 is a n-Brownian motion process and K 3 (x, t)N3 (t, u(t))d W (t) is the Itˆo 0

integral [14]. Numerical solution of deterministic integral equations of type (1) has been the subject of several studies, such as collocation method based on the radial basis functions [27], Legendre collocation method [26] and Chebyshev polynomials approach [3]. The outline of this paper is as follows. First, some concepts of stochastic calculus are given in Section 2. In Section 3, we derive operational matrix of integration P, the operational matrix of the product for Bernoulli polynomials and we introduce function approximation. Section 4, is devoted to the numerical solution of stochastic Volterra–Fredholm–Hammerstein integral equation. The convergence analysis and accuracy estimation are guaranteed in Section 5. In Section 6, we report our results on numerical examples to demonstrate the accuracy of the proposed method. Finally, we give a brief conclusion in Section 7. 2. Background 2.1. Stochastic concepts In this section, we review some basis properties of the stochastic calculus that are essential for this work. Let (Ω , F, P) denote a probability measure space, W (t, ω) be the random process and L 2 (Ω , F, P) be the space of all random variables u(t, ω). We will assume that for each t ∈ R+ , a minimal σ - algebra Ft , Ft ⊂ F, such that W (t, ω) is a measurable with respect of Ft . In addition, we will assume that the σ -algebra Ft is an increasing { }1/2 family such that u(t, ω) is Ft measurable and that E|u(t, ω)|2 < ∞, for each t ∈ R+ . Also we denote (∫ )1/2 { } 2 1/2 2 E|u(t, ω)| = ∥u(t, ω)∥ L 2 (Ω,F ,P) = |u(t, ω)| d P(ω) . (2) Ω

Definition 1 ([8]). We shall denote by Cc = Cc (R+ , L 2 (Ω , F, P)) the space of all continuous maps u(t, ω) from R+ into L 2 (Ω , F, P) such that for each t ≥ 0, u(t, ω) is Ft measurable. We define a topology in the space Cc by means of the following family of seminorms { }1/2 , n = 1, 2, . . . . (3) ∥u(t, ω)∥n = sup E|u(t, ω)|2 0≤t≤n

It is known that such topology is metrizable and that the metric space Cc is complete. { }1/2 Definition 2 ([8]). Let C ⊂ Cc be the space of maps u(t, ω) on R+ with E|u(t, ω)|2 < M, for some M > 0. The norm in the space C is defined by { }1/2 ∥u(t, ω)∥C = sup E|u(t, ω)|2 . (4) t≥0

240

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

Definition 3. {W (t), t ≥ 0} be the stochastic Brownian motion(SBM) with the corresponding properties as follows: 1. SBM has independent increments for 0 ≤ t0 ≤ t1 , · · · ≤ tn ≤ T . 2. W (t + h) − W (t) be normally distributed with mean 0 and variance h for all t ≥ 0, h > 0. 3. W (t) has continuous trajectory. Definition 4. Let M = M(S, T ) be the class of functions g(t, ω) : [0, ∞) −→ R such that: (1) The function g(t, ω) be B × F measurable, where B is the Borel σ -algebra of R+ . (2) The g(t,]ω) is Ft − adapted. [∫ function T 2 (3) E g (t, ω)dt < ∞. S

Theorem 1 (The Itô Isometry). Let g ∈ M(S, T ), then [∫ T ]2 [∫ T ] 2 g(t, ω)dW (t) = E g (t, ω)dt . E S

(5)

S

Proof of Theorem 1. See [14]. Lemma 1 (The Gronwall Inequality [14]). Let α, β ∈ [t0 , T ] −→ R be integrable with ∫ t 0 ≤ α(t) ≤ β(t) + L α(s)ds, t0

for t ∈ [t0 , T ] where L > 0. Then ∫ t α(t) ≤ β(t) + L e L(t−s) β(s)ds, t0

for t ∈ [t0 , T ]. 3. Bernoulli polynomials The Bernoulli polynomials Bn (x), which are usually defined by the exponential generating function [28] ∞

∑ te xt tn = B (x) n et − 1 n! n=0

(|t| < 2π ).

(6)

play an important role in different areas of mathematics, including number theory of finite differences, approximation theory, in the study of vertex algebras and other branches of mathematical analysis. They are constructed from the following relation ) n ( ∑ n+1 Bk (x) = (n + 1)x n , n = 0, 1, .... k k=0 and satisfy the following properties which will be of fundamental importance in this work. 1. Differentiation Bn′ (x) = n Bn−1 (x),

n ≥ 1,

(7)

2. Integral means conditions ∫ 1 Bn (x)d x = 0, n ≥ 1,

(8)

0

3. Differences Bn (x + 1) − Bn (x) = nx n−1 ,

n ≥ 1,

(9)

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

4. Series representation in terms of Bernoulli numbers n ( ) ∑ n Bn (x) = Bk (0)x n−k , n ≥ 1. k k=0

241

(10)

5. The polynomials Bn (x) satisfy the differential equation [25]. (1 ) Bn−1 (0) (n−1) B2 (0) ′′ Bn (0) (n) y (x) + y (x) + · · · y (x) + − x y ′ (x) + ny(x) = 0. n! (n − 1)! 2! 2 The Bernoulli polynomials form a complete basis over the interval [0, 1] [15].

(11)

3.1. Operational matrix of integration In this section, we derive the structure of deterministic operational matrix of Bernoulli polynomials. Let B(x) = [B0 (x), B1 (x), . . . , Bn (x)]T then B(x) = D −1 Tn (x),

(12)

where 0

0

0

...

0

1 2 2 0

1 2 2 1

0

...

0

1 3 3 0

1 3 3 1

1 3 3 2

()

...

0

.. .

.. .

.. .

.. .

.. .

...

1 n+1 n+1 n

(1 )

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ D=⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝

()

()

()

1 n+1 n+1 0

(

()

)

1 n+1 n+1 1

(

)

1 n+1 n+1 2

(

)

(

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ )

and Tn (x) = [1, x, x 2 , . . . , x n ]T . We can write ˆ n (x), B(x) = DT

(13)

where ⎛ (0) 0

B0 (0)

⎜ ⎜ (1) ⎜ ⎜ 1 B1 (0) ⎜ ⎜ (2) ⎜ D = ⎜ 2 B2 (0) ⎜ ⎜ ⎜ .. ⎜ . ⎜ ⎝ (n ) Bn (0) n

0 (1 ) 0

(2 ) 1

0

(2 ) 0

...

0

B0 (0)

...

0

.. .

.. .

.. . (

n n−1

)

Bn−1 (0)

.. . (

n n−2

0

0

B0 (0) B1 (0)

...

)

Bn−2 (0)

...

(n ) 0

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

B0 (0)

Since D is lower triangular matrix with nonzero diagonal elements, it is nonsingular and hence D −1 exists. By (12) and (13), it is easy to get Dˆ = D −1 . The dual matrix of B(x) is given by ∫ 1 Q= B(x)B T (x)d x = Dˆ H Dˆ T , (14) 0

where H is the Hilbert matrix of order (n + 1) [4]. To obtain the deterministic Bernoulli operational matrix of integration, we have ∫ x 1 (Bn+1 (x) − Bn+1 (0)), n ≥ 0 Bn (t)dt = n + 1 0

242

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

so x

∫

B(t)dt = P B(x) + 0

1 Bn+1 (x)en+1 , n+1

where en+1 = (0, 0, 0, . . . , 1) and P ⎛ −B1 (0) 1 0 ⎜ ⎜ −B2 (0) 1 ⎜ 0 2 2 ⎜ ⎜ ⎜ .. .. .. ⎜ P=⎜ . . . ⎜ ⎜ ⎜ −B (0) n ⎜ 0 0 ⎜ n ⎝ −Bn+1 (0) n+1

0

is the (n +1)×(n +1) Bernoulli deterministic operational matrix of integration, ⎞ ... 0 ⎟ ⎟ ... 0 ⎟ ⎟ ⎟ ⎟ .. .. ⎟ . . ⎟ ⎟ ⎟ ⎟ 1 ⎟ ... n ⎟ ⎠

0

...

0

The integration of the vector B(x) gives ∫ x B(t)dt ≃ P B(x),

(15)

0

and we approximate the integral of u(x) by ∫ x ∫ x U T B(t)dt = U T P B(x). u(t)dt =

(16)

0

0

For the evaluation of the product matrix of B(x) and B T (x) of Bernoulli polynomials basis (see [9]), we have π (x)C = B(x)B T (x)C ≃ Cˆ B(x),

(17)

where C = [c0 , c1 , . . . , cn ] is an arbitrary vector in R T

Cˆ = Dˆ C˜ T , C˜ = [ E˜ 0 , E˜ 1 . . . , E˜ n ]T ,

n+1

E˜ k = E k C, [e0k,i ,

, k = 0, 1, . . . , n,

e1k,i , . . . , enk,i ]T

and E k is an (n + 1) × (n + 1) matrix with ek,i = vector coefficients of x k Bi (x) ≃ B T (x)ek,i , and ∫ 1 ∫ ( 1 ) k x Bi (x)B(x)d x ≃ B(x)B T (x)d x ek,i = Qek,i , 0

as its columns, where ek,i are the Bernoulli

i, k = 0, 1, . . . , n

(18)

0

where ek,i are given by ∫ 1 ek,i ≃ Q −1 x k B(x)Bi (x)d x.

(19)

0

3.2. Function approximation Suppose that H = L 2 [0, 1] and {B0 (x), B1 (x), . . . , Bn (x)} ⊂ H be the set of Bernoulli polynomials and Y = span {B0 (x), B1 (x), . . . , Bn (x)}, let u ∈ H , since Y is a finite dimensional vector space, u has a unique best approximation that belongs to Y such as uˆ ∈ Y , that is ∀y ∈ Y,

∥u − u∥ ˆ ≤ ∥u − y∥,

(20)

and u ≃ uˆ =

n ∑

u k Bk (x) = B T (x)U,

(21)

k=0

where B(x) = [B0 (x), B1 (x), . . . , Bn (x)]T and U = [u 0 , u 1 , . . . , u n ]T . Since uˆ ∈ Y , there exist unique coefficients u 0 , u 1 , . . . , u n such that

243

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

Theorem 2. Assume that u ∈ H be an arbitrary function and also is approximated by the truncated Bernoulli n ∑ series u k Bk (x) then the coefficients u k for all k = 0, . . . , n can be calculated from the following relation k=0

( ) U T = u(x), B(x) Q −1 , (∫ 1 ∫ u(x)B0 (x)d x, where (u(x), B(x)) = 0

1

u(x)B1 (x)d x, . . .

1

∫

0

) u(x)Bn (x)d x .

0

Proof of Theorem 2. Since H = L [0, 1] is a Hilbert space with inner product (u, v) = 2

n ∑

u(x)v(x)d x, any 0

function u ∈ H can be expanded in Bernoulli basis as u ≃ uˆ =

1

∫

u k Bk (x) = B T (x)U = U T B(x),

k=0

we have u(x)Bi (x) =

n ∑

u k Bk (x)Bi (x),

i =0:n

k=0

so ∫

1

u(x)Bi (x)d x = 0

n ∑

∫

Bk (x)Bi (x)d x,

0

1

i = 0, . . . , n.

0

k=0

Consequently, we have (∫ 1 ∫ u(x)B0 (x)d x,

1

uk

u(x)B1 (x)d x, . . .

0

∫

)T

1

u(x)Bn (x)d x

= QU,

0

or (∫

1

1

∫ u(x)B0 (x)d x,

u(x)B1 (x)d x, . . .

Then, we get, (∫ T U =

1

∫ u(x)B0 (x)d x,

0

1

) u(x)Bn (x)d x

= U T Q T = U T Q.

0

0

0

∫

1

u(x)B1 (x)d x, . . .

0

∫

1

)

u(x)Bn (x)d x Q −1 = (u(x), B(x))Q −1 .

0

Theorem 3. Assume that Ψ (x, t) ∈ L 2 ([0, 1] × [0, 1]) be an arbitrary function and also is approximated by the two variables truncated Bernoulli series n ∑ n ∑ Pn [Ψ ](x, t) = Ψmk Bm (x)Bk (t) = B T (x)Ψ B(t), m=0 k=0

then the matrix Ψ can be calculated from the following relation [∫ 1 [∫ 1 ] ] Ψ = Q −1 B(x)Ψ (x, t)d x B T (t)dt Q −1 . 0

0

Proof of Theorem 3. Let Ψ (x, t) = B T (x)Ψ B(t) then we have B(x)Ψ (x, t)B T (t) = B(x)B T (x)Ψ B(t)B T (t). On the other hand, we get [∫ 1 ] [∫ B(x)Ψ (x, t)d x B T (t) = 0

1 0

] B(x)B T (x)d x Ψ B(t)B T (t) = QΨ B(t)B T (t),

244

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

Hence, ∫ 1 [∫

1

]

∫ 1 [∫

T

B(x)Ψ (x, t)d x B (t)dt = 0

0

0

∫

1

]

T

B(x)B (x)d x Ψ B(t)B T (t)dt

0

1

B(t)B T (t)dt = QΨ Q,

= QΨ 0

which implies Ψ=Q

−1

[∫ 1 [∫ 0

1

]

]

B(x)Ψ (x, t)d x B (t)dt Q −1 , T

0

This completes the proof. 4. Numerical solution of nonlinear stochastic Volterra–Fredholm–Hammerstein integral equation Let us give stochastic nonlinear Volterra–Fredholm–Hammerstein integral equation ∫ x ∫ 1 ∫ x u(x) = f (x) + K 1 (x, t)N1 (t, u(t))dt + K 2 (x, t)N2 (t, u(t))dt + K 3 (x, t)N3 (t, u(t))d W (t) 0

0

(22)

0

where f (x), K 1 (x, t), K 2 (x, t), K 3 (x, t), N1 , N2 and N3 are∫the stochastic processes defined on the same probability x K 3 (x, t)N3 (t, u(t))d W (t) is the Itˆo integral. space (Ω , F, P). W (t) is a Brownian motion process and 0

First, we define ⎧ ∫ x ∫ 1 ∫ x ⎪ ⎪ z (x) = N (x, f (x) + K (x, t)z (t)dt + K (x, t)z (t)dt + K 3 (x, t)z 3 (t)d W (t)) ⎪ 1 1 1 1 2 2 ⎪ ⎪ ⎪ ∫0 x ∫0 1 ∫0 x ⎨ K 3 (x, t)z 3 (t)d W (t)) K 2 (x, t)z 2 (t)dt + K 1 (x, t)z 1 (t)dt + z 2 (x) = N2 (x, f (x) + ⎪ ⎪ 0 0 0 ∫ ∫ ∫ ⎪ x 1 x ⎪ ⎪ ⎪ ⎩ z 3 (x) = N3 (x, f (x) + K 1 (x, t)z 1 (t)dt + K 2 (x, t)z 2 (t)dt + K 3 (x, t)z 3 (t)d W (t)). 0

0

(23)

0

We approximate the functions f (x), K 1 (x, t), K 2 (x, t), K 3 (x, t), z 1 (x), z 2 (x) and z 3 (x) in terms of Bernoulli polynomials as follows f (x) ≃ B T (x)F,

z 1 (x) ≃ B T (x)Z 1 ,

K 1 (x, t) ≃ B T (x)K 1 B(t),

z 2 (x) ≃ B T (x)Z 2 ,

K 2 (x, t) ≃ B T (x)K 2 B(t),

z 3 (x) ≃ B T (x)Z 3 ,

K 3 (x, t) ≃ B T (x)K 3 B(t),

where the vectors F, Z 1 , Z 2 , Z 3 and matrix K 1 , K 2 and K 3 are stochastic Bernoulli polynomial coefficients of f (x), z 1 (x), z 2 (x), z 3 (x) and K 1 (x, t), K 2 (x, t), K 3 (x, t) respectively. We have ∫ x ∫ x K 1 (x, t)z 1 (t)dt ≃ B T (x)K 1 B(t)B T (t)Z 1 dt (24) 0 0 ∫ x ∫ x = B T (x)K 1 π (t)Z 1 dt = B T (x)K 1 Zˆ 1 B(t)dt = B T (x)K 1 Zˆ 1 P B(x), 0

0

and ∫

1

∫

1

B T (x)K 2 B(t)B T (t)Z 2 dt ∫ 1 = B T (x)K 2 B(t)B T (t)Z 2 dt = B T (x)K 2 Q Z 2 ,

K 2 (x, t)z 2 (t)dt ≃ 0

(25)

0

0

where Q is defined in (14). For the stochastic term, we have ∫ x ∫ x K 3 (x, t)z 3 (t)d W (t) ≃ B T (x)K 3 B(t)B T (t)Z 3 d W (t) 0 0 ∫ x ∫ = B T (x)K 3 π (t)Z 3 d W (t) = B T (x)K 3 Zˆ 3 0

0

x

B(t)d W (t)

(26)

245

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

= B T (x)K 3 Zˆ 3

∫

x

∫

x

ˆ n (t)d W (t) = B T (x)K 3 Zˆ 3 Dˆ DT Tn (t)d W (t) 0 0 [∫ x ]T ∫ x ∫ x T n ˆ ˆ = B (x)K 3 Z 3 D d W (t), td W (t), . . . , t d W (t) 0

We can write Dˆ

x

[∫

td W (t), . . . ,

d W (t), 0

⎛ ∫

x

∫

x

d W (t)

⎜ ∫0 ⎜ x ⎜ td W (t) ⎜ ⎜ ∫0 ⎜ x ⎜ t 2 d W (t) ⎜ ⎜ 0 ⎜ ˙. ⎜ ⎜ ˙. ⎜ ∫ x ⎝ n t d W (t)

x

∫

0

t n d W (t)

0

0

]T in the form

0

⎞

⎛ ∫ x 0 ⎟ ⎜ ⎟ ⎜ ⎟ W (t)dt ⎜ ⎟ ∫ 0x ⎜ ⎟ ⎜ ⎟ t W (t)dt ⎜ ⎟ ⎟ = W (x)Tn (x) − ⎜ 0 ⎜ ⎟ ˙. ⎜ ⎟ ⎜ ⎟ ˙. ⎜ ∫ ⎟ x ⎟ ⎝ ⎠ n−1 t W (t)dt n

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ = An (t) = (ai )i=0,...,n ⎟ ⎟ ⎟ ⎟ ⎠

0

0 i

∫

where a0 = W (x) and ai = x W (x) − i

x

t

i−1

i = 1, . . . , n. For the integral

W (t)dt,

0

use composite trapezium rule as follows ) ( ∫ x x x x x i−1 W (x)dt ≃ 2( )i−1 W ( ) + x i−1 W (x) , i = 1, . . . , n. 4 2 2 0 so ( ) [ ] x x x i−1 x i i i i−1 ai = x W (x) − i 2( ) W ( ) + x W (x) = (1 − )W (x) − i W ( ) x i , 4 2 2 4 2 2

x

∫

t i−1 W (t)dt, we can 0

i = 1, . . . , n.

Also we approximate W (x) and W ( x2 ) for 0 ≤ x ≤ 1 by W (0.5) and W (0.25), then we obtain ⎛ ⎛ ⎜ ⎜ Dˆ An (x) = Dˆ ⎜ ⎜ ⎝

W (0.5) 0 ˙. ˙. ˙.

0 3 1 W (0.5) − W (0.25) 4 2 ˙. ˙. ˙.

0 0 ˙. ˙. ˙.

... ... ... ... ...

Hence Dˆ An (x) = Dˆ Ds Tn (x) = Dˆ Ds Dˆ −1 B(x) = Ps B(x), where ⎛ W (0.5) 0 0 ... 3 1 ⎜ 0 W (0.5) − W (0.25) 0 ... ⎜ 4 2 Ds = ⎜ ˙ . ˙ . ˙ . ... ⎜ ⎝ ˙. ˙. ˙. ... ˙. ˙. ˙. ...

0 0 ˙. ˙. (1 − n4 )W (0.5) −

0 0 ˙. ˙. (1 − n4 )W (0.5) −

⎞⎜ ⎜ ⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎠⎜ ⎜ n ⎜ W (0.25) 2n ⎝

⎞

n 2n

⎟ ⎟ ⎟, ⎟ ⎠ W (0.25)

and Ps = Dˆ Ds Dˆ −1 is (n + 1) × (n + 1) Bernoulli stochastic operational matrix. Then ∫ x B(t)d W (t) ≃ Ps B(x). 0

Consequently Eqs. (24), (25) and (26) yield ⎧ T ⎨ B (x)Z 1 = N1 (x, B T (x)F + B T (x)K 1 Zˆ 1T P B(x) + B T (x)K 2 Q Z 2 + B T (x)K 3 Zˆ 3T Ps B(x)), B T (x)Z 2 = N2 (x, B T (x)F + B T (x)K 1 Zˆ 1T P B(x) + B T (x)K 2 Q Z 2 + B T (x)K 3 Zˆ 3T Ps B(x)), ⎩ T B (x)Z 3 = N3 (x, B T (x)F + B T (x)K 1 Zˆ 1T P B(x) + B T (x)K 2 Q Z 2 + B T (x)K 3 Zˆ 3T Ps B(x)).

1 x x2 . . . . xn

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

246

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

Now, we collocate the previous system at n + 1 Newton–Cotes nodes as (2i − 1) , i = 1, . . . , n + 1, xi = 2(n + 1) we obtain ⎧ T ⎨ B (xi )Z 1 = N1 (xi , B T (xi )F + B T (xi )K 1 Zˆ 1T P B(xi ) + B T (xi )K 2 Q Z 2 + B T (xi )K 3 Zˆ 3T Ps B(xi )), B T (x )Z = N2 (xi , B T (xi )F + B T (xi )K 1 Zˆ 1T P B(xi ) + B T (xi )K 2 Q Z 2 + B T (xi )K 3 Zˆ 3T Ps B(xi )), ⎩ T i 2 B (xi )Z 3 = N3 (xi , B T (xi )F + B T (xi )K 1 Zˆ 1T P B(xi ) + B T (xi )K 2 Q Z 2 + B T (xi )K 3 Zˆ 3T Ps B(xi )), which corresponds to a nonlinear system with unknown Z 1 , Z 2 and Z 3 . This system can be solved by using an appropriate numerical method such as Newton’s method. After solving this nonlinear system of algebraic equations and finding unknown vectors Z 1 , Z 2 and Z 3 , the approximate solution of Eq. (22) can be obtained as u(x) ≃ B T (x)U = f (x) + B T (x)K 1 Zˆ 1T P B(x) + B T (x)K 2 Q Z 2 + B T (x)K 3 Zˆ 3T Ps B(x).

(27)

5. Error bounds In this section, we show that with an increase in the number of the Bernoulli basis, error of the approximation approaches 0. Theorem 4 ([2]). Suppose that u(x) be an enough smooth function in the interval [0, 1] and Pn [u](x) is the approximate polynomial of u(x) in terms of Bernoulli polynomials and Rn (x) is the remainder term. Then the associated formulae are stated as follows u(x) = Pn [u](x) + Rn (x), ∫ Pn [u](x) =

1

u(x)d x + 0

x ∈ [0, 1],

n ∑ B j (x) ( j−1) (u (1) − u ( j−1) (0)), j! j=1

∫ 1 1 ⋆ Rn (x) = − B (x − t)u (n) (t)dt n! 0 n where Bn⋆ (x) = Bn (x − [x]) and [x] denote the largest integer not greater than x. Theorem 5 ([30]). Suppose u ∈ C ∞ [0, 1] and Pn [u](x) is the approximate polynomial using Bernoulli polynomials. Then the error bound would be obtained as follows 1 E(u) = ∥u(x) − Pn [u](x)∥∞ ≤ Bn Un , n! where Bn and Un denote the maximum value of Bn (x) and u (n) (x) in the interval [0, 1] respectively. Theorem 6 ([10]). Suppose g(x, t) be a smooth enough function and P[g](x, t) is the approximate polynomial using Bernoulli polynomials. Then the error bound would be obtained as follows 1 E(g) = ∥g(x, t) − P[g](x, t)∥∞ ≤ B 2 G n,n , (n!)2 n ∂ 2n g(x, t) in the interval [0, 1] × [0, 1] respectively. ∂ x n ∂t n In the following theorem, we give an upper bound for the error and we investigate the convergence analysis of the proposed method. Let u(x) and u n (x) be the exact and approximate solutions of the considered equation.

where Bn and G n,n denote the maximum values of Bn (x) and

Theorem 7. Suppose u(x) and u n (x) be the exact and approximate solutions of Eq. (22), respectively, furthermore, we suppose that |N1 (x, u 1 (x)) − N1 (x, u 2 (x))| + |N2 (x, u 1 (x)) − N2 (x, u 2 (x))| + |N3 (x, u 1 (x)) − N3 (x, u 2 (x))| ≤ L|u 1 − u 2 |. (28)

|N1 (x, u(x))| + |N2 (x, u(x))| + |N3 (x, u(x))| ≤ L(1 + |u|).

(29)

247

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

Let 1. ∥u(x)∥ ≤ M, x ∈ [0, 1]. 2. ∥K i (x, ( t)∥ ≤ Mi , i )= 1, 2, 3. 3. 16L 2 M22 + 2λ22 (n) < 1.

x, t ∈ [0, 1].

∂ 2n K 2 (x, t) 1 2 2 , with Bn and Kˆ n,n Bn2 Kˆ n,n denote the maximum values of Bn (x) and in the interval 2 (n!) ∂ x n ∂t n [0, 1] respectively Then we have where λ2 (n) =

∥u(x) − u n (x)∥ −→ 0, where ( ) ∥u(x)∥2 = E |u(x)|2 . Proof of Theorem 4. Let zˆ i (x), i = 1, 2, 3 be the approximated form of z i (x) by Bernoulli polynomials i.e., zˆ i (x) = Nˆ i (x, u n (x)), i = 1, 2, 3, (30) and z in (x) = Ni (x, u n (x)), i = 1, 2, 3. (31) 1 1 1 i Let also η(n) = B 2 Kˆ i , γi (n) = Bn fˆn , λi (n) = Bn Zˆ ni , where Bn , fˆn , Zˆ n and Kˆ n,n denote the n! (n!)2 n n,n n! ∂ 2n K i (x, t) maximum values of Bn (x), f (n) (x), z i(n) (x) and in the interval [0, 1] respectively, for i = 1, 2, 3. ∂ x n ∂t n By Lipschitz condition (28) and Theorem 5, we have |ˆz i (t) − z i (t)|2 ≤ 2|ˆz i (t) − z in (t)|2 + 2|z in (t) − z i (t)|2 ≤ 2γi2 (n) + 2L 2 |u(t) − u n (t)|2 .

(32)

Let (u(x) − u n (x)) be the error function, then using Eq. (22) we have ∫ x (K 1 (x, t)z 1 (t) − Kˆ t (x, t)ˆz 1 (t))dt u(x) − u n (x) = f (x) − f n (x) + 0 ∫ 1 (K 2 (x, t)z 2 (t) − Kˆ 2 (x, t)ˆz 2 (t))dt + ∫0 x + (K 3 (x, t)z 3 (t) − Kˆ 3 (x, t)ˆz 3 (t))d W (t). 0

Via the inequality (a + b + c + d)2 ≤ 4a 2 + 4b2 + 4c2 + 4d 2 , we get ( ∫ x 2 ∥u(x) − u n (x)∥2 = E|u(x) − u n (x)|2 ≤ 4E | f (x) − f n (x)|2 + | (K 1 (x, t)z 1 (t) − Kˆ 1 (x, t)ˆz 1 (t))dt| 0

∫

1

2

(K 2 (x, t)z 2 (t) − Kˆ 2 (x, t)ˆz 2 (t))dt|

+| 0

∫ +|

x

(K 3 (x, t)z 3 (t) − Kˆ 3 (x, t)ˆz 3 (t))d W (t)|

2)

.

0

By using Schwartz inequality, Itˆo isometry and linearity of Itˆo integrals in their integrands, we get ( ∫ x 2 2 2 2 ∥u(x) − u n (x)∥ = E|u(x) − u n (x)| ≤ 4 E| f (x) − f n (x)| + E|(K 1 (x, t)z 1 (t) − Kˆ 1 (x, t)ˆz 1 (t))| dt 0 ∫ 1 2 + E|(K 2 (x, t)z 2 (t) − Kˆ 2 (x, t)ˆz 2 (t))| dt ) ∫0 x 2 + E|(K 3 (x, t)z 3 (t) − Kˆ 3 (x, t)ˆz 3 (t))| dt , 0

248

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

It follows from the inequality (a + b)2 ≤ 2a 2 + 2b2 that ∫ x ∫ x ( ) ( ) 2 2 ˆ E|(K 1 (x, t)z 1 (t) − K 1 (x, t)ˆz 1 (t))| dt = E|K 1 (x, t) z 1 (t) − zˆ 1 (t) − Kˆ 1 (x, t) − K 1 (x, t) zˆ 1 (t)| dt 0 0 ∫ x ( )2 ≤2 E|K 1 (x, t) z 1 (t) − zˆ 1 (t) | dt ∫0 x ( ) 2 E| Kˆ 1 (x, t) − K 1 (x, t) zˆ 1 (t)| dt +2 0

Thus from (28), (29) and (32) it follows that [∫ x ] [∫ x ] [∫ x ] |ˆz i (t)|2 dt ≤ 2E |ˆz i (t) − z i (t)|2 dt + 2E |z i (t)|2 dt E 0 0 ) ) ( ( 0 ∫ x 2 2 E|u(t) − u n (t)|2 + 2 L 2 (1 + M)2 , i = 1, 3. ≤ 2 2γi (n) + 2L

(33)

0

On the other hand, by using Theorems 5, 6, assumptions (i), (ii), inequality (33) and mean value theorem and noting that t, x ∈ [0, 1] we obtain ] [ ∫ x ∫ x 2 E|u(t) − u n (t)|2 dt (34) E|(K 1 (x, t)z 1 (t) − Kˆ 1 (x, t)ˆz 1 (t))| dt ≤ 2M12 2γ12 (n) + 2L 2 0 0 [ ] ∫ x ( )2 + 2λ21 (n) 4γ12 (n) + 4L 2 E|u(t) − u n (t)|2 dt + 4L 2 1 + M , 0

and ∫ 0

1

[ ] 2 2 2 2 2 ˆ E|(K 2 (x, t)z 1 (t) − K 2 (x, t)ˆz 2 (t))| dt ≤ 2M2 2γ2 (n) + 2L E|u(x) − u n (x)| [ ] ( )2 2 2 2 2 2 + 2λ2 (n) 4γ2 (n) + 4L E|u(x) − u n (x)| + 4L 1 + M .

(35)

Also it can be shown that [ ] ∫ x ∫ x 2 E|(K 3 (x, t)z 1 (t) − Kˆ 3 (x, t)ˆz 3 (t))| dt ≤ 2M32 2γ32 (n) + 2L 2 E|u(t) − u n (t)|2 dt (36) 0 0 [ ] ∫ x ( )2 2 2 2 2 2 + 2λ3 (n) 4γ3 (n) + 4L E|u(t) − u n (t)| dt + 4L 1 + M , 0

Combining (34), (35) and (36), we obtain ∫ x ∥u(x) − u n (x)∥2 ≤ β(x) + L ′ ∥u(t) − u n (t)∥2 dt, 0

where ( ) 16 M12 γ12 (n) + M22 γ22 (n) + M32 γ32 (n) + 4η(n)2 ( ) β(x) = ) ( 1 − 16L 2 M22 + 2λ22 (n) ( ) ( 2 ) ( 2 ) ( 2 ) 2 2 2 2 2 2 2 2 2 32 λ1 (n) γ1 (n) + L (1 + M) + λ2 (n) γ2 (n) + L (1 + M) + λ3 (n) γ3 (n) + L (1 + M) ( ) + ( 2 ) 2 2 1 − 16L M2 + 2λ2 (n) 16L L′ =

2

(

M12

+

M32 (

+

2λ21 (n)

+

2λ23 (n)

1 − 16L 2 M22 + 2λ22 (n)

)

) .

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

249

Table 1 Computed errors for Example 1 for different values of n. x

n=1

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

6.4285 2.3112 9.0049 3.2928 4.9436 5.9236 6.2965 6.1202 5.4469 4.3238

E-2 E-2 E-3 E-2 E-2 E-2 E-2 E-2 E-2 E-2

n=3

n=5

n=7

n = 13

1.6312 E-1 5.6182 E-2 8.6497E−3 3.9639 E-2 4.4967 E-2 3.2738E−2 1.0991 E-2 1.2294 E-2 2.9196 E-2 3.1841 E-2

1.1390 E-3 2.0833 E-3 1.3769 E-3 2.0626 E-4 7.3314 E-4 1.1363 E-3 9.8885 E-4 4.6763 E-4 1.1363 E-3 6.2285E−4

9.4532 E-5 1.1466 E-4 8.6272E−5 3.1485 E-5 2.2314 E-5 5.5158 E-5 5.9614 E-5 4.0261 E-5 9.4649 E-6 1.862 E-5

6.0902 E-8 2.6170 E-7 4.3312 E-7 4.1989 E-7 2.6011 E-7 3.6961 E-8 1.5895 E-7 2.6304E−7 2.5498 E-7 1.5792E−7

Table 2 Computed errors for Example 1 for different values of n at Newton–Cotes nodes. x

n=3

x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14

3.6351 4.5596 5.0038 3.2945

n=5 E-2 E-2 E-3 E-2

2.0807 7.9580 8.4103 1.0458 1.5019 6.7107

n=7 E-3 E-4 E-4 E-3 E-4 E-4

1.1340 9.1839 2.4192 3.7775 6.1205 4.3603 5.5388 2.6163

n=9 E-4 E-5 E-5 E-5 E-5 E-5 E-6 E-5

4.1887 4.7083 3.3587 1.0815 1.0759 2.3566 2.4906 1.6797 4.2265 7.2598

n = 13 E-6 E-6 E-6 E-6 E-6 E-6 E-6 E-6 E-7 E-7

6.6342 E-8 2.7981 E-7 4.1198 E-7 4.4866 E-7 3.9559 E-7 2.7474 E-7 1.1797 E-7 4.0308 E-8 1.6994 E-7 2.5021 E-7 2.7247 E-7 2.4022 E-7 1.668 E-7 7.1584 E-8

The Gronwall inequality for this inequality results ( ) ∫ x ′ ∥u(x) − u n (x)∥2 ≤ β(x) 1 + L ′ e L (x−t) dt , 0

by increasing n, it implies ∥u(x) − u n (x)∥2 −→ 0 as n −→ ∞, and so u n (x) converges to u(x) in L 2 . 6. Numerical examples To illustrate the method stated in section (4), we consider the following examples. In order to analyze the error method we introduce the absolute error which is defined as en (x) = |u(x) − u n (x)|, where u(x) and u n (x) denote exact and approximate solutions of order n respectively, Also, let xi , i = 1, . . . , n be the Newton–Cotes nodes. The programs have been provided by Matlab. Example 1. Let us give the following first kind Volterra integral equation ∫ x e(x+t) u(t)dt = xe x , x ∈ [0, 1], 0

the exact solution is u(x) = e−x . The numerical results for different values of n are given in Tables 1–3.

250

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254 Table 3 Comparison with the results obtained here with the results obtained in [7] for Example 1. x

Exact solution

BPFs method [7]

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.000000000000000 0.9048374180359600 0.8187307530779820 0.7408182206817179 0.6703200460356394 0.6065306597126334 0.5488116360940264 0.4965853037914095 0.4493289641172216 0.4065696597405991

Present method

m = 32

m = 64

N =9

0.989584 0.891689 0.820394 0.739236 0.680130 0.600213 0.540837 0.497594 0.448370 0.412520

0.994792 0.905768 0.824711 0.735426 0.669613 0.603372 0.549376 0.500213 0.446058 0.406141

1.0000000000000 0.904842141415613 0.818734969497211 0.740820488279243 0.670319976336825 0.60652880480573 0.548809072216355 0.496583127497036 0.449327893353760 0.406569855489044

Table 4 Computed errors for Example 2 for λ = −100, µ = 10, X 0 = 0.5 and M = 500 for different values of n at Newton–Cotes nodes. t

n=2

n=5

x1 x2 x3 x4 x5 x6 x7 x8 x9 x10

8.4296 E-2 2.4853 E-1 6.6464 E-2

1.6188 3.3984 1.1702 1.7605 3.9589 2.9200

n=9 E-2 E-2 E-2 E-2 E-2 E-2

1.1456 2.6473 3.5329 2.9812 1.1877 1.4339 3.9653 6.0913 7.9642 1.0846

E-2 E-2 E-2 E-2 E-2 E-2 E-2 E-2 E-2 E-1

Table 5 Computed errors for Example 2 for λ = −100, µ = 10, X 0 = 0.001 and M = 500 for different values of n. t

n=5

n=7

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.00 E-3 4.1575 E-5 6.7593 E-5 6.652 E-5 4.039 E-5 1.3454 E-5 4.5812 E-5 7.779 E-5 8.693 E-5 6.9115 E-5

1.0495 4.5245 9.4555 1.0767 7.5726 3.5219 6.8055 1.3656 2.5213 4.6582

n=9 E-3 E-5 E-5 E-4 E-5 E-5 E-5 E-4 E-4 E-4

1.0126 E-3 4.068 E-5 7.2105 E-5 7.4680 E-5 4.2611 E-5 3.132 E-4 1.0555 E-4 2.1409 E-4 3.6629 E-4 6.1828 E-4

Example 2 (The Basic Black–Scholes Model). Let us give the following linear stochastic equation d X (t) = λX (t)dt + µX (t)d W (t), X (0) = X 0 , t ∈ [0, 1],

(37)

where the exact solution is given by X (t) = ex p((λ − 12 µ2 )t + µW (t)). We take M simulations. The results obtained for this example are given in Tables 4–7. Example 3. problems).

Consider the nonlinear stochastic Itˆo–Volterra integral equation as follows (population growth t

∫ X (t) = 0.5 +

∫ X (s)(1 − X (s))ds +

0

t

X (s)d W (s). 0

(38)

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

251

Table 6 Computed errors for Example 2 for λ = −100, µ = 0 and n = 10 (left), λ = −100, µ = 10, n = 10 and M = 100 (right). t

X 0 = 0.1

X 0 = 0.001

X 0 = 0.1

X 0 = 0.001

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.0000000 8.4548 E-2 3.6769 E-2 2.4986 E-2 7.7114 E-2 9.9807 E-2 8.4470 E-2 3.6911 E-2 2.4815 E-2 7.7150 E-2

0.0000000 8.4548 E-4 3.6769 E-4 2.496 E-4 7.7114 E-4 9.9807 E-4 8.447 E-4 3.6911 E-4 2.4815 E-4 7.715 0E-4

1.041 E-1 3.7489 E-2 1.6409 E-2 1.1905 E-2 3.5135 E-2 4.5271 E-2 3.7757 E-2 1.5193 E-2 1.4992 E-2 4.114 E-2

1.041 E-3 3.7489E−4 1.6409 E-4 1.1905 E-4 3.5135 E-4 4.5271 E-4 3.7757 E-4 1.5193 E-4 1.4992 E-4 4.114 E-4

Table 7 Computed errors for Example 2 for λ = −5, µ = 1, n = 10 and M = 100. t

X 0 = 0.01

X 0 = 0.001

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

8.2717 E-2 2.8235 E-2 7.2744 E-3 1.4271 E-2 1.5735 E-2 1.700 E-2 1.6031 E-2 1.3352 E-2 9.4216 E-3 6.5464 E-3

8.2717 E-3 2.8235 E-3 7.2744 E-4 1.4271 E-3 1.5735 E-3 1.700 E-3 1.6031 E-3 1.3352 E-3 9.4216 E-4 6.5464 E-4

Table 8 Computed errors for Example 3 for different values of n and for M = 10 simulations. t

n=3

n=5

n=7

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

7.9519 E-2 2.1665 E-2 1.7132E−1 2.7696 E-1 3.1796 E-1 3.2587E−1 3.0767 E-1 1.0084 E-1 1.7871 E-1 1.6489 E-1

3.5834 E-1 3.5197 E-1 2.8920 E-1 1.2466 E-1 1.2362 E-1 1.8244 E-1 1.5324E−2 5.7701 E-3 1.0707 E-2 4.1180E−2

1.2203 E-1 8.4684 E-2 1.0078E−1 2.4080E−1 1.1590E−1 2.7728E−1 5.2755 E-2 1.2892 E-2 3.4658 E-1 4.4299E−1

with the exact solution X (t) =

e(0.5t+W (t)) , ∫ t (0.5t+W (t)) 2+ e ds 0

where X (t) is an unknown stochastic process defined on the probability space (Ω , F, P) and W (t) is a Brownian motion process. The numerical results are shown in Table 8. Example 4.

Consider the following linear stochastic Volterra integral equation. ∫ t ∫ t X (t) = 1 + s 2 X (s)ds + s X (s)d W (s) s, t ∈ [0, 0.5] 0

0

(39)

252

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254 Table 9 Computed errors for Example 4 for n = 5 and n = 7. Exact solution 1.000000 0.999695 0.998964 0.998601 1.008343 1.003699 1.019657 9.33198 1.048115 1.091446 1.074498

t 0 0.01 0.02 0.03 0.04 0.05 0.1 0.2 0.3 0.4 0.5

n=5

n=7

present method 1.010910 1.011107 1.011293 1.011469 1.011638 1.011801 1.012590 1.014707 1.018764 1.025983 1.037320

present method 1.002284 1.006794 1.011327 1.015885 1.020470 1.025083 1.011611 1.012479 1.013565 1.022765 1.033300

Table 10 The absolute errors of the approximate solution for Example 5 for n = 3. t

Absolute errors

t

Absolute errors

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

0.1032 0.0995 0.2754 0.1450 0.1008 0.0388 0.3199 0.1779 0.2292 0.4424

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1032 E-3 0.0005 0.0011 0.0007 0.0001 0.0003 0.0012 0.0008 0.0009 0.0016

E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3

) ∫ t t3 with the exact solution X (t) = exp + sd W (s) for t ∈ [0, 0.5], where X (t) is an unknown stochastic 6 0 process defined in probability space (Ω , F, P). The numerical results are shown in Table 9. (

Example 5.

Consider the following problem [14] ∫ t ∫ t 2 3 sin2 (X (s))d W (s). X (t) = X 0 + a cos(X (s)) sin (X (s))ds − a 0

(40)

0

The exact solution of this problem is X (t) = ar ccot(aW (s)+cot(X 0 )). This problem is now solved by the proposed method for a = 1/8 and X 0 = π/32 for n = 3. The absolute errors of the approximate solution at some different points t ∈ [0, 1] are given in Table 10. Example 6.

Consider the following nonlinear stochastic integral equation [14]

d X t = a 2 X t (1 + X t2 )dt + a(1 + X t2 )d Wt ,

(41)

The exact solution of this problem is X (t) = tan(aWt + arctan(X 0 )). We have solved this problem by the proposed method for a = 0.001, X 0 = 1 and n = 5. The absolute errors of the approximate solution at some different points t ∈ [0, 1] are shown in Table 11 and the results for a = 1/104 , X 0 = 0.01 are given in Table 12. 7. Conclusion and comments This paper suggests a numerical method to solve stochastic Volterra–Fredholm–Hammerstein integral equations (SVFHIE) by using Bernoulli polynomials and their operational matrix that was introduced in [9]. The stochastic operational matrix of Itˆo integration of Bernoulli polynomials was derived and applied for solving SVFHIE. The main advantage of this method is its efficiency and simple applicability, and based upon the numerical results,

253

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254 Table 11 The absolute errors of the approximate solution for Example 6 for n = 5, a = 0.001 and X 0 = 1. t

Absolute errors

t

Absolute errors

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

0.1702 0.1670 0.4583 0.2428 0.1651 0.0674 0.5328 0.2985 0.3835 0.7354

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.1702 E-3 0.0009 0.0018 0.0011 0.0002 0.0006 0.0021 0.0013 0.0016 0.0027

E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3

Table 12 The absolute errors of the approximate solution for Example 6 for n = 5, a = 1/104 and X 0 = 0.01. t

Absolute errors

t

Absolute errors

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

0.0852 0.0834 0.2292 0.1214 0.0827 0.0336 0.2664 0.1492 0.1917 0.3678

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.0852 0.0448 0.0909 0.0568 0.0077 0.0290 0.1027 0.0656 0.0790 0.1347

E-4 E-4 E-4 E-4 E-4 E-4 E-4 E-4 E-4 E-4

E-4 E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3 E-3

when the solution is sufficiently smooth, a small number of basis functions is enough to obtain a high accuracy solution. At the end, we note that the method can be easily extended and applied to stochastic multi-dimensional integral equations. We also believe to extend this method to linear and nonlinear Volterra–Fredholm–Hammerstein integro-differential systems which will be the subject of future research. Acknowledgments The author would like to thank the referees for their valuable comments and suggestions which have led to a much improved version of this paper. References [1] M.A. Abadon, Fredholm-Volterra integral equation with singular kernel, App. Math. Comput. 137 (2003) 231–243. [2] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, National Bureau of Standards, Wiley, New York, 1972. [3] A. Akyuz-Dascioglu, A Chebychev polynomial approach for linear Fredholm-Volterra integro-differential equations in most general form, Appl. Math. Comput. Sci. 3 (2011) 3826-389. [4] Alexandru Aleman, Montes-Rodriguez, The eigenfunctions of the Hilbert matrix, in: Constructive Approximation, vol. 36(3), Springer, 2012, pp. 353–374. [5] L. Arnold, Stochastic Differential Equations, Wiley, New York, 1974. [6] Mahnaz Asari, E. Hashemizadeh, M. Khodabin, M. Maleknejad, Numerical solution of nonlinear stochastic integral equations by stochastic operational matrix based on Bernstein polynomials, Bull. Math. Soc. Sci. Math. Roumanie. 57 (1) (2014) 3–12. [7] E. Babolian, Z. Masouri, Direct method to solve volterra integral equation of the first kind using operational matrix with block-pulse functions, J. Comput. Appl. Math. 220 (1–2) (2008) 51–57. [8] K. Balachandran, J.H. Kim, Existence of solutions on nonlinear stochastic volterra fredholm integral equations of mixed type, Int. J. Math. Math. Sci. (2010) http://dx.doi.org/10.1155/2010/603819. [9] S. Bazm, Bernoulli polynomials for the numerical solution of some classes of linear and nonlinear integral equations, J. Comput. Appl. Math. 275 (2015) 44–60. [10] A.H. Bhrawy, E. Tohidi, F. Soleymani, A new Bernoulli matrix method for solving high-order linear and nonlinear Fredholm integro-differential equations with piecewise intervals, Appl. Math. Comput. 219 (2) (2012) 482–497.

254

R. Zeghdane / Mathematics and Computers in Simulation 165 (2019) 238–254

[11] H. Brunner, J.P. Kauthen, The numerical solution of two-dimensional volterra integral equations by collocation and iterated collocation, IMA J. Numer. Anal. 9 (1989) 47–59. [12] H. Guoqiang, Asymptotic error expansion variation of a collocation method for Volterra-Hammerstein equations, J. Appl. Numer. Math. 13 (1993) 357–369. [13] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (3) (2001) 525–546. [14] P.E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 1992. [15] E. Kreyszig, Introductory Functional Analysis with Applications, John Wiley and Sons Press, New York, 1978. [16] S. Kumar, I.H. Shoan, A new collocation-type method for Hammerstein integral equations, J. Math. Comput. 48 (1987) 123–129. [17] D.H. Lehmer, A new approach to Bernoulli polynomials, Amer. Math. Monthly 95 (10) (1998) 905–911. [18] K. Maleknejad, E. Hashemizadeh, B. Basirat, Computational method based on Bernstein operational matrices for nonlinear Volterra–Fredholm– Hammerstein integral equations, Commun. Nonlinear Sci. Numer. Simul. 17 (1) (2012) 52–61. [19] K. Maleknejad, R. Mollapourasl, M. Alizadeh, Numerical solution of the Volterra type integral equation of the first kind with wavelet basis, Appl. Math. Comput. 194 (2) (2007) 400–405. [20] K. Maleknejad, B. Rahimi, Modification of block pulse functions and their application to solve numerically Volterra integral equation of the first kind, Commun. Nonlinear Sci. Numer. Simul. 16 (6) (2011) 2469–2477. [21] G.N. Milstein, Numerical Integration of Stochastic Differential Equations, Mathematics and Its Application, Kluwer, Dordrecht, 1995. [22] F. Mirzaee, Elham Hadadiyan, A collocation technique for solving nonlinear Stochastic Itô-Volterra integral equation, Appl. Math. Comput. 247 (2014) 1011–1020. [23] F. Mirzaee, Nasim Samadyar, Seyede Fatemeh Hoseini, Euler polynomial solutions of nonlinear stochastic Itô-Volterra integral equations, J. Comput. Appl. Math. 330 (2018) 574–585. [24] N.I. MuskHelishvili, Some Basic Problems of Mathematical Theory of Elasticity, Noord hoff, Groningen, 1953. [25] P. Natalini, A. Bernaridini, A generalization of the Bernoulli polynomials, J. Appl. Math. 3 (2003) 155–163. [26] S. Nemati, Numerical solution of Volterra-Fredholm integral equations using Legendre collocation method, J. Comput. Appl. Math. 278 (2015) 29–36. [27] K. Par, J.A. Rad, Numerical solution of nonlinear Volterra-Fredholm-Hammerstein integral equations via collocation method based on radial basis functions, Appl. Math. Comput. 218 (2012) 5292–5309. [28] H.M. Srivastava, J. Choi, Series Associated with the Zeta and Related Functions, Kluwer Academic, Dordreeht, 2001. [29] T.H. Tian, K. Burrage, Implicit Taylor methods for stiff stochastic differential equations, App. Numer. Math. 38 (2001) 167–185. [30] E. Tohidi, A.H. Bhrawy, K. Erfani, A collocation method based on Bernoulli operational matrix for numerical solution of generalized pantograph equation, Appl. Math. Model. 37 (6) (2013) 4283–4294. [31] F. Toutounian, E. Tohidi, A new Bernoulli matrix method for solving second order linear partial differential equations with the convergence analysis, Appl. Math. Comput. 223 (2013) 298–310.