[SOLVED] scheme algorithm This article appeared in a journal published by Elsevier. The attached

This article appeared in a journal published by Elsevier. The attached
copy is furnished to the author for internal non-commercial research
and education use, including for instruction at the authors institution

and sharing with colleagues.

Other uses, including reproduction and distribution, or selling or
licensing copies, or posting to personal, institutional or third party

websites are prohibited.

In most cases authors are permitted to post their version of the
article (e.g. in Word or Tex form) to their personal website or
institutional repository. Authors requiring further information

regarding Elseviers archiving and manuscript policies are
encouraged to visit:

http://www.elsevier.com/copyright

http://www.elsevier.com/copyright

Authors personal copy

Journal of Computational and Applied Mathematics 233 (2009) 938950

Contents lists available at ScienceDirect

Journal of Computational and Applied
Mathematics

journal homepage: www.elsevier.com/locate/cam

Spectral methods for weakly singular Volterra integral equations with
smooth solutions
Yanping Chen a,, Tao Tang b
a School of Mathematical Sciences, South China Normal University, Guangzhou 510631, China
b Department of Mathematics, Hong Kong Baptist University, Kowloon Tong, Hong Kong & Faculty of Science, Beijing University of Aeronautics and Astronautics,
Beijing, China

a r t i c l e i n f o

Article history:
Received 23 July 2008
Received in revised form 19 May 2009

MSC:
35Q99
35R35
65M12
65M70

Keywords:
Volterra integral equations
Weakly singular kernels
Smooth solutions
Spectral methods
Error analysis

a b s t r a c t

Wepropose and analyze a spectral Jacobi-collocation approximation for the linear Volterra
integral equations (VIEs) of the second kind with weakly singular kernels. In this work,
we consider the case when the underlying solutions of the VIEs are sufficiently smooth.
In this case, we provide a rigorous error analysis for the proposed method, which shows
that the numerical errors decay exponentially in the infinity norm and weighted Sobolev
space norms. Numerical results are presented to confirm the theoretical prediction of the
exponential rate of convergence.

Crown Copyright 2009 Published by Elsevier B.V. All rights reserved.

1. Introduction

We consider the linear Volterra integral equations (VIEs) of the second kind, with weakly singular kernels

y(t) = g(t)+
t
0
(t s)K(t, s)y(s)ds, t I, (1.1)

where I = [0, T ], the function g C(I), y(t) is the unknown function, (0, 1) and K C(I I) with K(t, t) 6= 0 for
t I . Several numerical methods have been proposed for (1.1) (see, e.g., [111]).
The numerical treatment of the VIEs (1.1) is not simple, mainly due to the fact that the solutions of (1.1) usually have a

weak singularity at t = 0, even when the inhomogeneous term g(t) is regular. As discussed in [4], the first derivative of
the solution y(t) behaves like y(t) t. In [12], a Jacobi-collocation spectral method is developed for (1.1). To handle the
non-smoothness of the underlying solutions, both function transformation and variable transformation are used to change
the equation into a new Volterra integral equation defined on the standard interval [1, 1], so that the solution of the new
equation possesses better regularity and the Jacobi orthogonal polynomial theory can be applied conveniently. However,
the function transformation (see also [9]) generally makes the resulting equations and approximations more complicated.
We point out that for (1.1) without the singular kernel (i.e., = 0), spectral methods and the corresponding error analysis

Corresponding author.
E-mail addresses: [email protected] (Y. Chen), [email protected] (T. Tang).

0377-0427/$ see front matter Crown Copyright 2009 Published by Elsevier B.V. All rights reserved.
doi:10.1016/j.cam.2009.08.057

Authors personal copy

Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950 939

have been provided recently [13,14]; see also [15,16] for spectral methods to pantograph-type delay differential equations.
In both cases, the underlying solutions are smooth.
In this work, we will consider a special case, namely, the exact solutions of (1.1) are smooth. This case may occur when

the source function g in (1.1) is non-smooth; see, e.g., Theorem 6.1.11 in [4]. In this case, the Jacobi-collocation spectral
method can be applied directly; and the main purpose of this work is to carry out an error analysis for the spectral method.
It is known that most systems of weakly singular VIEs that arise in many application areas are of large dimension. There
have been also some recent developments for solving systems of weakly singular VIEs of high dimensions, see, e.g., [58]. It
is pointed out that although problem (1.1) under consideration is scalar, the proposed methods can be applied to systems
of large dimension in a quite straightforward way. We will demonstrate this for a two-dimensional case in Section 2.
This paper is organized as follows. In Section 2, we introduce the spectral approaches for the Volterra integral equations

with weakly singular kernels. We present some lemmas useful for establishing the convergence results in Section 3. The
convergence analysis is provided in Section 4. Numerical experiments are carried out in Section 5, which will be used to
verify the theoretical results obtained in Section 4.

2. Jacobi-collocation methods

Throughout the paper C will denote a generic positive constant that is independent of N but which will depend on the
length T of the interval I = [0, T ] and on bounds for the given functions f , K which will be defined in (2.5), and the index.
Let ,(x) = (1 x)(1+ x) be a weight function in the usual sense, for , > 1. As defined in [1720], the set of

Jacobi polynomials {J,n (x)}n=0 forms a complete L
2
,

(1, 1)-orthogonal system, where L2
,

(1, 1) is a weighted space
defined by

L2
,

(1, 1) =
{
v : v is measurable and v, <},equipped with the normv, =( 11|v(x)|2,(x)dx) 12,and the inner product(u, v), = 11u(x)v(x),(x)dx, u, v L2,(1, 1).For a given positive integer N , we denote the collocation points by {xi}Ni=0, which is the set of (N + 1) JacobiGauss, orJacobiGaussRadau, or JacobiGaussLobatto points, and by {wi}Ni=0 the corresponding weights. Let PN denote the spaceof all polynomials of degree not exceeding N . For any v C[1, 1], we can define the Lagrange interpolating polynomialI,N v PN , satisfyingI,N v(xi) = v(xi), 0 i N, (2.1)see, e.g., [17,18,20]. The Lagrange interpolating polynomial can be written in the formI,N v(x) =Ni=0v(xi)Fi(x),where Fi(x) is the Lagrange interpolation basis function associated with {xi}Ni=0.In this paper, we deal with the special case that = , = 0, ,0(x) = (1 x).2.1. Numerical scheme in one dimensionFor the sake of applying the theory of orthogonal polynomials, we use the change of variablet =12T (1+ x), x =2tT 1,to rewrite the weakly singular VIEs problem (1.1) as followsu(x) = f (x)+ T (1+x)/20(T2(1+ x) s)K(T2(1+ x), s)y(s)ds, (2.2)Author’s personal copy940 Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950where x [1, 1], andu(x) = y(T2(1+ x)), f (x) = g(T2(1+ x)). (2.3)Furthermore, to transfer the integral interval [0, T (1 + x)/2] to the interval [1, x], we make a linear transformation:s = T (1+ )/2, [1, x]. Then, Eq. (2.2) becomesu(x) = f (x)+ x1(x )K(x, )u( )d , x [1, 1], (2.4)whereK(x, ) =(T2)1K(T2(1+ x),T2(1+ )). (2.5)Firstly, Eq. (2.4) holds at the collocation points {xi}Ni=0 on [1, 1]:u(xi) = f (xi)+ xi1(xi )K(xi, )u( )d , 0 i N. (2.6)In order to obtain high order accuracy for the VIEs problem, the main difficulty is to compute the integral term in (2.6). Inparticular, for small values of xi, there is little information available for u( ). To overcome this difficulty, we first transferthe integral interval [1, xi] to a fixed interval [1, 1] xi1(xi )K(xi, )u( )d =(1+ xi2)1 11(1 )K(xi, i())u(i())d, (2.7)by using the following variable change = i() =1+ xi2 +xi 12, [1, 1]. (2.8)Next, using a (N + 1)-point Gauss quadrature formula relative to the Jacobi weight {wi}Ni=0, the integration term in (2.6) canbe approximated by 11(1 )K(xi, i())u(i())d Nk=0K(xi, i(k))u(i(k))wk, (2.9)where the set {k}Nk=0 coincides with the collocation points {xi}Ni=0 on [1, 1]. We use ui, 0 i N to approximate thefunction value u(xi), 0 i N , and useuN(x) =Nj=0ujFj(x) (2.10)to approximate the function u(x), namely, u(xi) ui, u(x) uN(x), andu(i(k)) Nj=0ujFj(i(k)). (2.11)Then, the Jacobi-collocation method is to seek uN(x) such that {ui}Ni=0 satisfies the following collocation equations:ui = f (xi)+(1+ xi2)1 Nj=0uj(Nk=0K(xi, i(k))Fj(i(k))wk), 0 i N. (2.12)It follows from (2.3) that the exact solution of the VIEs problem (1.1) can be written asy(t) = y(T2(1+ x))= u(x), t [0, T ] and x [1, 1]. (2.13)Thus, we can defineyN(t) = yN(T2(1+ x))= uN(x), t [0, T ] and x [1, 1], (2.14)as the approximated solution of the VIEs problem (1.1). It is obvious to see that(y yN)(t) = (u uN)(x) := e(x), t [0, T ] and x [1, 1]. (2.15)Author’s personal copyY. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950 9412.2. Two-dimensional extensionProblem (1.1) is considered to be scalar; however, many applications involve systems of weakly singular VIEs with highdimensions. The spectral methods proposed in the last subsection are generalizable to large systems of VIEs and to higherdimensions. To demonstrate this, we briefly outline how to solve the second-kind VIEs in two dimension:y(s, t) = g(s, t)+ s0 t0(s )(t )K(s, t, , )y( , )dd , (2.16)where (s, t) [0, T ]2. By using some linear transformations as in Section 2.1, Eq. (2.16) becomesu(x, y) = f (x, y)+ x1 y1(x )(y ) K(x, y, , )u( , )dd, (2.17)where (x, y) [1, 1]2 andK(x, y, , ) =(T2)(2)K(T2(1+ x),T2(1+ y),T2(1+ ),T2(1+ )).For the weight function ,0(x), we denote the collocation points by {xi}Ni=0, which is the set of (N + 1) JacobiGauss,or JacobiGaussRadau, or JacobiGaussLobatto points, and by {w(1)i }Ni=0 the corresponding weights. Similarly, for theweight function ,0(y), we denote the collocation points by {yj}Nj=0, which is the set of (N + 1) JacobiGauss, orJacobiGaussRadau, or JacobiGaussLobatto points, and by {w(2)j }Nj=0 the corresponding weights. Assume that Eq. (2.17)holds at the Jacobi-collocation point-pairs (xi, yj). Using the linear transformation and tricks used in one-dimensional caseyieldsui,j = f (xi, yj)+(1+ xi2)1 (1+ yj2)1 Nk=0Nl=0uk,lak,l, (2.18)whereak,l =Nm=0Nn=0K(xi, yj, i(m), j(n))Fk(i(m))Fl(j(n))w(1)m w(2)n ,i(m) =1+ xi2m +xi 12, j(n) =1+ yj2n +yj 12,for any 0 i, j N .3. Some useful lemmasWe first introduce some weighted Hilbert spaces. For simplicity, denote xv(x) = (/x)v(x), etc. For non-negativeintegerm, defineHm,(1, 1) :={v : kx v L2,(1, 1), 0 k m},with the semi-norm and the norm as|v|m,, = mx v, , vm,, =(mk=0|v|2k,,)1/2,respectively. Let (x) = 12 ,12 (x) denote the Chebyshev weight function. In bounding some approximation error ofChebyshev polynomials, only some of the L2-norms appearing on the right-hand side of above norm enter into play. Thus,it is convenient to introduce the semi-norms|v|Hm;N (1,1)=(mk=min(m,N+1)kx v2L2(1,1)) 12.For bounding some approximation error of Jacobi polynomials, we need the following non-uniformly weighted Sobolevspaces:Hm, ,(1, 1) :={v : kx v L2+k,+k(1, 1), 0 k m},Author’s personal copy942 Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950equipped with the inner product and the norm as(u, v)m,, , =mk=0(kxu, kx v)+k,+k , vm,, , =(v, v)m,, ,.Furthermore, we introduce the orthogonal projection N,, : L2,(1, 1) PN , which is a mapping such that for anyv L2,(1, 1),(v N,, v, ), = 0, PN .It follows from Theorem 1.8.1 in [20] and (3.18) in [18] thatLemma 3.1. Let , > 1. Then for any function v Hm
, ,

(1, 1) and any non-negative integer m, we have

k
x (v N,, v)+k,+k CN

km

m
x v+m,+m , 0 k m. (3.1)

In particular,

v N,, v, CN
1
|v|1,+1,+1 . (3.2)

Applying Theorem 1.8.4 in [20] Theorem 4.3, 4.7, and 4.10 in [21], we have the following optimal error estimate for the
interpolation polynomials based on the JacobiGauss points, JacobiGaussRadau points, and GaussLobatto points.

Lemma 3.2. For any function v satisfying xv Hm, ,(1, 1), we have, for 0 k m,

k
x (v I

,

N v)+k,+k CN
km

m
x v+m,+m , (3.3)

|v I,N v|1,, C (N(N + + ))
1m/2

m
x v+m,+m . (3.4)

Let us define a discrete inner product. For any u, v C[1, 1], define

(u, v)N =
N
j=0

u(xj)v(xj)wj. (3.5)

Due to (5.3.4) in [17], and Lemmas 3.1 and 3.2, we can have the integration error estimates from the JacobiGauss polyno-
mials quadrature.

Lemma 3.3. Let v be any continuous function on [1, 1] and be any polynomial of PN . For the JacobiGauss and JacobiGauss
Radau integration, we have

|(v, ), (v, )N | v I
,

N v,,

CNmmx v+m,+m, . (3.6)

For the JacobiGaussLobatto integration, we have

|(v, ), (v, )N | C
(
v N1,, v, + v I

,

N v,

)
,

CNmmx v+m,+m, . (3.7)

From [22], we have the following result on the Lebesgue constant for the Lagrange interpolation polynomials associated
with the zeros of the Jacobi polynomials.

Lemma 3.4. Assume = , = 0 and assume that Fj(x) is the corresponding Nth Lagrange interpolation polynomials
associated with the Gauss, or GaussRadau, or GaussLobatto points of the Jacobi polynomials. Then

I,0N := max
x(1,1)

N
j=0

|Fj(x)| = O(

N). (3.8)

The following generalization of Gronwalls Lemma for singular kernels, whose proof can be found, e.g. in [23] (Lemma
7.1.1), will be essential for establishing our main results.

Authors personal copy

Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950 943

Lemma 3.5. Suppose L 0, 0 < < 1 and v(t) is a non-negative, locally integrable function defined on [0, T ] satisfyingu(t) v(t)+ L t0(t s)u(s)ds. (3.9)Then there exists a constant C = C() such thatu(t) v(t)+ CL t0(t s)v(s)ds for 0 t < T . (3.10)From now on, for r 0 and [0, 1], C r,([1, 1])will denote the space of functions whose rth derivatives are Hldercontinuous with exponent , endowed with the usual normvr, = max0krmaxx[1,1]|kx v(x)| + max0krsupx,y[1,1]x6=y|kx v(x) kx v(y)||x y|.When = 0, C r,0([1, 1]) denotes the space of functions with r continuous derivatives on [1, 1], which is also commonlydenoted by C r([1, 1]), and with norm r .We shall make use of a result of Ragozin [24,25] (see also [26]), which states that, for non-negative integers r and (0, 1), there exists a constant Cr, > 0 such that for any function v C r,([1, 1]), there exists a polynomial function
TNv PN such that

v TNv Cr,N
(r+)

vr, . (3.11)

Actually, as stated in [24,25], TN is a linear operator from C r,([1, 1]) into PN .
We further define a linear, weakly singular integral operatorM:

Mv =

x
1
(x )K(x, )v( )d . (3.12)

Below we will show thatM is compact as an operator from C([1, 1]) to C0,([1, 1]) provided that the index satisfies
0 < < 1. A similar result can be found in Theorem 3.4 of [27]. The proof of the following lemma can be found in [12].Lemma 3.6. Let (0, 1) andM be defined by (3.12) under the assumption that 0 < < 1 . Then, for any function v C([1, 1]), there exists a positive constant C, which is dependent of K0, , such that|Mv(x)Mv(x)||x x| C maxx[1,1]|v(x)|, (3.13)for any x, x [1, 1] and x 6= x. This implies thatMv0, Cv, (3.14)where is the standard norm in C([1, 1]).4. Convergence analysis4.1. Error estimate in LTheorem 4.1. Let u be the exact solution to the Volterra integral equation (2.4), which is assumed to be sufficiently smooth. Letthe approximated solution uN be obtained by using the spectral collocation scheme (2.12) togetherwith a polynomial interpolation(2.10). If associated with the weakly singular kernel satisfies 0 < < 12 and u Hm (1, 1) Hm,0,(1, 1), thenu uN CN1m|u|Hm;N (1,1)+ CN1/2m K u, (4.1)for N sufficiently large, where i() is defined by (2.8) andK := max0iNm K(xi, i())m,m .Proof. First, we use the weighted inner product to rewrite (2.6) asu(xi) = f (xi)+(1+ xi2)1 (K(xi, i()), u(i())),0, 0 i N. (4.2)Author’s personal copy944 Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950By using the discrete inner product (3.5), we set(K(xi, i()), (i()))N=Nk=0K(xi, i(k))(i(k))wk.Then, the numerical scheme (2.12) can be written asui = f (xi)+(1+ xi2)1 (K(xi, i()), uN(i()))N, 0 i N, (4.3)where uN is defined by (2.10). Subtracting (4.3) from (4.2) gives the error equation:u(xi) ui =(1+ xi2)1 (K(xi, i()), e(i())),0+(1+ xi2)1Ii,2= xi1(xi )K(xi, )e( )d +(1+ xi2)1Ii,2, (4.4)for 0 i N , where e(x) = u(x) uN(x) is the error function,Ii,2 =(K(xi, i()), uN(i())),0(K(xi, i()), uN(i()))N,and the integral transformation (2.7) was used here. Using the integration error estimates from JacobiGauss polynomialsquadrature in Lemma 3.3, we have(1+ xi2)1Ii,2 CNmm K(xi, i())m,muN(i()),0 . (4.5)Multiplying Fi(x) on both sides of the error equation (4.4) and summing up from i = 0 to i = N yieldI,0N u uN= I,0N( x1(x )K(x, )e( )d)+Ni=0(1+ xi2)1Ii,2Fi(x). (4.6)Consequently,e(x) = x1(x )K(x, )e( )d + I1 + I2 + I3, (4.7)whereI1 = u I,0N u, I2 =Ni=0(1+ xi2)1Ii,2Fi(x), (4.8a)I3 = I,0N( x1(x )K(x, )e( )d) x1(x )K(x, )e( )d . (4.8b)It follows from the Gronwall inequality (Lemma 3.5)e C(I1 + I2 + I3). (4.9)Let IcNu PN denote the interpolant of u at any of the three families of ChebyshevGauss points. From (5.5.28) in [17], theinterpolation error estimate in the maximum norm is given byu IcNu CN1/2m|u|Hm;N (1,1). (4.10)Note thatI,0N p(x) = p(x), i.e., (I,0N I)p(x) = 0, p(x) PN . (4.11)By using (4.11), Lemma 3.4 and (4.10), we obtained thatI1 = u I,0N u= u IcNu+ I,0N (IcNu) I,0N u u IcNu + I,0N (IcNu u)Author’s personal copyY. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950 945(1+ I,0N )u IcNu(1+N) N1/2m|u|Hm;N (1,1) CN1m|u|Hm;N (1,1). (4.12)Next, using the estimate (4.5) and Lemma 3.4, we haveI2 CN1/2m max0iNm K(xi, i())m,m max0iNuN(i()),0 CN1/2m max0iNm K(xi, i())m,m (e + u)13e + CN1/2m max0iNm K(xi, i())m,m u, (4.13)provided that N is sufficiently large. We now estimate the third term I3. It follows from (3.11) and (4.11), and Lemma 3.4thatI3 = (I,0N I)Me = (I,0N I)(Me TNMe) (1+ I,0N ) Me TNMe CN Me TNMe CN12Me0, CN12e, (4.14)where in the last step we have used Lemma 3.6 under the assumption 0 < < 1 . It is clear that if > 12 (which is
equivalent to < 12 ), thenI3 13e, (4.15)provided that N is sufficiently large. Combining (4.9), (4.12), (4.13) and (4.15) gives the desired estimate (4.1). 4.2. Error estimate in weighted L2 normTo prove the error estimate in weighted L2 norm, we need the generalized Hardys inequality with weights (see, e.g.,[2830]).Lemma 4.1. For all measurable function f 0, the following generalized Hardys inequality( ba|(Tf )(x)|qu(x)dx)1/q C( ba|f (x)|pv(x)dx)1/pholds if and only ifsupa 0. Then, we can see that

I3,0 I3,1,0 + I3,2,0 , (4.23)

where

I3,1 =
(
I,0N I

) x
1

(x )K(x, )e( )d ,

I3,2 =
(
I,0N I

) x
x
(x )K(x, )e( )d .

Here I denotes the identical operator. It follows from Lemma 3.2 that

I3,1,0 CN
1
x

( x
1

(x )K(x, )e( )d
)

1,1

= CN1
I(1)3,1 + I(2)3,2

1,1
, (4.24)

Authors personal copy

Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950 947

where I(1)3,1 :=
K(x, x )e(x ) and

I(2)3,2 :=
x
1

[
(x )1K(x, )+ (x )Kx(x, )

]
e( )d . (4.25)

Extending e 0 for x 0, we can easily obtain that

I(1)3,11,1

K(, )0,e( )1,1 C

K(, )0,e,0 . (4.26)

It follows from (4.25) that

I(2)3,11,1 C

K(, )0,e,0 + CK(, )1,e,0 . (4.27)

Combining the above two estimates leads to

I3,1,0 C(N
1

+ N1)e,0 . (4.28)

We choose N2, so that N1 N(12) which can be sufficiently small due to the assumption that 0 < < 12 .Therefore, for sufficiently large N , we haveI3,1,0 14e,0 . (4.29)On the other hand, we haveI3,2,0 I(1)3,2,0 + I(2)3,2,0 , (4.30)whereI(1)3,2 := I,0N xx(x )K(x, )e( )d ,I(2)3,2 := xx(x )K(x, )e( )d .It follows from (2.25) of [18] thatI(1)3,22,0=Ni=0 xixi(xi )K(xi, )e( )d2 wi, (4.31)where {xi}Ni=0 is the set of (N + 1) JacobiGauss, or JacobiGaussRadau, or JacobiGaussLobatto points, and {wi}Ni=0 arethe corresponding weights. From (2.13) of [18] and noting that < 1/2, we haveI(1)3,22,0 CN1Ni=0 xixi(xi )K(xi, )e( )d2 (1 xi)+1/2(1+ xi)1/2 CN1Ni=0 xixi(xi )K(xi, )e( )d2 .From Cauchy inequality and using the fact that K is bounded, we obtainI(1)3,22,0 CN1Ni=0 xixi(1 )(xi )2d xixi(1 )e2( )d CN1e2,0, (4.32)where we have used the assumptions N2 1 and 1 2 > 0. Again, using Cauchy inequality and the boundedness
of K gives

|I(2)3,2|
2
C

x
x
(1 )(x )2d

x
x
(1 )e2( )d C12e2

,0
.

Consequently,

I(2)3,2
2
,0

CN2(12)e2
,0

, (4.33)

Authors personal copy

948 Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950

where we also used N2 and < 1/2. It follows from (4.32) and (4.33) thatI3,2,0 I3,2,1,0 + I3,2,2,0 14e,0 , (4.34)provided that N is sufficiently large. Hence, by (4.23), (4.29) and (4.34), we haveI3,0 12e,0 , (4.35)for sufficiently large N . The desired estimate (4.17) is obtained by combining (4.18), (4.19), (4.22) and (4.35). 5. Algorithm implementation and numerical experimentsDenoting UN = [u0, u1, . . . , uN ]T and FN = [f (x0), f (x1), . . . , f (xN)]T, we can obtain an equation of the matrix form:UN = FN + AUN , (5.1)where the entries of the matrix A is given byaij =(1+ xi2)1 Nk=0K(xi, i(k))Fj(i(k))wk.Here, we simply introduce the computation of GaussJacobi quadrature rule nodes and weights (see the detailedalgorithm and download related codes in [32]). The GaussJacobi quadrature formula is used to numerically calculate theintegral 11(1 x)(1+ x) f (x)dx, f (x) [1, 1], , > 1,

by using the formula 1
1
(1 x)(1+ x) f (x)dx

N
i=0

wif (xi).

With the help of a change in the variables (which changes both weights wi and nodes xi), we can get onto the arbitrary
interval [a, b].

Example 5.1. We consider the following the linear Volterra integral equations of second kind with weakly singular kernels

y(t) = b(t)
t
0
(t s)y(s)ds, 0 t T , (5.2)

with

b(t) = tn+ + tn+1+B(n+ 1+ , 1 ),

where 0 1, and B(, ) is the Beta function defined by

B(x, y) =
1
0
tx1(1 t)y1dt.

This problem has an unique solution: y(t) = tn+ . Obviously, y(t) belongs to Hn+1 (I). It follows from the theoretical
results obtained in this work, the numerical errors will decay with a rate ofO(Nn1). The weighted function is chosen as
= (1 x) with = 0.35. The solution interval is t [0, 6], and the exact solution is y(t) = t3.6, indicating that n = 3
and = 0.6. In Fig. 1, numerical errors are plotted for 2 N 16 in both L and L2-norms. We also present in Table 1
the corresponding numerical errors. As expected, the errors decay algebraically as the exact solution for this example is not
sufficiently smooth.

Example 5.2. Consider the following the nonlinear Volterra integral equations of second kind with weakly singular kernels

y(t) = b(t)+
t
0
(t s)1/3y2(s)ds, 0 t T , (5.3)

where

b(t) = (t + 2)3/2
243
440
t11/3

81
20
t8/3

54
5
t5/3 12t2/3.

Authors personal copy

Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938950 949

1 2 3 4 5

y
(t

)
=

tn
+

Numerical Solution

ExactSolution

101

10
-5

10
-4

10
-3

10
-2

10
-1

10
0

10
1

2 N 16

0

100

200

300

400

500

600

700

0 6

0 t 6

Fig. 1. Left: Numerical and exact solution of y(t) = t3+0.6 with n = 3 and = 0.6. Right: The L error and L2 error versus N .

Table 1
Example 5.1: L error and L2 errors, for the solution interval t [0, 6].

N 2 4 6 8

L error 6.7887e+01 2.4594e01 1.3307e02 1.9500e03
L2 error 2.3692e+01 6.5956e02 1.7042e03 1.4331e04

N 10 12 14 16

L error 4.4826e04 1.3478e04 4.8583e05 1.9980e05
L2 error 2.2145e05 7.8788e06 6.6148e06 6.5174e06

Table 2
Example 5.2: L error and L2 errors, for the solution interval t [0, 10].

N 6 8 10 12

L error 1.0571e03 1.1271e04 1.3630e05 1.7825e06
L2 error 3.8912e04 3.6562e05 3.9691e06 4.7263e07

N 14 16 18 20

L error 2.4594e07 3.5287e08 5.2114e09 7.9513e010
L2 error 6.0029e08 7.9972e09 1.1046e09 1.5821e010

This example has a smooth solution y(t) = (2 + t)3/2. As a result, we can expect an exponential rate of convergence.
In Fig. 2, numerical errors are plotted for 2 N 24 in both L- and L2-norms. The weighted function is chosen as
= (1 x) with = 1/3. The solution interval is t [0, 10]. We also present in Table 2 the corresponding numerical
errors. As expected, the errors decay exponentially which confirmed our theoretical predictions.

Acknowledgments

The authors are grateful to Mr. Xiang Xu of Fudan University for his assistance in the numerical work. The first author
is supported by Guangdo