www.i-csrs.org
Available free online at http://www.geman.in
Numerical Solution of Volterra-Fredholm Integro-Differential Equation by Block Pulse
Functions and Operational Matrices
∗Leyla Rahmani1, Bijan Rahimi1 and Mohammad Mordad2
1Department of Mathematics, Faculty of science, Islamic Azad University, Takestan branch, Takestan, Iran
E-mail: [email protected] E-mail: [email protected]
2Department of Mathematics, Faculty of science, Islamic Azad University, Islamshahr branch, Islamshahr, Tehran, Iran
E-mail: [email protected]
(Received:11-5-11/Accepted:20-6-11) Abstract
In this paper, Block Pulse Functions and their operational matrices are used to solve Volterra-Fredholm integro-differential equation (VFIDE). First the equation is integrated over interval [0, x] and then Block Pulse functions are used to obtain numerical solution. Some theorems will prove convergence of the method. Some numerical examples are included to illustrate accuracy of the method.
Keywords: Volterra-Fredholm integro-differential equation, Block Pulse Functions, Operational matrix.
1 Introduction
Integro-differential equations have been discussed in many applied fields, such as biological, physical and engineering problems. Therefore, their numerical treatment is desired. Many numerical methods used to solve such equations [1, 3, 4, 6, 7, 8, 10, 12, 14]. In this paper, we propose a new numerical method
to solve the linear VFIDE:
f0(x) =g(x) + Z x
0
k(x, t)f(t)dt+ Z 1
0
l(x, t)f(t)dt 0≤x, t <1, f(0) =a0,
(1)
wheref(x) andg(x) are inL2([0,1)) andk(x, t) andl(x, t) belong toL2([0,1)×
[0,1)). Moreover k(x, t),l(x, t) andg(x) are known and f(x) is unknown. We assume Eq. (1.1) has a unique solution.
The paper is organized as follows: In section 2, we describe Block Pulse functions, their properties and their operational matrices [5, 15]. In section 3, we propose a new numerical method for solving linear VFIDE. In section 4 we analysis the error. Finally in section 5 we apply the proposed method on some examples to show the accuracy and efficiency of the method.
2 Review of Block Pulse Functions
A set of Block Pulse functions (BPFs) are usually defined in [0,1) as:
φi(t) =
1, mi ≤t < (i+1)m , 0, otherwise,
(2)
wherei = 0,1,2, . . . , m−1 and m is an arbitrary positive integer. According to (2.1), the unit interval [0,1) is divided into m equidistant subintervals and theith Block Pulse function φi has only one rectangular pulse of unit hight in the subinterval [mi,i+1m ).
One of the important properties of BPFs is the disjointness of them, which can directly be obtained from the definition of BPFs. Indeed:
φi(t)φj(t) =
φi(t), i=j, 0, i6=j,
(3) wherei, j = 0,1, . . . , m−1.
The orthogonality of BPFs is derived immediately from:
Z 1 0
φi(t)φj(t)dt= h δij, (4) whereδij is the Kronecker delta and h= m1.
The set of BPFs is complete,i.e.for everyf ∈L2([0,1)), Parseval’s identity holds:
Z 1 0
f2(t)dt=
∞
X
i=0
fi2kφi(t)k2, (5) where
fi = 1 h
Z 1 0
f(t)φi(t)dt, i= 0,1, . . . , m−1. (6)
An arbitrary function can be expanded in vector form as:
f(t)'FTΦ(t) = ΦT(t)F. (7) whereF = [f0, f1, ..., fm−1]T and Φ(t) = [φ0(t), φ1(t), ..., φm−1(t)]T.
Now let k(x, t) be arbitrary L2 function of two variables on [0,1)×[0,1).
It can be expanded by BPFs as:
k(x, t)'ΦT(x)KΨ(t), (8)
where Φ(x) and Ψ(t) are m1 and m2 dimensional BPF vectors, respectively, and K is the (m1 ×m2) Block Pulse coefficient matrix with entries kij, i = 0,1, ..., m1−1, j = 0,1, ..., m2−1, as follows:
kij =m1m2 Z 1
0
Z 1 0
k(x, t)φi(x)φj(t)dx dt. (9) In this paper for convenience, we put m1 =m2 =m.
Moreover, from the disjointness property of BPFs, follows:
Φ(t)ΦT(t) =
φ0(t) 0 · · · 0 0 φ1(t) · · · 0 ... ... . .. ... 0 0 · · · φm−1(t)
, (10)
Z 1 0
Φ(t)ΦT(t)dt= h I, (11)
ΦT(t)Φ(t) = 1, (12)
Φ(t)ΦT(t)V = ˜VΦ(t), (13)
whereV is an m dimensional vector and ˜V =diag(V). Moreover:
ΦT(t)BΦ(t) = ΦT(t)B,b (14) where B is an (m×m) matrix, ˆB is an m dimensional vector with elements equal to the diagonal entries of matrixB. As represented in [2, 9]:
Z t 0
Φ(τ)dτ 'PΦ(t), (15)
whereP is the following (m×m) matrix:
P = h 2
1 2 2 · · · 2 0 1 2 · · · 2 0 0 1 · · · 2 ... ... ... . .. ...
0 0 0 · · · 1
, (16)
called operational matrix of integration. So, the integral of every function f can be approximated as follows:
Z t 0
f(τ)dτ 'FTPΦ(t). (17)
3 A Method to Solve VFIDE
Consider the following VFIDE:
f0(x) = g(x) + Z x
0
k(x, t)f(t)dt+ Z 1
0
l(x, t)f(t)dt 0≤x, t <1, f(0) =a0,
(18) First, we integrate from (3.1) in the interval [0, x], we will have:
f(x) =a0+ Z x
0
g(u)du+ Z x
0
Z u 0
k(u, t)f(t)dtdu+
Z x 0
Z 1 0
l(u, t)f(t)dtdu. (19) Now, we approximate functionsf, g,k and l with respect to BPFs as follows:
f(x)'FTΦ(x) = ΦT(x)F, g(x)'ΦT(x)G,
k(x, t)'ΦT(x)KΦ(t), l(x, t)'ΦT(x)LΦ(t),
(20)
where vectorsF andG and matrices K andL are BPFs coefficients off, g,k and l, respectively.
Substituting (3.3) into (3.2) we will have:
ΦT(x)F 'a0+ ΦT(x)PTG+ Z x
0
ΦT(u)KZ u 0
Φ(t)ΦT(t)F dt du+
Z x 0
ΦT(u)LZ 1 0
Φ(t)ΦT(t)dt
F du. (21) Substituting (2.10), (2.11), (2.12) and (2.14) into (3.4), results:
ΦT(x)F 'a0ΦT(x) Φ(x) + ΦT(x)PT G+
Z x 0
ΦT(u)KF P˜ Φ(u)du+hΦT(x)PTLF. (22) PutB =KF P˜ and apply (2.13), we obtain:
Z x 0
ΦT(u)BΦ(u)du= Z x
0
ΦT(u)Bdub 'ΦT(x)PTB,b puttingAb=PTBb and using (2.13), above relation converts to:
ΦT(x)PTBb = ΦT(x)Ab= ΦT(x)AΦ(x), (23) whereA=diag(A). Substituting (3.6) into (3.5), results:b
ΦT(x)F 'ΦT(x)(a0I+A)Φ(x) + ΦT(x)(PTG+hPTLF).
Now, withC =a0I+A and using(2.13), we have:
ΦT(x)F 'ΦT(x)Cb+ ΦT(x)(PTG+hPTLF), or
F 'Cb+PTG+hPTLF, (24)
whereCb can be written as follows:
Cb =a01+h2K0F (25)
where1= (1,1,· · ·,1)T and
K0 =
1
4k11 0 0 · · · 0 0
2
X
i=1 00
ki1 14k22 0 · · · 0 0
3
X
i=1 00
ki1
3
X
i=2 00
ki2 14k33 · · · 0 0 ... ... ... . .. ... ...
m
X
i=1 00
ki1
m
X
i=2 00
ki2
m
X
i=3 00
ki3 · · ·
m
X
i=m−1 00
ki,m−1 14kmm
and X00
means that the first and last terms have factor 12.
Substituting (3.8) into (3.7) and replacing ' with =, we will have:
(I−hPTL−h2K0)F = (PTG+a01). (26) Solving the system of equations (3.9) we obtain m unknowns f0, f1,· · ·, fm−1. Also, structure ofK0 shows that we do not need to evaluate kij for j > i.
4 Error Analysis
In this section, we analyse the error when a differentiable functionf(x) is rep- resented in a series of BPFs over the intervalI = [0,1). We need the following theorem.
Theorem: Suppose f is continuous in I, is differentiable in (0,1), and there is a number M such that |f0(x)| ≤M, for every x∈I. Then
|f(b)−f(a)| ≤M|b−a|, for all a, b∈I.
Proof. See [11].
Now, we assume that f(x) is a differentiable function on I such that
|f0(x)| ≤M. We define the error between f(x) and its BPFs expansion over every subinterval Ii as follows:
ei(x) =fi−f(x), x∈Ii whereIi = [mi,i+1m ).
It can be shown that:
keik2 = Z i+1m
i m
ei2(x)dx = Z i+1m
i m
fi−f(x)
2dx= 1 m
fi−f(η)
2, η ∈Ii,
(27) where we used mean value theorem for integral. Using Eq.(2.5) and the mean value theorem, we have:
fi =m Z i+1m
i m
f(x)dx=m 1
mf(ζ) =f(ζ), ζ ∈Ii. (28) Substituting (4.2) into (4.1) and using Theorem, we will have:
keik2 = 1 m
f(ζ)−f(η)
2 ≤ M2
m |ζ−η|2 ≤ M2
m3. (29)
This leads to:
ke(x)k2 = Z 1
0
e2(x)dx= Z 1
0
m−1X
i=0
ei(x)
2dx
= Z 1
0
m−1X
i=0
ei2(x)
dx+ 2 X
i≤j
Z 1 0
ei(x)ej(x)dx.
Since fori6=j, Ii∩Ij =∅, then ke(x)k2 =
m−1
X
i=0
Z 1 0
ei2(x)dx
=
m−1
X
i=0
keik2. (30)
Substituting (4.3) into(4.4), we have
ke(x)k2 ≤ M2 m2,
hence,ke(x)k=O(m1), where e(x) = fm(x)−f(x) and fm(x) =
m−1
X
i=0
fiφi(x).
5 Numerical Examples
In this section, we present some examples and their numerical results. Ex- amples 2 and 3 were used in [13] and [1], respectively. All computations are performed by the Maple 9.5 software package.
Example 1: Consider the following VIDE:
f0(x) = Z x
0
f(t)dt, f(0) = 1
with the exact solution f(x) = cosh(x). Indeed, in this example, we have g(x) = 0, l(x, t) = 0 and k(x, t) = 1.
The numerical results are shown in table 1.
Table 1: Numerical results for Example 1
x Exact solution Approximate solution Approximate solution
m = 32 m = 64
0 1 1.000244 1.000061
0.1 1.005004 1.006111 1.005193
0.2 1.020067 1.020829 1.019166
0.3 1.045339 1.044527 1.046811
0.4 1.081072 1.077413 1.080467
0.5 1.127626 1.136067 1.131772
0.6 1.185465 1.191664 1.186506
0.7 1.255169 1.257743 1.251676
0.8 1.337435 1.334886 1.341668
0.9 1.433086 1.423773 1.431547
Example 2: Consider the following VIDE [13]:
f0(x) = 1− Z x
0
f(t)dt,
f(0) = 0
with the exact solution f(x) = sin(x). In this example, we have g(x) = 1, l(x, t) = 0 and k(x, t) = −1.
Table 2: Numerical results for Example 2
x Exact solution Approximate solution Approximate solution
m = 64 m = 128
0 0 0.007812 0.003906
0.1 0.099833 0.101383 0.097500
0.2 0.198669 0.194063 0.197901
0.3 0.295520 0.299980 0.296263
0.4 0.389418 0.387959 0.391571
0.5 0.479426 0.486243 0.482844
0.6 0.564642 0.565904 0.562700
0.7 0.644218 0.640595 0.646312
0.8 0.717356 0.720580 0.717892
0.9 0.783327 0.782319 0.784773
The numerical results are shown in table 2.
Example 3: Consider the following linear FIDE [1]:
f0(x) =−e−x+e−1−1 + Z 1
0
f(t)dt, f(0) = 1
with the exact solutionf(x) =e−x. We haveg(x) = −e−x+e−1−1,l(x, t) = 1 and k(x, t) = 0. This example was used in [1] by using spline functions with m = 1 and m = 2, where m equals the order of spline functions. The abso- lute values of error in [1] with stepsize h = 0.1 at the point 0.4 are given as 1.5×10−2 and 1.8×10−3 for m= 1 and m= 2, respectively.
The numerical results are shown in table 3.
Example 4: Consider the following linear VFIDE:
f0(x) =−2sin(x)−x2sin(2x) + 2sin(2x)−2xcos(2x)−2ex +5ex−1+ 2x+
Z x 0
cos(x+t)f(t) dt+ Z 1
0
ex−tf(t)dt, f(0) = 0
with the exact solutionf(x) =x2. In this example, we haveg(x) = −2sin(x)−
x2sin(2x) + 2sin(2x) − 2xcos(2x) − 2ex + 5ex−1 + 2x, l(x, t) = ex−t and
Table 3: Numerical results for Example 3
x Exact solution Approximate solution Approximate solution
m = 16 m = 64
0 1 0.969719 0.992248
0.1 0.904837 0.910993 0.903455
0.2 0.818731 0.804005 0.822608
0.3 0.740818 0.755324 0.737384
0.4 0.670320 0.666636 0.671399
0.5 0.606530 0.588375 0.601842
0.6 0.548812 0.552766 0.547986
0.7 0.496585 0.487894 0.498951
0.8 0.449329 0.458378 0.447261
0.9 0.406570 0.404606 0.407240
k(x, t) = cos(x+t).
The numerical results are shown in table 4.
Table 4: Numerical results for Example 4
x Exact solution Approximate solution Approximate solution
m = 32 m = 64
0 0 0.000493 0.000123
0.1 0.01 0.012240 0.010384
0.2 0.04 0.041571 0.038226
0.3 0.09 0.088486 0.092926
0.4 0.16 0.152985 0.158856
0.5 0.25 0.266336 0.257994
0.6 0.36 0.371862 0.362012
0.7 0.49 0.494971 0.483610
0.8 0.64 0.635662 0.647692
0.9 0.81 0.793932 0.807378
6 Conclusion
The proposed method for solving linear VFIDE was based on BPFs and their operational matrix. This method converts equation to a linear system of algebraic equations. It,s accuracy is shown on some examples. As examples
show error will decrease asm increases. The method is easy to apply and it,s accuracy is comparable with other methods.
Acknowledgements
The authors are highly grateful to the Islamic Azad University, Takestan Branch, Iran, for giving all types of support in conducting this research
References
[1] A. Ayad, Spline approximation for first-order Fredholm integro- differential equation, Studia Universitatis Babes-Bolyai, Mathematica, 41(3) (1996), 1-8.
[2] E. Babolian and Z. Masouri, Direct method to solve Volterra integral equation of the first kind using operational matrix with block-pulse func- tions, J. Comput. Appl. Math., 220(2008), 51-57.
[3] G.M. Phillips, Analysis of numerical iterative methods for solving integral and integro-differential equations,Comput. J., 13(1970), 297-300.
[4] G. Micala and G.F. Weather, Direct numerical spline methods for first order Fredholm integro-differential equation,Analyse Numerical et theorie de L,Approximation, Revue D, 53(1) (1993), 59-66.
[5] G.P. Rao, Piecewise Constant Orthogonal Functions and Their Applica- tion to Systems and Control, Springer-Verlag, (1983).
[6] K. Maleknejad, F. Mirzaee and S. Abbasbandy, Solving linear integro differential equations system by using rationalized Haar functions method, Appl. Math. Comput., 155(2004), 317-328.
[7] K. Maleknejad and M. Hadizadeh, New computational method for Volterra Fredholm integral equations, Computers and Mathematics with Applications, 37(1999), 1-8.
[8] K. Maleknejad and M.R. Fadaei Yami, A computational method for system of Volterra Fredholm integral equations, Appl. Math. Comput., 183(2006), 589-595.
[9] K. Maleknejad and Y. Mahmoudi, Numerical solution of linear Fredholm integral equation by using hybrid Taylor and Block Pulse functions,Appl.
Math. Comput., 149(2004), 799-806.
[10] K. Maleknejad and Y. Mahmoudi, Taylor polynomial solution of high order nonlinear Volterra Fredholm integro-differential equations, Appl.
Math. Comput., 145(2003), 641–653.
[11] R.G. Bartle,The Elements of Real Analysis,(Second edition), John Wiley and Sons, (1976).
[12] S.H. Behiry and H. Hashish, Wavelet methods for the numerical solution of Fredholm integro-differential equations, International Journal of Applied Mathematics, 11(1) (2002), 27-35.
[13] S. Shahmorad, Numerical solution of the general from linear Fredholm- Volterra integro-differential equations by the Tau method with an error es- timation, Applied Mathematics and Computation, 167(2005), 1418–1429.
[14] W. Volk, The numerical solution of linear integro-differential equations by projection methods,J. InEq., 9(1985), 171-190.
[15] Z.H. Jiang and W. Schaufelberger, Block Pulse Functions and Their Ap- plications in Control Systems, Springer-Verlag, (1992).
