Sinc numerical solution for the Volterra integro-differential equation

Commun Nonlinear Sci Numer Simulat 15 (2010) 700–706

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Sinc numerical solution for the Volterra integro-differential equation M. Zarebnia * Department of Mathematics, University of Mohaghegh Ardabili, 56199-11367 Ardabil, Iran

a r t i c l e

i n f o

Article history: Received 12 June 2008 Received in revised form 16 March 2009 Accepted 17 April 2009 Available online 3 May 2009 PACS: 02.60.Nm 02.60.x 04.25.g

a b s t r a c t In this paper, numerical solution of Volterra integro-differential equation by means of the Sinc collocation method is considered. Convergence  pffiffiffi  analysis is given, it is shown that the Sinc solution produces an error of order O ek N where k > 0 is a constant. This approximation reduces the Volterra integro-differential equation to an explicit system of algebraic equations. The method is applied to a few test examples to illustrate the accuracy and the implementation of the method. Ó 2009 Elsevier B.V. All rights reserved.

Keywords: Volterra integro-differential equation Sinc method Exponential convergence

1. Introduction We consider the following Volterra integro-differential equation of the form:

u0 ðxÞ ¼ k

Z

x

kðx; tÞuðtÞdt þ lðxÞuðxÞ þ f ðxÞ;

x; t 2 C ¼ ½a; b;

a

uðaÞ ¼ ua ;

ð1Þ

where the functions f ðxÞ, lðxÞ and the kernel kðx; tÞ are known and uðxÞ is the solution to be determined [1]. Several numerical methods for approximating the linear Volterra integro-differential equation are known. Ref. [2] presented a solution to the Volterra integro-differential equations by the Tau method with arbitrary polynomial bases. Berenguer et al. [3] introduced a numerical method to approximate the solution of the linear Volterra integro-differential equation. By making use of the expression of a function of the Banach space C½0; 1 in terms of a Schauder basis, they defined a sequence of functions which approximate (in the uniform sense) the solution of such equation. A method for the solution of Volterra integro-differential equations by using single-term Walsh series has been proposed by Sepehrian et al. [4]. Sinc methods have increasingly been recognized as powerful tools for solving a wide range of linear and nonlinear problems arising from scientific and engineering applications including, heat transfer [5], population growth [6], fluid mechanics [7], inverse problems [8] and medical imaging [9]. In particular, they have become very popular in solving initial and boundary value problems of ordinary and partial differential equations including those with Dirichlet, Neuman and other boundary conditions [10,11]. Sinc methods has also been employed as forward solvers in the solution of integral equations [12,13]. There are several advantages to using approximations based on Sinc numerical methods. Unlike most numerical techniques, * Tel.: +98 4515514702; fax: +98 4515514701. E-mail address: [email protected] 1007-5704/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2009.04.021

701

M. Zarebnia / Commun Nonlinear Sci Numer Simulat 15 (2010) 700–706

it is now well-established that they are characterized by exponentially decaying errors [14], and also, they are highly efficient and adaptable in handling problems with singularities [15]. Finally, due to their rapid convergence, Sinc numerical methods do not suffer from the common instability problems associated with other numerical methods [16]. In this paper a global approximation for the solution of the Eq. (1) using the Sinc functions is developed. The method consists of reducing the solution of Eq. (1) to a set of algebraic equations. The properties of Sinc function are then utilized to evaluate the unknown coefficients. The layout of the paper is as follows. In Section 2, we briefly review the concept and some properties of the Sinc function and Sinc method that are necessary for the formulation of the discrete system. Section 3, is concerned with the Sinc collocation discretizations for Volterra integro-differential equation. In Section 4, the convergence analysis of the method has been discussed. It is shown that the Sinc procedure converges to the solution at an exponential rate. Finally, numerical experiments are given in Section 5 to verify the exponential convergence rate, and to demonstrate the efficiency and accuracy of the proposed numerical scheme. 2. A survey of some properties of the Sinc function The goal of this section is to recall notations and definition of the Sinc function, state some known results, and derive useful formulas that are important for this paper. These are discussed thoroughly in [17]. The Sinc function is defined on the whole real line by

(

SincðxÞ ¼

sinðpxÞ px ;

x–0;

1;

x ¼ 0:

ð2Þ

For any h > 0, the translated Sinc functions with evenly spaced nodes are given as follows:

  x  jh Sðj; hÞðxÞ ¼ Sinc ; h

j ¼ 0; 1; 2; . . .

ð3Þ

which are called the jth Sinc functions. The Sinc function form for the interpolating point xk ¼ kh is given by ð0Þ

Sðj; hÞðkhÞ ¼ djk ¼ Let ð1Þ

dkj

¼

1 þ 2

Z

kj

0



1; k ¼ j;

ð4Þ

0; k–j:

sinðptÞ dt; pt

ð5Þ ð1Þ

then define a matrix whose ðk; jÞth entry is given by dkj series

Cðu; hÞðxÞ ¼

ð1Þ

as Ið1Þ ¼ ½dkj . If u is defined on the real line, then for h > 0 the

  x  jh uðjhÞSinc h j¼1 1 X

ð6Þ

is called the Whittaker cardinal expansion of u, whenever this series converges, u is approximated by using the finite number of terms in (6). For positive integer N, we defined

C N ðu; hÞðxÞ ¼

  x  jh : uðjhÞSinc h j¼N N X

ð7Þ

For further explanation of the procedure, we consider the following definitions and theorem in [17]. Definition 1. Let D be a simply-connected domain in the complex plane(z ¼ x þ iy) having boundary @D. Let a and b denote two distinct points of @D and / denote a conformal map of D onto Dd , where Dd denote the region fw 2 C : jIwj < d; d > 0g such that /ðaÞ ¼ 1 and /ðbÞ ¼ 1. Let w ¼ /1 denote the inverse map, and let C be defined by C ¼ fz 2 C : z ¼ wðuÞ; u 2 Rg. Given /, w and a positive number h, let us set zk ¼ zk ðhÞ ¼ wðkhÞ; k ¼ 0; 1; 2; . . ., let us also define q by qðzÞ ¼ e/ðzÞ . Definition 2. Let La ðDÞ be the set of all analytic functions u, for which there exists a constant, C, such that

juðzÞj 6 C

jqðzÞja ½1 þ jqðzÞj2a

;

z 2 D;

0 < a 6 1:

ð8Þ ð1Þ

u Theorem 1. Let /0 2 La ðDÞ, with 0 < a 6 1, and 0 < d 6 p, let dkj constant, c1 , which is independent of N, such that

 Z   zk N 1 X  ð1Þ uðzj Þ  2 uðtÞdt  h dkj  6 c1 eðpdaNÞ :  0   a ðz Þ / j j¼N

1

be defined as in (5), and let h ¼ ðapNd Þ2 . Then there exists a

ð9Þ

702

M. Zarebnia / Commun Nonlinear Sci Numer Simulat 15 (2010) 700–706

3. The approximate solution of Volterra integro-differential equation Let uðxÞ be the exact solution of the integral Eq. (1) and uðxÞ 2 La ðDÞ. By considering Eq. (1) and integrating (1) from a to x, we get:

uðxÞ ¼

Z

x

 Z k

a

n

kðn; tÞuðtÞdt þ lðnÞuðnÞ þ f ðnÞ dn þ uðaÞ;

x 2 C ¼ ½a; b:

ð10Þ

a

For simplicity, we write

Z

KðnÞ ¼

n

kðn; tÞuðtÞdt:

ð11Þ

a

l K Now, we let /0 2 La ðDÞ, /0 u 2 La ðDÞ and /0f 2 La ðDÞ. By setting x ¼ xk ; k ¼ N; . . . ; N and applying Theorem 1 for the right-hand side of (10), we obtain; N X

uðxk Þ ’ kh

ð1Þ

dk;l

l¼N

N N X X Kðxl Þ ð1Þ lðxl Þ ð1Þ f ðxl Þ þh uðxl Þ þ h þ ua : dk;l dk;l 0 0 ðx Þ / ðxl Þ / /0 ðxl Þ l l¼N l¼N

ð12Þ

For the first term on the right-hand side of above relation, by considering relation (11) and using Theorem 1, we have;

h

N X

ð1Þ

dk;l

l¼N

N N X X Kðxl Þ kðxl ; t j Þ 2 ð1Þ ð1Þ ’h uðt j Þ: dk;l dl;j 0 0 / ðxl Þ / ðxl Þ/0 ðt j Þ l¼N j¼N

ð13Þ

Having replaced the first term on the right-hand side of (12) with the Eq. (13), we obtain; N N X X

2

uk  kh

ð1Þ ð1Þ

dk;l dl;j

l¼N j¼N

where uk ¼ uðxk Þ, kl;j ¼ kðxl ; tj Þ,

x  a ; /ðxÞ ¼ ln bx /0 ðxÞ ¼

N N X X kl;j ð1Þ ll ð1Þ fl dk;l dk;l þ ua : 0 0 uj  h 0 ul ’ h /l /j / /0l l l¼N l¼N

ð14Þ

ll ¼ lðxl Þ, fl ¼ f ðxl Þ and

/ðaÞ ¼ 1;

/ðbÞ ¼ þ1; kh

ba ; ðx  aÞðb  xÞ

xk ¼ wðkhÞ ¼ /1 ðkhÞ ¼

a þ be : 1 þ ekh

There are ð2N þ 1Þ unknowns uj ; j ¼ N; N þ 1; . . . ; N  1; N; to be determined in (14). In order to determine these ð2N þ 1Þ unknowns, we rewrite these system which is the linear system of equations in matrix form. Corresponding to a given funcðmÞ

ðmÞ

tion u defined on C, we use the notation DðuÞ ¼ diag ðuðxN Þ; . . . ; uðxN ÞÞ. We set IðmÞ ¼ ½dkj ; m ¼ 1; 0, where dkj denotes the ðk; jÞth element of the matrix I

ðmÞ

. Now, we can simplify the system (14) in the matrix form

AU ¼ P;

ð15Þ

where

A ¼ Ið0Þ  hI

ð1Þ



D

 

l

/0

 kh

2

      1 1 e Ið1Þ D 0 Ið1Þ  KD ; / /0



e ¼ kðxi ; t j Þ ; i; j ¼ N; . . . ; N; K " #T   N N X X ð1Þ fl ð1Þ fl dN;l 0 þ ua ; . . . ; h dN;l ð 0 Þ þ ua ; P¼ h /l /l l¼N l¼N U ¼ ½uN ; . . . ; uN T : The notation ‘‘” denotes the Hadamard matrix multiplication. Having used the solution uj ; j ¼ N; N þ 1; . . . ; N, in the system (15), we employ a method similar to the Nyström’s idea for the Volterra integro-differential equation [17], i.e., we use

uN ðxÞ ¼ kh

2

N N X X

N N X X kðnj ; tl Þ lðnl Þ f ðnl Þ Xh;l ðnj ÞXh;j ðxÞul þ h Xh;l ðxÞul þ h Xh;l ðxÞ þ ua ; 0 0 / ðnl Þ /0 ðnl Þ / ðnj Þ/ ðtl Þ l¼N l¼N 0

l¼N j¼N

where

Xh;l ðxÞ ¼

1 þ 2

Z a

x

Sðl; hÞ  /ðtÞdt:

ð16Þ

M. Zarebnia / Commun Nonlinear Sci Numer Simulat 15 (2010) 700–706

703

4. Error analysis Now, we discuss the convergence of the Sinc method for the linear Volterra integro-differential Eq. (1). For each N, we can find uj which is the solution of the linear system (15), also by using uj we obtain the approximate solution uN ðxÞ. In order to e  P, where U e is a vector defined by derive a bound for juðxÞ  uN ðxÞj we need to estimate the norm of the vector A U

e ¼ ðuðxN Þ; . . . ; uðxN ÞÞT ; U where uðxj Þ is the value of the exact solution of integro-differential equation at the Sinc points xj . For this purpose we need the following lemma. 1 e be the exact solution of the given integral Eq. (1) and let h ¼ pd 2 ; kðx;Þ Lemma 1. Let U 2 La ðDÞ for all x 2 C. Then there exists a aN /0 constant c2 independent of N such that 1

1

e  Pk 6 c2 N2 expfðpdaNÞ2 g; kA U

ð17Þ

the norm ðk  kÞ is norm two. e  P. Using the Theorem 1 with the optimal mesh Proof. We derive a bound for the kth component v k of the vector m ¼ A U 1 size h ¼ apNd 2 and by the assumption on the kernel, we have the following bound on v k :

e  PÞ j jv k j ¼ jðA U k   N N X X kðxl ; t j Þ  2 ð1Þ ð1Þ uðt j Þ dk;l dl;j ¼ uðxk Þ  kh 0  / ðxl Þ/0 ðt j Þ l¼N j¼N h

N X

ð1Þ dk;l

l¼N

lðxl Þ /0 ðxl Þ

uðxl Þ  h

N X

ð1Þ dk;l

l¼N

  f ðxl Þ   ua  0  / ðxl Þ

1 2

6 c3 expfðpdaNÞ g; and it thus follows that

e  Pk ¼ kA U

N X

!12

1

1

jv k j2

6 c2 N2 expfðpdaNÞ2 g:



k¼N

We can also obtain a bound on the error uðxÞ  uN ðxÞ in the maximum norm, where uðxÞ is the exact solution and uN  ðxÞ is the  1=2 approximate solution (16). From Lemma 1 we can show that the Sinc method converges at rate of O expðkN Þ , where k > 0. Theorem 2. Let us consider all assumptions of Lemma 1. Let U N ðxÞ be the approximate solution of integro-differential Eq. (1) given by (16). Then there exists a constant c4 independent of N such that 1

1

sup juðxÞ  uN ðxÞj 6 c4 lN2 expfðpdaNÞ2 g;

ð18Þ

x2ða;bÞ

where

l ¼ kA1 k.

Proof. By considering the given Eq. (1) and all assumptions we obtain

 Z  juðxÞ  uN ðxÞj ¼ k

x a

Z

n

kðn; tÞuðtÞdtdn þ a 2

 kh

Z

x

lðnÞuðnÞdn þ

a N N X X l¼N j¼N

kðnj ; t l Þ Xh;l ðnj ÞXh;j ðxÞul /0 ðnj Þ/0 ðt l Þ

Z

x

f ðnÞdn

a

  N X f ðnl Þ  X ðxÞ  h;l 0 0  ðn Þ ðn Þ / / l l l¼N l¼N    N N  X  X kðnj ; t l Þ  2 6 kh Xh;l ðnj ÞXh;j ðxÞjuðtl Þ  ul j  0 0   ðn Þ/ ðt Þ / l j l¼N j¼N   N X  lðnl Þ    þh /0 ðn Þ Xh;l ðxÞjuðnl Þ  ul j ¼ EN ðSayÞ: l l¼N

h

Note that

N X

lðnl Þ

Xh;l ðxÞul  h

704

M. Zarebnia / Commun Nonlinear Sci Numer Simulat 15 (2010) 700–706

0

2 112   N X N X kðnj ; t l Þ   Xh;l ðnj ÞXh;j ðxÞ A 6 S1 ;  0 0   / ðnj Þ/ ðt l Þ l¼N j¼N

2@

kh

h

2 !12 N  X   lðnl Þ Xh;l ðxÞ 6 S2  /0 ðn Þ l

l¼N

hold for x 2 ½a; b, then using the Schwarz inequality, we get

0 EN 6 kh

2@

2 112  !12  N X N N X kðnj ; t l Þ  A X  2 X ðn Þ X ðxÞ uðt Þ  u j j   h;l j h;j l l   /0 ðnj Þ/0 ðtl Þ l¼N j¼N l¼N

!12 2 !12 X N  N X  lðnl Þ 2   þh juðnl Þ  ul j /0 ðn Þ Xh;l ðxÞ l¼N

l

l¼N

e  Uk; 6 Sk U where S ¼ S1 þ S2 . Since from (15)

U ¼ A1 P; then

e  Uk ¼ k U e  A1 Pk 6 kA1 kkA U e  Pk kU

ð19Þ

holds, we have from Lemma 1

e  Pk EN 6 SkA1 kkA U 1

1

6 Slc2 N 2 expfðpdaNÞ2 g:

ð20Þ

Therefore from the above relation we conclude that 1

1

sup juðxÞ  uN ðxÞj 6 c4 lN2 expfðpdaNÞ2 g:



x2ða;bÞ

5. Numerical examples In order to illustrate the performance of the Sinc method in solving Volterra integro-differential equation and justify the accuracy and efficiency of the method, we consider the following examples. The solutions of the given examples are obtained for different values of a; 0 < a 6 1 and also for different values of N. Let uðxj Þ denote the exact solution of the given examples, and let uN ðxj Þ be the computed solutions by the presented method. The error is reported on set of the Sinc grid points

S ¼ fxN ; . . . ; x0 ; . . . ; xN g; ð21Þ

kh

xk ¼

a þ be ; 1 þ ekh

k ¼ N; . . . ; N:

The maximum absolute error on the grid points Sinc is defined as

  kEs ðhÞk ¼ max uðxj Þ  un ðxj Þ:

ð22Þ

N6j6N

Table 1 Results for Example 1. N 5 10 20 30 40 50 60 70 80 90 100

h 1.404963 0.993459 0.702481 0.573574 0.496729 0.444288 0.405578 0.375492 0.351241 0.331153 0.314159

condðAÞ ¼ kAkkA1 k

kEs ðhÞk1 3

4:61435  10 1:95377  104 6:54692  106 4:49157  107 3:72779  108 3:73076  109 4:08840  1010 6:00222  1011 1:46743  1011 3:66818  1012 9:31033  1013

11.195 21.2442 41.3027 61.3445 81.3789 101.409 121.436 141.46 161.483 181.504 201.525

705

M. Zarebnia / Commun Nonlinear Sci Numer Simulat 15 (2010) 700–706 Table 2 Results for Example 2. N

h

5 10 20 30 40 50 60 70 80 90 100

condðAÞ ¼ kAkkA1 k

kEs ðhÞk1 3

1.404963 0.993459 0.702481 0.573574 0.496729 0.444288 0.405578 0.375492 0.351241 0.331153 0.314159

4:61435  10 2:92118  104 6:95317  106 3:49034  107 2:71282  108 2:78475  109 3:49452  1010 5:12093  1011 8:47589  1012 1:55120  1012 3:09086  1013

26.8193 45.0391 77.3449 107.082 135.418 162.835 189.585 215.826 241.658 267.153 292.365

Table 3 Results for Example 3. N

h

5 10 20 30 40 50 60 70 80 90 100

condðAÞ ¼ kAkkA1 k

kEs ðhÞk1 3

1.404963 0.993459 0.702481 0.573574 0.496729 0.444288 0.405578 0.375492 0.351241 0.331153 0.314159

1:75196  10 1:05678  104 1:84725  106 8:03194  108 5:87463  109 5:80425  1010 7:09318  1011 1:02003  1011 1:67444  1012 3:08975  1013 6:261666  1014

14.3287 25.8357 48.0291 69.733 91.1767 112.452 133.607 154.671 175.662 196.593 217.474

The maximum absolute errors in numerical results are tabulated in Tables 1–3. Example 1. Consider the linear Volterra integro-differential equation

u0 ðxÞ ¼

Z

x

xð1 þ 2xÞetðxtÞ uðtÞdt  uðxÞ þ 1 þ 2x;

6x61

0

uð0Þ ¼ 1; 2

with the exact solution uðxÞ ¼ ex [3]. This problem is solved by the presented method. The approximate solution is calculated for different values of N, a ¼ 12, 1 and d ¼ p2 , which yield the optimal Sinc mesh size h ¼ pðN1 Þ2 . Maximum absolute errors in numerical solution of Example 1 are tabulated in Table 1. These results show the efficiency and applicability of the Sinc method. Example 2. We consider the integro-differential equation

u0 ðxÞ ¼

Z

x

exþt uðtÞdt þ uðxÞ þ 0

ex  e3x ; 2

uð0Þ ¼ 1; with the exact solution uðxÞ ¼ ex . 1 We solved Example 2 by the Sinc method for different values of N and h ¼ pðN1 Þ2 , since d ¼ p2 and a ¼ 12. The maximum absolute errors on the Sinc grid S is tabulated in Table 2. Example 3. Consider the Volterra integro-differential equation

u0 ðxÞ ¼

Z

x 0

x  x 1 uðtÞdt þ uðxÞ  lnð1 þ xÞ lnð1 þ xÞ þ 1 þ ; tþ1 2 1þx

06x61

uð0Þ ¼ 0; with the exact solution uðxÞ ¼ lnðx þ 1Þ. 1

The approximate solution is calculated for different values of N and the optimal Sinc mesh size h ¼ p N1 2 d ¼ p2 ; a ¼ 12 . The maximum absolute errors in computed solutions are tabulated in Table 3. These results show the accuracy and efficiency of the Sinc method.

706

M. Zarebnia / Commun Nonlinear Sci Numer Simulat 15 (2010) 700–706

Acknowledgements The author thank the referees for many valuable comments in the revision of the paper. References [1] Delves LM, Mohamed JL. Computational methods for integral equations. Cambridge: Cambridge University Press; 1985. [2] Pour-Mahmoud J, Rahimi-Ardabili MY, Shahmorad S. Numerical solution of Volterra integro-differential equations by the Tau method with the Chebyshev and Legendre bases. Appl Math Comput 2005;170:314–38. [3] Berenguer MI, Fortes MA, Garralda Guillem AI, Ruiz Galán M. Linear Volterra integro-differential equation and Schauder bases. Appl Math Comput 2004;159:495–507. [4] Sepehrian B, Razzaghi M. Single-term Walsh series method for the Volterra integro-differential equations. Eng Anal Boundary Elem 2004;28:1315–9. [5] Narasimhan S, Chen K, Stenger F. A Harmonic-Sinc solution of the laplace equation for problems with singularities and semi-infinite domains. Numer Heat Transfer B 1998;33:433–50. [6] Al-Khaled K. Numerical approximations for population growth models. Appl Math Comput 2005;160:865–73. [7] Winnter DF, Bowers KL, Lund J. Wind-driven currents in a sea with variable eddy viscosity calculated via a Sinc–Galerkin technique. J Numer Methods Fluids 2000;33:1041–73. [8] Smith R, Bowers K. A Sinc–Galerkin estimation of diffusivity in parabolic problems. Inverse Probl 1993;9:113–45. [9] Stenger F, O’Reilly MJ. Computing solution to medical problems via Sinc convolution. IEE Trans Automat Control 1998;43:843–8. [10] Smith RC, Bogar GA, Bowers KL, Lund J. The Sinc–Galerkin method for fourth-order differential equations. SIAM J Numer Anal 1991;28:760–88. [11] Smith RC, Bowers KL, Lund J. Fully Sinc–Galerkin method for Euler–Bernoulli beam models. Numer Methods Partial Differential Equations 1992;8:171–202. [12] Rashidinia J, Zarebnia M. Solution of a Volterra integral equation by the Sinc-collocation method. J Comput Appl Math 2007;206:801–13. [13] Rashidinia J, Zarebnia M. Convergence of approximate solution of system of Fredholm integral equations. J Math Anal Appl 2007;333:1216–27. [14] Keinert F. Uniform approximation jxjb by Sinc functions. J Approx Theory 1991;66:44–52. [15] Stenger F. Sinc pack-summary of basic Sinc methods. Salt Lake City, UT: Department of Computer Science, University of Utah; 1995. [16] Sababheh MS, Al-Khaled AMN. Some convergence results on Sinc interpolation. J Inequal Pure Appl Math 2003;4:32–48. [17] Stenger F. Numerical methods based on Sinc and analytic functions. New York: Springer Verlag; 1993.