[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 Elsevier’s archiving and manuscript policies are
encouraged to visit:

http://www.elsevier.com/copyright

http://www.elsevier.com/copyright

Author’s personal copy

Journal of Computational and Applied Mathematics 233 (2009) 938–950

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., [1–11]).
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

Author’s personal copy

Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938–950 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., [5–8]. 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 [17–20], 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 norm‖v‖ωα,β =(∫ 1−1|v(x)|2ωα,β(x)dx) 12,and the inner product(u, v)ωα,β =∫ 1−1u(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) Jacobi–Gauss, orJacobi–Gauss–Radau, or Jacobi–Gauss–Lobatto 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) =N∑i=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) 938–950where 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)+∫ x−1(x− τ)−µK̃(x, τ )u(τ )dτ , x ∈ [−1, 1], (2.4)whereK̃(x, τ ) =(T2)1−µK(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)+∫ xi−1(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]∫ xi−1(xi − τ)−µK̃(xi, τ )u(τ )dτ =(1+ xi2)1−µ ∫ 1−1(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∫ 1−1(1− θ)−µK̃(xi, τi(θ))u(τi(θ))dθ ∼N∑k=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) =N∑j=0ujFj(x) (2.10)to approximate the function u(x), namely, u(xi) ≈ ui, u(x) ≈ uN(x), andu(τi(θk)) ≈N∑j=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−µ N∑j=0uj(N∑k=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) 938–950 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(σ , τ )dσdτ , (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)+∫ x−1∫ y−1(x− ξ)−α(y− η)−β K̃(x, y, ξ , η)u(ξ , η)dξdη, (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) Jacobi–Gauss,or Jacobi–Gauss–Radau, or Jacobi–Gauss–Lobatto 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) Jacobi–Gauss, orJacobi–Gauss–Radau, or Jacobi–Gauss–Lobatto 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−β N∑k=0N∑l=0uk,lak,l, (2.18)whereak,l =N∑m=0N∑n=0K̃(xi, yj, ξi(θm), ηj(θn))Fk(ξi(θm))Fl(ηj(θn))w(1)m w(2)n ,ξi(θm) =1+ xi2θm +xi − 12, ηj(θn) =1+ yj2θn +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‖ωα,β , ‖v‖m,ωα,β =(m∑k=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)=(m∑k=min(m,N+1)‖∂kx v‖2L2ω(−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) 938–950equipped with the inner product and the norm as(u, v)m,ωα,β ,∗ =m∑k=0(∂kxu, ∂kx v)ωα+k,β+k , ‖v‖m,ωα,β ,∗ =√(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

k−m
‖∂
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 Jacobi–Gauss points, Jacobi–Gauss–Radau points, and Gauss–Lobatto 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
k−m
‖∂
m
x v‖ωα+m,β+m , (3.3)

|v − Iα,βN v|1,ωα,β ≤ C (N(N + α + β))
1−m/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 Jacobi–Gauss polyno-
mials quadrature.

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

|(v, φ)ωα,β − (v, φ)N | ≤ ‖v − I
α,β

N v‖ωα,β‖φ‖ωα,β

≤ CN−m‖∂mx v‖ωα+m,β+m‖φ‖ωα,β . (3.6)

For the Jacobi–Gauss–Lobatto integration, we have

|(v, φ)ωα,β − (v, φ)N | ≤ C
(
‖v − πN−1,ωα,β v‖ωα,β + ‖v − I

α,β

N v‖ωα,β

)
‖φ‖ωα,β

≤ CN−m‖∂mx 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 Gauss–Radau, or Gauss–Lobatto 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 Gronwall’s 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.

Author’s personal copy

Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938–950 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 Höldercontinuous with exponent κ , endowed with the usual norm‖v‖r,κ = max0≤k≤rmaxx∈[−1,1]|∂kx v(x)| + max0≤k≤rsupx,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+κ)

‖v‖r,κ . (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 ‖K̃‖0,κ , 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 that‖Mv‖0,κ ≤ C‖v‖∞, (3.14)where ‖ · ‖∞ is the standard norm in C([−1, 1]).4. Convergence analysis4.1. Error estimate in L∞Theorem 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), then‖u− uN‖∞ ≤ CN1−m|u|Hm;Nω (−1,1)+ CN1/2−m · K ∗‖u‖∞, (4.1)for N sufficiently large, where τi(θ) is defined by (2.8) andK ∗ := max0≤i≤N‖∂mθ 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) 938–950By using the discrete inner product (3.5), we set(K̃(xi, τi(·)), φ(τi(·)))N=N∑k=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)1−µIi,2=∫ xi−1(xi − τ)−µK̃(xi, τ )e(τ )dτ +(1+ xi2)1−µIi,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 Jacobi–Gauss polynomialsquadrature in Lemma 3.3, we have∣∣∣∣∣(1+ xi2)1−µIi,2∣∣∣∣∣ ≤ CN−m‖∂mθ K̃(xi, τi(·))‖ωm−µ,m‖uN(τ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(∫ x−1(x− τ)−µK̃(x, τ )e(τ )dτ)+N∑i=0(1+ xi2)1−µIi,2Fi(x). (4.6)Consequently,e(x) =∫ x−1(x− τ)−µK̃(x, τ )e(τ )dτ + I1 + I2 + I3, (4.7)whereI1 = u− I−µ,0N u, I2 =N∑i=0(1+ xi2)1−µIi,2Fi(x), (4.8a)I3 = I−µ,0N(∫ x−1(x− τ)−µK̃(x, τ )e(τ )dτ)−∫ x−1(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 Chebyshev–Gauss points. From (5.5.28) in [17], theinterpolation error estimate in the maximum norm is given by‖u− IcNu‖∞ ≤ CN1/2−m|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 that‖I1‖∞ = ‖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) 938–950 945≤(1+ ‖I−µ,0N ‖∞)‖u− IcNu‖∞≤(1+√N)· N1/2−m|u|Hm;Nω (−1,1)≤ CN1−m|u|Hm;Nω (−1,1). (4.12)Next, using the estimate (4.5) and Lemma 3.4, we have‖I2‖∞ ≤ CN1/2−m max0≤i≤N‖∂mθ K̃(xi, τi(·))‖ωm−µ,m · max0≤i≤N‖uN(τi(·))‖ω−µ,0≤ CN1/2−m max0≤i≤N‖∂mθ K̃(xi, τi(·))‖ωm−µ,m · (‖e‖∞ + ‖u‖∞)≤13‖e‖∞ + CN1/2−m max0≤i≤N‖∂mθ 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.4that‖I3‖∞ = ‖(I−µ,0N − I)Me‖∞ = ‖(I−µ,0N − I)(Me− TNMe)‖∞≤ (1+ ‖I−µ,0N ‖∞) · ‖Me− TNMe‖∞≤ C√N · ‖Me− TNMe‖∞ ≤ CN12−κ‖Me‖0,κ≤ CN12−κ‖e‖∞, (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 ), then‖I3‖∞ ≤13‖e‖∞, (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 Hardy’s inequality with weights (see, e.g.,[28–30]).Lemma 4.1. For all measurable function f ≥ 0, the following generalized Hardy’s 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

= CN−1
∥∥∥I(1)3,1 + I(2)3,2∥∥∥

ω1−µ,1
, (4.24)

Author’s personal copy

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

where I(1)3,1 := δ
−µK̃(x, x− δ)e(x− δ) and

I(2)3,2 :=
∫ x−δ
−1

[
−µ(x− τ)−µ−1K̃(x, τ )+ (x− τ)−µK̃x(x, τ )

]
e(τ )dτ . (4.25)

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

‖I(1)3,1‖ω1−µ,1 ≤ δ
−µ
‖K̃(·, ·)‖0,∞‖e(· − δ)‖ω1−µ,1 ≤ Cδ

−µ
‖K̃(·, ·)‖0,∞‖e‖ω−µ,0 . (4.26)

It follows from (4.25) that

‖I(2)3,1‖ω1−µ,1 ≤ Cδ
−µ
‖K̃(·, ·)‖0,∞‖e‖ω−µ,0 + C‖K̃(·, ·)‖1,∞‖e‖ω−µ,0 . (4.27)

Combining the above two estimates leads to

‖I3,1‖ω−µ,0 ≤ C(N
−1
δ
−µ
+ N−1)‖e‖ω−µ,0 . (4.28)

We choose δ ≤ N−2, so that N−1δ−µ ≤ N−(1−2µ) which can be sufficiently small due to the assumption that 0 < µ < 12 .Therefore, for sufficiently large N , we have‖I3,1‖ω−µ,0 ≤14‖e‖ω−µ,0 . (4.29)On the other hand, we have‖I3,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] that‖I(1)3,2‖2ω−µ,0=N∑i=0∣∣∣∣∫ xixi−δ(xi − τ)−µK̃(xi, τ )e(τ )dτ∣∣∣∣2 · wi, (4.31)where {xi}Ni=0 is the set of (N + 1) Jacobi–Gauss, or Jacobi–Gauss–Radau, or Jacobi–Gauss–Lobatto points, and {wi}Ni=0 arethe corresponding weights. From (2.13) of [18] and noting that µ < 1/2, we have‖I(1)3,2‖2ω−µ,0≤ CN−1N∑i=0∣∣∣∣∫ xixi−δ(xi − τ)−µK̃(xi, τ )e(τ )dτ∣∣∣∣2 (1− xi)−µ+1/2(1+ xi)1/2≤ CN−1N∑i=0∣∣∣∣∫ xixi−δ(xi − τ)−µK̃(xi, τ )e(τ )dτ∣∣∣∣2 .From Cauchy inequality and using the fact that K̃ is bounded, we obtain‖I(1)3,2‖2ω−µ,0≤ CN−1N∑i=0∫ xixi−δ(1− τ)µ(xi − τ)−2µdτ ·∫ xixi−δ(1− τ)−µe2(τ )dτ≤ CN−1‖e‖2ω−µ,0, (4.32)where we have used the assumptions δ ≤ N−2 ≤ 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− τ)−2µdτ

∫ x
x−δ
(1− τ)−µe2(τ )dτ ≤ Cδ1−2µ‖e‖2

ω−µ,0
.

Consequently,

‖I(2)3,2‖
2
ω−µ,0

≤ CN−2(1−2µ)‖e‖2
ω−µ,0

, (4.33)

Author’s personal copy

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

where we also used δ ≤ N−2 and µ < 1/2. It follows from (4.32) and (4.33) that‖I3,2‖ω−µ,0 ≤ ‖I3,2,1‖ω−µ,0 + ‖I3,2,2‖ω−µ,0 ≤14‖e‖ω−µ,0 , (4.34)provided that N is sufficiently large. Hence, by (4.23), (4.29) and (4.34), we have‖I3‖ω−µ,0 ≤12‖e‖ω−µ,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−µ N∑k=0K̃(xi, τi(θk))Fj(τi(θk))wk.Here, we simply introduce the computation of Gauss–Jacobi quadrature rule nodes and weights (see the detailedalgorithm and download related codes in [32]). The Gauss–Jacobi quadrature formula is used to numerically calculate theintegral∫ 1−1(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
tx−1(1− t)y−1dt.

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(N−n−1). 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.

Author’s personal copy

Y. Chen, T. Tang / Journal of Computational and Applied Mathematics 233 (2009) 938–950 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.4594e−01 1.3307e−02 1.9500e−03
L2ω error 2.3692e+01 6.5956e−02 1.7042e−03 1.4331e−04

N 10 12 14 16

L∞ error 4.4826e−04 1.3478e−04 4.8583e−05 1.9980e−05
L2ω error 2.2145e−05 7.8788e−06 6.6148e−06 6.5174e−06

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.0571e−03 1.1271e−04 1.3630e−05 1.7825e−06
L2ω error 3.8912e−04 3.6562e−05 3.9691e−06 4.7263e−07

N 14 16 18 20

L∞ error 2.4594e−07 3.5287e−08 5.2114e−09 7.9513e−010
L2ω error 6.0029e−08 7.9972e−09 1.1046e−09 1.5821e−010

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