A M\"untz-Collocation spectral method for weakly singular volterra integral equations

TL;DR

This paper introduces a fractional Jacobi-collocation spectral method for solving weakly singular Volterra integral equations, achieving exponential convergence for smooth solutions after a variable transformation, with demonstrated numerical efficiency.

Contribution

The paper develops a new fractional Jacobi polynomial family and constructs an efficient spectral method with proven exponential convergence for weakly singular VIEs.

Findings

Method achieves exponential convergence for smooth solutions.

Numerical examples confirm high efficiency and accuracy.

New fractional Jacobi polynomials facilitate spectral approximation.

Abstract

In this paper we propose and analyze a fractional Jacobi-collocation spectral method for the second kind Volterra integral equations (VIEs) with weakly singular kernel $(x-s)^{-\mu},0<\mu<1$. First we develop a family of fractional Jacobi polynomials, along with basic approximation results for some weighted projection and interpolation operators defined in suitable weighted Sobolev spaces. Then we construct an efficient fractional Jacobi-collocation spectral method for the VIEs using the zeros of the new developed fractional Jacobi polynomial. A detailed convergence analysis is carried out to derive error estimates of the numerical solution in both $L^{\infty}$- and weighted $L^{2}$-norms. The main novelty of the paper is that the proposed method is highly efficient for typical solutions that VIEs usually possess. Precisely, it is proved that the exponential convergence rate can be achieved for solutions which are smooth after the variable change $x\rightarrow x^{1/\lambda}$ for a suitable real number $\lambda$. Finally a series of numerical examples are presented to demonstrate the efficiency of the method.

Equations

$$u(x)=g(x)+{\int }_0^x(x-s{)}^{-\mu }K(x,s)u(s)ds,0<\mu <1,x\in I:=[0,1],$$

$$P^{\lambda}_{n}(\mathbb{R}^{+}):=span\big{\{}1,x^{\lambda},x^{2\lambda},\cdots,x^{n\lambda}\big{\}},$$

$$p_n^{\lambda }(x):=k_n{x}^{n\lambda }+k_{n-1}{x}^{(n-1)\lambda }+\cdots +k_1{x}^{\lambda }+k_0,k_n\ne 0,x\in {\mathbb{R}}^{+},$$

$$(p_n^{\lambda },p_m^{\lambda }{)}_{\omega }={\int }_0^1{p}_n^{\lambda }(x)p_m^{\lambda }(x)\omega (x)dx={\gamma }_n{\delta }_{m,n},$$

$$P^{\lambda}_{n}(I):=span\big{\{}p^{\lambda}_{0},p^{\lambda}_{1},\cdots,p^{\lambda}_{n}\big{\}}.$$

$${\bar{p}}_0^{\lambda }=1,{\bar{p}}_1^{\lambda }=x^{\lambda }-{\alpha }_0,{\bar{p}}_{n+1}^{\lambda }=(x^{\lambda }-{\alpha }_n){\bar{p}}_n^{\lambda }-{\beta }_n{\bar{p}}_{n-1}^{\lambda },n\ge 1,$$

$$\alpha_{n}=\displaystyle\frac{(x^{\lambda}\bar{p}^{\lambda}_{n},\bar{p}^{\lambda}_{n})_{\omega}}{\|\bar{p}^{\lambda}_{n}\|^{2}_{{\color[rgb]{0,0,0}0,}\omega}},\ \ n\geq 0;\ \ \ \ \beta_{n}=\displaystyle\frac{\|\bar{p}_{n}\|^{2}_{{\color[rgb]{0,0,0}0,}\omega}}{\|\bar{p}^{\lambda}_{n-1}\|^{2}_{{\color[rgb]{0,0,0}0,}\omega}},\ \ n\geq 1.$$

$$J_n^{\alpha ,\beta ,\lambda }(x)=J_n^{\alpha ,\beta }(2x^{\lambda }-1),\forall x\in I,$$

$$\displaystyle J^{\alpha,\beta}_{n}(x)=\displaystyle\frac{\Gamma(n+\alpha+1)}{n!\Gamma(n+\alpha+\beta+1)}\sum_{k=0}^{n}\binom{n}{k}\frac{\Gamma(n+k+\alpha+\beta+1)}{\Gamma(k+\alpha+1)}\big{(}\frac{x-1}{2}\big{)}^{k}.$$

$$J_n^{\alpha ,\beta ,\lambda }(x)=\frac{\Gamma (n+\alpha +1)}{n!\Gamma (n+\alpha +\beta +1)}\sum _{k=0}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)\frac{\Gamma (n+k+\alpha +\beta +1)}{\Gamma (k+\alpha +1)}(x^{\lambda }-1{)}^k.$$

$$D_{\lambda }^{1}v(x)=\frac{d}{d{x}^{\lambda }}v(x)=\frac{x^{1-\lambda }}{\lambda }{v}^{\prime }(x),$$

$$D_{\lambda }^{k}v(x)=\stackrel{k}{\overbrace{D_{\lambda }^1\cdot D_{\lambda }^1\cdots D_{\lambda }^1}}v(x),$$

$$\begin{array}[]{r@{}l}{}^{+}D^{1}_{\lambda}v(x)&:=\displaystyle\lim_{\Delta x\rightarrow 0^{+}}\frac{v(x+\Delta x)-v(x)}{(x+\Delta x)^{\lambda}-x^{\lambda}},\\[9.0pt] {}^{-}D^{1}_{\lambda}v(x)&:=\displaystyle\lim_{\Delta x\rightarrow 0^{-}}\frac{v(x+\Delta x)-v(x)}{(x+\Delta x)^{\lambda}-x^{\lambda}}.\\ \end{array}$$

$${\omega }^{\alpha ,\beta ,\lambda }(x):=\lambda (1-x^{\lambda }{)}^{\alpha }{x}^{(\beta +1)\lambda -1}.$$

$${\hat{\omega }}^{\alpha ,\beta ,\lambda }(x):=(1-x^{\lambda }{)}^{\alpha }{x}^{\beta \lambda }={\lambda }^{-1}{x}^{1-\lambda }{\omega }^{\alpha ,\beta ,\lambda }(x).$$

$${\int }_0^1{\omega }^{\alpha ,\beta ,\lambda }(x)J_n^{\alpha ,\beta ,\lambda }(x)J_m^{\alpha ,\beta ,\lambda }(x)dx={\hat{\gamma }}_n^{\alpha ,\beta }{\delta }_{m,n},$$

$${\hat{\gamma }}_n^{\alpha ,\beta }=\frac{\Gamma (n+\alpha +1)\Gamma (n+\beta +1)}{(2n+\alpha +\beta +1)n!\Gamma (n+\alpha +\beta +1)}.$$

$$\begin{array}[]{r@{}l}\mathcal{L}^{\alpha,\beta}_{\lambda}v(x)&=-(\widehat{\omega}^{\alpha,\beta,\lambda}(x))^{-1}D^{1}_{\lambda}\big{\{}(1-x^{\lambda})^{\alpha+1}x^{(\beta+1)\lambda}D^{1}_{\lambda}v(x)\big{\}}.\end{array}$$

$${\mathcal{L}}_{\lambda }^{\alpha ,\beta }{J}_n^{\alpha ,\beta ,\lambda }(x)={\sigma }_n^{\alpha ,\beta }{J}_n^{\alpha ,\beta ,\lambda }(x).$$

$$-(\omega^{\alpha,\beta,\lambda}(x))^{-1}\frac{d}{dx}\Big{\{}\lambda^{-1}(1-x^{\lambda})^{\alpha+1}x^{\beta\lambda+1}\frac{d}{dx}J^{\alpha,\beta,\lambda}_{n}(x)\Big{\}}=\sigma^{\alpha,\beta}_{n}J^{\alpha,\beta,\lambda}_{n}(x),$$

$${\int }_0^1{\omega }^{\alpha +k,\beta +k,\lambda }(x)D_{\lambda }^k{J}_n^{\alpha ,\beta ,\lambda }(x)D_{\lambda }^k{J}_m^{\alpha ,\beta ,\lambda }(x)dx={\hat{h}}_{n,k}^{\alpha ,\beta }{\delta }_{m,n},$$

$${\hat{h}}_{n,k}^{\alpha ,\beta }=\frac{\Gamma (n+\alpha +1)\Gamma (n+\beta +1)\Gamma (n+k+\alpha +\beta +1)}{(2n+\alpha +\beta +1)(n-k)!{\Gamma }^2(n+\alpha +\beta +1)}.$$

$$D_{\lambda }^k{J}_n^{\alpha ,\beta ,\lambda }(x)={\hat{d}}_{n,k}^{\alpha ,\beta }{J}_{n-k}^{\alpha +k,\beta +k,\lambda }(x),$$

$${\hat{d}}_{n,k}^{\alpha ,\beta }=\frac{\Gamma (n+k+\alpha +\beta +1)}{\Gamma (n+\alpha +\beta +1)}.$$

$${\int }_0^1{\omega }^{\alpha +1,\beta +1,\lambda }(x)D_{\lambda }^1{J}_n^{\alpha ,\beta ,\lambda }(x)D_{\lambda }^1{J}_m^{\alpha ,\beta ,\lambda }(x)dx=(J_n^{\alpha ,\beta ,\lambda },{\mathcal{L}}_{\lambda }^{\alpha ,\beta }{J}_m^{\alpha ,\beta ,\lambda }{)}_{{\omega }^{\alpha ,\beta ,\lambda }}={\sigma }_n^{\alpha ,\beta }{\hat{\gamma }}_n^{\alpha ,\beta }{\delta }_{m,n}.$$

$$D_{\lambda }^1{J}_n^{\alpha ,\beta ,\lambda }(x)={\hat{d}}_{n,1}^{\alpha ,\beta }{J}_{n-1}^{\alpha +1,\beta +1,\lambda }(x),$$

$${\hat{d}}_{n,1}^{\alpha ,\beta }=\frac{n{k}_n^{\alpha ,\beta }}{k_{n-1}^{\alpha +1,\beta +1}}=n+\alpha +\beta +1.$$

$$D_{\lambda }^k{J}_n^{\alpha ,\beta ,\lambda }={\hat{d}}_{n,k}^{\alpha ,\beta }{J}_{n-k}^{\alpha +k,\beta +k,\lambda }(x),$$

$${\hat{d}}_{n,k}^{\alpha ,\beta }=\frac{\Gamma (n+k+\alpha +\beta +1)}{\Gamma (n+\alpha +\beta +1)}.$$

$${\int }_0^1{\omega }^{\alpha +k,\beta +k,\lambda }(x)D_{\lambda }^k{J}_n^{\alpha ,\beta ,\lambda }(x)D_{\lambda }^k{J}_m^{\alpha ,\beta ,\lambda }(x)dx={\hat{h}}_{n,k}^{\alpha ,\beta }{\delta }_{m,n},$$

Full text

∗This research is partially supported by NSF of China (Grant numbers 11971408, 51661135011, 11421110001, and 91630204). The third author has received financial support from the French State in the frame of the “Investments for the future” Programme Idex Bordeaux, reference ANR-10-IDEX-03-02 and NSFC/ANR joint program 51661135011/ANR-16- CE40-0026-01.

1School of Mathematics and Statistics, Jiangsu Normal University, 221116 Xuzhou, China.

2School of Mathematical Sciences and Fujian Provincial Key Laboratory of Mathematical Modeling and High Performance Scientific Computing, Xiamen University, 361005 Xiamen, China.

3Bordeaux INP, Laboratoire I2M UMR 5295, 33607 Pessac, France.

4Corresponding author. Email: cjxu@xmu.edu.cn

1. Introduction

Among various methods, spectral methods have proven to be one of the most efficient approaches for solving partial differential equations. The efficiency of spectral methods, however, depends crucially on the regularity of the solution. The convergence of the spectral method is exponential, provided that the solution of the problem is sufficiently smooth. On the other word, the traditional spectral methods will lose high accuracy when facing problems with solutions of limited regularity. Such problems include the second kind Volterra integral equations (VIEs), which we consider in this paper. It has been well known that the solution of the VIEs with weakly singular kernel is of low regularity at the boundaries of the domain.

Volterra integral equations model evolutionary problems with memory, which can be found in a number of disciplines, such as electromagnetic scattering, demography, viscoelastic materials, insurance mathematics, etc. They have been subject of many theoretical and numerical investigations. In this paper we consider the following integral equation:

[TABLE]

where the source function $g(x)\in C(I)$, and $K\in C(I\times I)$ with $K(x,x)\neq 0$ for $x\in I$. Our aim is to design a numerical method for this equation, which will be shown more efficient than the existing methods in capturing the singularities of typical solutions of this kind of equations.

In fact there exists a considerable amount of numerical methods in the literature for the equation under consideration; see, e.g., [26, 31, 28, 18, 7, 8, 9, 11, 12, 13, 14, 20, 21] and the references therein. In the framework of high order methods for example, a Legendre-collocation method, along with rigorous error analysis, was proposed in [30, 29] for VIEs with smooth kernel (i.e., $\mu=0$) and regular source term; There have also been spectral methods [2, 3], used to approximate smooth solutions of delay differential or integral equations with two or more vanishing delays. In [12, 11], Chen et al. proposed and analyzed a spectral Jacobi-collocation method for VIEs under assumption that the underlying solutions are sufficiently smooth. However, as demonstrated in [9], VIEs with weakly singular kernel typically have solutions whose derivatives are unbounded at the left end point, even when the source term $g(x)$ is regular. Partially for this reason, the numerical treatment of the VIEs is not simple. In particular, numerical methods based on approximations of high order polynomials, such as spectral methods or p-version finite element methods, fail to yield high order convergence since increasing the polynomial degree does not improve the approximation accuracy for low regular functions. Nevertheless, efforts have been made to handle this problem in some different circumstances. Let us mention, among others, Chen and Tang [13] for a transformation method to transform 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. Recently, Li, Tang, and Xu [21] have constructed a Chebyshev and Legendre pseudo-spectral Galerkin method for weakly singular VIEs based on a new function transformation to make the solution smooth. Allaei et al. [4] used a transformation of the independent variable for the Volterra integral equation with weakly singular kernel, then constructed a Jacobi collocation method for the transformed equation. It is worthwhile to emphasize that the efficiency of the aforementioned approaches depends on specific form of the exact solution and the assumption of sufficient smoothness of the source term. Very recently, Hou and Xu [17] proposed a general framework using fractional polynomials for some weakly singular integro-differential equations and fractional differential equations. They showed that the convergence rate of the fractional spectral methods is closely related to the regularity of the transformation of the exact solution by the variable change $x\rightarrow x^{1/\lambda}$, where $\lambda$ is a suitable parameter appearing in the fractional approximation space.

This paper aims at designing, developing, and testing a fractional Jacobi spectral method for the weakly singular VIEs, which has the capability to capture the limited regular solution in a more efficient way. The new method will make use of the fractional Jacobi polynomial $J^{\alpha,\beta,\lambda}_{N+1}(x)$, recently introduced in [17, 16] to deal with some singular problems. The advantage of the proposed approach is that the exponential convergence can be guaranteed for solutions, which are smooth after the variable change $x\rightarrow x^{1/\lambda}$ for suitable parameter $\lambda$. It is worth to emphasize that the present method is different from the one proposed in [17]. Firstly, the method in [17] was for the integral equation: $u_{x}=a_{1}u(x)+a_{2}~{}{{}_{0}}I_{x}^{\mu}u(x)+f(x),$ where $a_{1}$ and $a_{2}$ are constants, which were assumed to satisfy some conditions in order to guarantee the well-posedness of the discrete problem; see Theorem 4.1 in [17]. Secondly, the present method is of Collocation type while the method in [17] was constructed under Galerkin or Petrov-Galerkin framework. It is known that the method of Collocation type is easier to implement (usually more difficult to analyze) than the Galerkin type.

Below are the main ingredients of the paper.

$\bullet$ First, one key point of the study is the new $\lambda-$polynomial space constructed as the approximation space. Some new derivatives are defined such that the fractional Jacobi polynomials can inherit some desirable properties from classical Jacobi polynomials. We also set up the relationship between the new defined derivatives and the classical ones. Then we derive the approximating results for some projection and interpolation operators in weighted Sobolev spaces, This is the main content of Section 2.

$\bullet$ In Section 3, the fractional Jacobi spectral-collocation method is constructed for weakly singular VIEs. A rigorous convergence analysis is carried out for the proposed method and some optimal error estimates in $L^{\infty}-$ and $L^{2}_{\omega^{\alpha,\beta,\lambda}}-$norms are derived.

$\bullet$ Some numerical examples are presented in Section 4 to validate the theoretical prediction. Finally we give some concluding remarks in the last section.

2. Preliminaries

In this section, we will define the fractional Jacobi polynomials and study their fundamental properties. Especially we will introduce some weighted projection and interpolation operators and derive optimal error estimates for these operators in different weighted Sobolev spaces. These results play a key role in the convergence analysis for the numerical method that we design for the VIEs later on.

2.1. Some basic properties of fractional Jacobi polynomials

We begin by defining the $\lambda-$polynomial space as follows:

[TABLE]

where $\mathbb{R}^{+}=[0,+\infty),0<\lambda\leq 1.$

A $\lambda-$polynomial of degree $n$ is denoted by

[TABLE]

where \big{\{}k_{i}\big{\}}_{i=0}^{n} are real constants and $k_{n}$ is called the leading coefficient of $p^{\lambda}_{n}.$ Hereafter, we simply denote the degree of $\lambda-$polynomial $p^{\lambda}_{n}(x)$ by $\deg(p^{\lambda}_{n})$.

Let $I:=[0,1]$ and $\omega(x)\in L^{1}(I)$ be a positive weight function. A sequence of $\lambda-$polynomials \big{\{}p^{\lambda}_{n}\big{\}}_{n=0}^{\infty} with $\deg(p^{\lambda}_{n})=n$ is said to be orthogonal in $L^{2}_{\omega}(I)$ if

[TABLE]

where \gamma_{n}=\|p_{n}^{\lambda}\|^{2}_{{\color[rgb]{0,0,0}0},\omega}{\color[rgb]{0,0,0}:=(p^{\lambda}_{n},p^{\lambda}_{n})_{\omega}}, and $\delta_{m,n}$ is the Kronecker delta.

We define the space

[TABLE]

The following two lemmas can be easily proved by following the standard way to prove the existence of the classical orthogonal polynomials; see, e.g., [27, p48-50].

Lemma 2.1**.**

$p^{\lambda}_{n+1}(x)$* is $\omega(x)$-weighted orthogonal to any $\lambda-$polynomial $q(x)\in P^{\lambda}_{n}(I)$.*

Lemma 2.2**.**

For any given positive weight function $\omega\in L^{1}(I)$, there exists a unique sequence of monic orthogonal $\lambda-$polynomials \big{\{}\bar{p}^{\lambda}_{n}\big{\}}_{n=0}^{\infty} with $\deg(\bar{p}^{\lambda}_{n})=n$. This unique sequence can be obtained through the following recurrence relation

[TABLE]

where

[TABLE]

Now we turn to define the fractional Jacobi polynomials.

Definition 2.1**.**

The fractional Jacobi polynomials of degree $n$ are defined by

[TABLE]

where $J^{\alpha,\beta}_{n}(x)$ denotes the Jacobi polynomial of degree n, and $\alpha,\beta>-1,0<\lambda\leq 1$.

When $\lambda=1,$ the polynomials $\{J^{\alpha,\beta,1}_{n}(x)\}_{n=0}^{\infty}$ are called shifted Jacobi polynomials up to a constant, which are orthogonal polynomials with the weight $(1-x)^{\alpha}x^{\beta}.$

It has been well known that the classical Jacobi polynomial $J^{\alpha,\beta}_{n}(x)$ has the following representation

[TABLE]

Consequently, we have

[TABLE]

In order to inherit important properties from the classical Jacobi polynomials, we modify the definition of the derivatives as follows.

Definition 2.2**.**

The first order new derivative is defined by

[TABLE]

and the new defined derivative of order $k,k\geq 1$ is denoted by

[TABLE]

where $0<\lambda\leq 1$ and $v(x)$ is defined on $\mathbb{R}^{+}$.

Furthermore, the left-side (right-side) limit definitions of the new derivatives are denoted by

[TABLE]

$D^{1}_{\lambda}v(x)$ exists if and only if ${}^{+}D^{1}_{\lambda}v(x)={}^{-}D^{1}_{\lambda}v(x)$, and $D^{1}_{\lambda}v(x)={}^{+}D^{1}_{\lambda}v(x)={}^{-}D^{1}_{\lambda}v(x).$ Obviously when $\lambda=1$, the new defined derivatives become the classical ones.

We set two weight functions as follows:

[TABLE]

[TABLE]

Lemma 2.3** (see [17, 16]).**

The fractional Jacobi polynomials $J^{\alpha,\beta,\lambda}_{n}(x)$ are orthogonal with respect to the weight function $\omega^{\alpha,\beta,\lambda}(x)$, $\alpha,\beta>-1,0<\lambda\leq 1$, i.e.,

[TABLE]

where

[TABLE]

The special case $\alpha=0,\beta=\frac{1}{\lambda}-1$ yields the Müntz Legendre polynomials, which have been investigated in a different way in [23, 6].

Next, we shall show that the fractional Jacobi polynomials are the eigenfunctions of a singular Sturm-Liouville operator $\mathcal{L}^{\alpha,\beta}_{\lambda}$ defined by

[TABLE]

Lemma 2.4** (see [17, 16]).**

The fractional Jacobi polynomials \big{\{}J^{\alpha,\beta,\lambda}_{n}\big{\}}_{n=0}^{\infty} satisfy the following singular Sturm-Liouville problem:

[TABLE]

That is

[TABLE]

where $\sigma^{\alpha,\beta}_{n}=n(n+\alpha+\beta+1).$

Lemma 2.5**.**

The new defined k-th order derivatives of the fractional Jacobi polynomials are orthogonal with respect to the weight $\omega^{\alpha+k,\beta+k,\lambda}(x)$, i.e.,

[TABLE]

where

[TABLE]

Moreover, we have

[TABLE]

where

[TABLE]

Proof.

We start with the case $k=1$. Using integration by parts, Lemma 2.4, and the orthogonality of \big{\{}J^{\alpha,\beta,\lambda}_{n}\big{\}}_{n=0}^{\infty}, we obtain

[TABLE]

This means that \big{\{}D^{1}_{\lambda}J^{\alpha,\beta,\lambda}_{n}\big{\}}_{n=1}^{\infty} are orthogonal with respect to the weight $\omega^{\alpha+1,\beta+1,\lambda}$. Furthermore it follows from Lemma 2.2 that $D^{1}_{\lambda}J^{\alpha,\beta,\lambda}_{n}$ must be proportional to $J^{\alpha+1,\beta+1,\lambda}_{n-1},$ namely,

[TABLE]

where $\widehat{d}^{\alpha,\beta}_{n,1}$ is a constant.

In virtue of (2.2), the leading coefficient of $J^{\alpha,\beta,\lambda}_{n}$ is $\displaystyle k^{\alpha,\beta}_{n}=\frac{\Gamma(2n+\alpha+\beta+1)}{n!\Gamma(n+\alpha+\beta+1)}.$ Thus comparing the leading coefficients on both sides of (2.12) gives:

[TABLE]

This proves (2.10) for $k=1$. Applying (2.12) recursively, we obtain

[TABLE]

where

[TABLE]

This proves (2.10) for all $k\geq 1$. Finally, using Lemma 2.3 gives

[TABLE]

where

[TABLE]

The proof is completed. ∎

2.2. $L^{2}_{\omega^{\alpha,\beta,\lambda}}(I)$-orthogonal projector with $\alpha,\beta>-1$

Let $\pi_{N,\omega^{\alpha,\beta,\lambda}}:L^{2}_{\omega^{\alpha,\beta,\lambda}}(I)\rightarrow P^{\lambda}_{N}(I)$ be the $L^{2}_{\omega^{\alpha,\beta,\lambda}}$-orthogonal projection operator defined by: for all $v\in L^{2}_{\omega^{\alpha,\beta,\lambda}}(I)$, $\pi_{N,\omega^{\alpha,\beta,\lambda}}v\in P^{\lambda}_{N}(I)$ such that

[TABLE]

Equivalently, $\pi_{N,\omega^{\alpha,\beta,\lambda}}$ can be characterized by:

[TABLE]

where $J^{\alpha,\beta,\lambda}_{n}(x)$ are the fractional Jacobi polynomials defined in (2.1), and

[TABLE]

Immediately, we have

[TABLE]

To measure the truncation error, we introduce the non-uniformly fractional Jacobi-weighted Sobolev space:

[TABLE]

equipped with the inner product, norm, and semi-norm as follows:

[TABLE]

The special case $\lambda=1$ gives the classical non-uniformly Jacobi-weighted Sobolev space:

[TABLE]

The functions in $B^{m,\lambda}_{\alpha,\beta}(I)$ and the ones in $B^{m,1}_{\alpha,\beta}(I)$ are linked through the following lemma.

Lemma 2.6**.**

A function $v$ belongs to $B^{m,\lambda}_{\alpha,\beta}(I)$ if and only if $v(x^{1/\lambda})$ belongs to $B^{m,1}_{\alpha,\beta}(I).$

Proof.

Using the variable change $x=t^{1/\lambda}$, we have

[TABLE]

This proves the desired result. ∎

Proposition 2.1**.**

The orthogonal projector $\pi_{N,\omega^{\alpha,\beta,\lambda}}$ admits the following error estimate: for any $v(x^{\frac{1}{\lambda}})\in B^{m,1}_{\alpha,\beta}(I),$ and $0\leq l\leq m\leq N+1,$

[TABLE]

For a fixed $m$, the above estimate can be simplified as

[TABLE]

where $c\approx 1$ for $N\gg 1.$ In particular, for $l=0,1$, we have

[TABLE]

[TABLE]

where $\tilde{\omega}^{\alpha,\beta,\lambda}=\lambda^{-1}(1-x^{\lambda})^{\alpha+1}x^{\beta\lambda+1}$.

Proof.

For $v(x^{\frac{1}{\lambda}})\in B^{m,1}_{\alpha,\beta}(I)$, Lemma 2.6 gives $v(x)\in B^{m,\lambda}_{\alpha,\beta}(I)$. Thanks to (2.13) and the orthogonality (2.8), we have

[TABLE]

Further estimation on the right hand side term gives

[TABLE]

In virtue of (2.9), we find that for $0\leq l\leq m\leq N+1$,

[TABLE]

where we have used the facts: $\alpha+\beta+2>0.$ Then the estimate (2.15) follows from (2.19), (2.20), and (2.14).

Next, we prove (2.16). Recall the property of the Gamma function (see [1, (6.1.38)]):

[TABLE]

Moreover it can be shown that for any constant $a,b\in\mathbb{R},n+a>1$ and $n+b>1$ (see [32, Lemma 2.1],[10]),

[TABLE]

where

[TABLE]

Hence, we obtain

[TABLE]

where $\nu^{2-m,2-l}_{N}\approx 1$ for fixed $m$ and $N\gg 1.$

Combining (2.15) and (2.21) gives (2.16). Furthermore, using (2.16) and the definition of $D^{1}_{\lambda}$ we obtain (2.17) and (2.18). ∎

2.3. $I^{\alpha,\beta}_{N,\lambda}$-interpolation on fractional Jacobi-Gauss-type points with $\alpha,\beta>-1$

Let $h^{\alpha,\beta}_{j,\lambda}(x)$ be the generalized Lagrange basis function:

[TABLE]

where $x_{0}<x_{1}<\dots<x_{N-1}<x_{N}$ are zeros in $I$ of $J^{\alpha,\beta,\lambda}_{N+1}(x)$. It is clear that the functions $h^{\alpha,\beta}_{j,\lambda}(x)$ satisfy

[TABLE]

Let $z(x)=x^{\lambda}$. Then $z_{i}:=z(x_{i})=x_{i}^{\lambda},~{}~{}0\leq i\leq N$, are zeros of $J^{\alpha,\beta,1}_{N+1}(x)$, and

[TABLE]

We define the generalized interpolation operator $I^{\alpha,\beta}_{N,\lambda}$ by

[TABLE]

Then

[TABLE]

Lemma 2.7** (Case $\lambda=1$, [27] Lemma 3.8).**

The interpolation operator $I^{\alpha,\beta}_{N,1}$ satisfies following stability estimate: for any $v\in B^{1,1}_{\alpha,\beta}(I)$,

[TABLE]

Proposition 2.2**.**

For any v(x^{\frac{1}{\lambda}})\in{\color[rgb]{0,0,0}B^{1,1}_{\alpha,\beta}(I)}, we have

[TABLE]

Proof.

Making the variable change $z=x^{\lambda}$ in the definition of $I^{\alpha,\beta}_{N,\lambda}v(x)$, we get

[TABLE]

Then it follows from Lemma 2.7 and (2.26)

[TABLE]

The desired result (2.25) follows from combining (2.27) and (2.14). ∎

We will also need some inverse inequalities for the $\lambda-$polynomials, which are given in the following proposition.

Proposition 2.3**.**

For any $\phi\in P^{\lambda}_{N}(I)$, we have

[TABLE]

and

[TABLE]

where $\sigma^{\alpha,\beta}_{N}$ is defined in (2.7) and $c\approx 1$ for fixed $k$ and $N\gg 1$.

Proof.

For any $\phi\in P^{\lambda}_{N}(I)$, we express it under the form

[TABLE]

where $\hat{\gamma}_{i}^{\alpha,\beta}$ was defined in (2.5). Then we deduce from the orthogonality of $\{J^{\alpha,\beta,\lambda}_{i}(x)\}_{i=0}^{N}$

[TABLE]

Applying the differential operator $D^{1}_{\lambda}$ to the both sides of (2.30), and using the orthogonality (2.11), we get

[TABLE]

This results in (2.28).

Applying the operators $D^{1}_{\lambda}$ $k$ times to (2.30) and using (2.8), we obtain

[TABLE]

Similar to the discussion for (2.21), we have

[TABLE]

where $\nu^{k+\alpha+\beta+1,\alpha+\beta+1}_{N}\nu^{1,1-k}_{N}\approx 1$ for fixed $k$ and $N\gg 1.$

Then the inequality (2.29) follows from combining (2.31) and (2.32). ∎

We now present the error estimation for the fractional Jacobi-Gauss interpolation operator based on the roots of $J^{\alpha,\beta,\lambda}_{N+1}(x)$.

Proposition 2.4**.**

It holds for any $v(x^{1/\lambda})\in B^{m,1}_{\alpha,\beta}(I),m\geq 1$, and any $0\leq l\leq m\leq N+1$,

[TABLE]

If $m$ is fixed, then the above estimate can be simplified as

[TABLE]

where $c\approx 1$ for $N\gg 1.$ In particular, for $l=0,1$ we have

[TABLE]

[TABLE]

Proof.

From (2.15) and Proposition 2.2, we deduce

[TABLE]

Furthermore, using the inverse inequality (2.29), we obtain

[TABLE]

The above inequality, together with the triangle inequality and (2.15) gives

[TABLE]

This proves the first estimate. In case $m$ is fixed, the estimates (2.33)-(2.35) can be derived in a similar way as in Proposition 2.1. The proof is completed. ∎

Lemma 2.8** (Weighted Sobolev inequality).**

For all $v(x)\in B^{1,\lambda}_{\alpha,\beta}(I),v(\xi)=0$ for some $\xi\in[0,1]$, and any $-1<\alpha,\beta<-\frac{1}{2}$, it holds

[TABLE]

Proof.

Using the Cauchy-Schwarz inequality, we have, for all $x\in[\xi,1]$,

[TABLE]

The condition $-1<\alpha,\beta<-\dfrac{1}{2}$ guarantees that it holds for any $x\in I$,

[TABLE]

Bringing (2.39) into (2.38) gives

[TABLE]

In a similar way we can prove

[TABLE]

This completes the proof. ∎

Proposition 2.5** (Interpolation error in $L^{\infty}-$norm).**

If $-1<\alpha,\beta\leq-\dfrac{1}{2},$ we have

[TABLE]

Proof.

Obviously if $v(x^{1/\lambda})\in B^{m,1}_{\alpha,\beta}(I),m\geq 1$, then $v(x)\in B^{1,\lambda}_{\alpha,\beta}(I)$. By definition, the function $v(x)-I^{\alpha,\beta}_{N,\lambda}v(x)$ vanishes at all the roots of $J^{\alpha,\beta,\lambda}_{N+1}(x)$. Thus we can use Lemma 2.8 and Proposition 2.4 to conclude

[TABLE]

This proves the proposition. ∎

We use the Jacobi Gauss points to define the discrete inner product: for any $u,v\in C(\bar{I})$,

[TABLE]

where $\{z_{i}\}_{i=0}^{N}$ are the zeros of the shifted Jacobi polynomial $J^{\alpha,\beta,1}_{N+1}(x)$ and $\{\omega_{i}\}_{i=0}^{N}$ are the corresponding weights.

Lemma 2.9** ([27, 5]).**

For all $v\in B^{m,1}_{\alpha,\beta}(I),m\geq 1$ and all $\phi\in P^{1}_{N}(I)$, we have

[TABLE]

The following result, which can be found in [22], concerns the Lebesgue constant of the Lagrange interpolation polynomials associated with the zeros of shifted Jacobi polynomials.

Lemma 2.10**.**

Let $\{h^{\alpha,\beta}_{j,1}(x)\}_{j=0}^{N}$ be the Lagrange interpolation polynomials associated with the Gauss points of the shifted Jacobi polynomial $J^{\alpha,\beta,1}_{N+1}(x).$ Then

[TABLE]

The above classical results can be readily extended to the generalized Lagrange interpolation polynomials associated with the zeros of fractional Jacobi polynomials, which is given in the following lemma.

Lemma 2.11**.**

Let $\{h^{\alpha,\beta}_{j,\lambda}(x)\}_{j=0}^{N}$ be the generalized Lagrange interpolation basis functions associated with the Gauss points of the fractional Jacobi polynomial $J^{\alpha,\beta,\lambda}_{N+1}(x).$ Then

[TABLE]

From now on, for $r\geq 0$ and $\kappa\in[0,1],$ $C^{r,\kappa}(I)$ will denote the space of functions whose $r-$th derivatives are Hölder continuous with exponent $\kappa$, endowed with the usual norm:

[TABLE]

When $\kappa=0,C^{r,0}(I)$ denotes the space of functions with $r$ continuous derivatives on $I$, which is also commonly denoted by $C^{r}(I)$, endowed with the norm $\|\cdot\|_{r}$.

A well-known result from Ragozin [24, 25] will be useful, which states that, for non-negative integer $r$ and real number $\kappa\in(0,1)$, there exists a linear operator $\mathcal{T}_{N}$ from $C^{r,\kappa}(I)$ into $P_{N}^{1}(I)$, such that

[TABLE]

where $C_{r,\kappa}$ is a constant which may depend on $r$ and $\kappa$.

We further define the linear, weakly singular integral operator $\mathcal{K}$:

[TABLE]

where $K\in C(I\times I)$ with $K(x,x)\neq 0$ for $x\in I$. We will prove that $\mathcal{K}$ is a compact operator from $C(I)$ to $C^{0,\kappa}(I)$ for any $0<\kappa<1-\mu.$

Lemma 2.12**.**

For any function $v\in C(I)$, $K\in C(I\times I)$, and $K(\cdot,s)\in C^{0,\kappa}(I)$ with $0<\kappa<1-\mu$, we have

[TABLE]

This implies that

[TABLE]

Proof.

Without loss of generality, assume $0\leq y<x\leq 1$. We have

[TABLE]

where

[TABLE]

For $M_{1}$, we use the triangle inequality to deduce

[TABLE]

where

[TABLE]

$M^{(1)}$ and $M^{(2)}$ can be bounded respectively by

[TABLE]

For $M_{2},$ we have

[TABLE]

Gathering all these estimates thus gives (2.44). ∎

We can also prove the following result.

Lemma 2.13**.**

For any function $v(x)\in C(I),$ $K\in C(I\times I)$, and $K(\cdot,s)\in C^{0,\kappa}(I)$ with $0<\kappa<1-\mu$, there exists a positive constant $c$ such that

[TABLE]

Thus

[TABLE]

Proof.

It follows from (2.44) that

[TABLE]

Noticing that, for any $x,y\in I,x\neq y,0<\lambda\leq 1$, we have

[TABLE]

Thus we obtain

[TABLE]

Combining with (2.47) and (2.48) gives (2.46). Thus we complete the proof. ∎

3. Fractional Jacobi Spectral-Collocation Method and Convergence Analysis

This section is devoted to developing and analyzing an efficient method for the following equation: $0<\mu<1$,

[TABLE]

where $\mathcal{K}(\cdot)$ is the weakly singular integral operator defined in (2.43), $g(x)\in C(I)$, $K\in C(I\times I),K(x,x)\neq 0$ for $x\in I$.

3.1. Fractional Jacobi Spectral-Collocation Method for VIEs

We consider the fractional Jacobi spectral-collocation method as follows: find fractional polynomial $u_{N}^{\lambda}\in P^{\lambda}_{N}(I)$, such that

[TABLE]

where the collocation points $\{x_{i}\}_{i=0}^{N}$ are roots of $J^{\alpha,\beta,\lambda}_{N+1}(x)$.

Since the exact evaluation of \big{(}\mathcal{K}\varphi\big{)}(x_{i}) is not realizable in pratical cases, we need to find an efficient way to approximate this term for all points $\{x_{i}\}_{i=0}^{N}$ without losing the high order accuracy of the scheme. For a given collocation point $x_{i}$, we rewrite the integral term by using the variable change $s=\tau_{i}(\theta):=x_{i}\theta^{\frac{1}{\lambda}}$:

[TABLE]

Let

[TABLE]

Then

[TABLE]

We choose to approximate the integral $(\bar{K}(x_{i},\tau_{i}(\theta)),\varphi(\tau_{i}(\theta)))_{\omega^{-\mu,1/\lambda-1,1}}$ by the numerical quandrature $(\bar{K}(x_{i},\tau_{i}(\theta)),\varphi(\tau_{i}(\theta)))_{N,\omega^{-\mu,1/\lambda-1,1}}$, which is defined by

[TABLE]

where $\{\theta_{j}\}_{j=0}^{N}$ are the zeros of $J^{-\mu,1/\lambda-1,1}_{N+1}(\theta)$ and $\{\omega_{j}\}_{j=0}^{N}$ are the corresponding weights. For the sake of simplification, we will denote

[TABLE]

This leads us to consider the following discrete problem: find $u_{N}^{\lambda}(x):=\sum_{i=0}^{N}u_{i}h_{i,\lambda}^{\alpha,\beta}(x)\in P^{\lambda}_{N}(I)$, such that

[TABLE]

where $\{h_{i,\lambda}^{\alpha,\beta}\}_{i=0}^{N}$ are the fractional Lagrange polynomials defined in (2.22).

3.2. Convergence analysis

The purpose of this subsection is to analyze the discrete problem (3.5) and derive error estimates for the discrete solution. First we derive the error estimate in the $L^{\infty}-$norm.

Theorem 3.1**.**

Let $u(x)$ be the exact solution to the Volterra integral equation (3.1) and $u_{N}^{\lambda}(x)$ is the numerical solution of the fractional Jacobi spectral-collocation problem (3.5). Assume $0<\mu<1,-1<\alpha,\beta\leq-\frac{1}{2}$, $K(x,s)\in C^{m}(I,I)$ and $u(x^{\frac{1}{\lambda}})\in B^{m,1}_{\alpha,\beta}(I),m\geq 1$. Then we have

[TABLE]

where

[TABLE]

Proof.

Let $e(x)=u(x)-u_{N}^{\lambda}(x)$ be the error function. Subtracting (3.5) from (3.1) gives the error equation:

[TABLE]

where $e_{i}=u(x_{i})-u_{i}$, and $q_{i}$ is the quadrature error term

[TABLE]

which can be bounded by using Lemma 2.9

[TABLE]

Multiplying both sides of (3.8) by $h^{\alpha,\beta}_{i,\lambda}(x)$ and summing up the resulting equation from $i=0$ to $i=N$ gives

[TABLE]

Rearranging this equation leads to

[TABLE]

where

[TABLE]

which correspond respectively the interpolation error, numerical quadrature error, and the interpolation error for the integral operator $\mathcal{K}$. Using Gronwall lemma gives

[TABLE]

where $K_{0}=\displaystyle\max_{0\leq s<x\leq 1}|K(x,s)|.$ We then obtain

[TABLE]

where c=1+\big{(}K_{0}/(1-\mu)\big{)}\exp(K_{0}/(1-\mu)\big{)}. Now we bound the right hand side term by term. For the first term, it follows from Proposition 2.5:

[TABLE]

For $I_{2}$, we deduce from (3.9)

[TABLE]

This, together with Lemma 2.11, gives

[TABLE]

It remains to estimate the third term $I_{3}$. It follows from (2.24), Lemma 2.13, Lemma 2.11, and (2.42) that

[TABLE]

The desired result then follows from bringing the estimates (3.15), (3.17), and (3.18) into (3.14). This completes the proof. ∎

Remark 3.1**.**

*Note that the error estimate (3.6) involves the quantity $K^{*}$ defined in (3.7). The boundedness of this quantity depends on the regularity of the function $\bar{K}$, which is a transformation of $K$ through (3.2). Thus it is important to observe that if we set $\lambda=1/p$ with $p$ being a positive integer, then under the initial regularity assumption on $K(x,s)$, i.e., $K(x,s)\in C^{m}(I,I)$, we have for all $x_{i}\in I$, $h(\theta):=K(x_{i},x_{i}\theta^{1/\lambda})=K(x_{i},x_{i}\theta^{p})\in C^{m}(I)$. Furthermore, noticing that \Big{(}\dfrac{1-\theta^{1/\lambda}}{1-\theta}\Big{)}^{-\mu}\in C^{\infty}(I), we then have $\bar{K}(x_{i},\tau_{i}(\theta))\in C^{m}(I)$. This guarantees the boundedness of $K^{*}$ in (3.7).

It is also notable that the constant $c$ in front of the estimate (3.14) for $\|u-u_{N}^{\lambda}\|_{\infty}$ will blow up when $\mu\rightarrow 1$. We have tried but unfortunately found no way to improve the estimate given in the proof of Theorem 3.1. Possible improvement for this constant certainly requires new estimation technique, and is worth further investigation.*

To derive the error estimate in the weighted $L^{2}-$norm, we will make use of the following generalized Hardy’s inequality with weights (see [19, 15]).

Lemma 3.1**.**

For all measurable function $f\geq 0$, weight functions $u$ and $v$, $1<p\leq q<\infty$, $-\infty\leq a<b\leq\infty$, the following generalized Hardy’s inequality holds

[TABLE]

if and only if

[TABLE]

where $T$ is an operator of the form

[TABLE]

with $\rho(x,s)$ being a given kernel.

From Lemma 4.2 in [11] and (2.26), we have the following weighted mean convergence result for the fractional interpolation operator $I^{\alpha,\beta}_{N,\lambda}$.

Lemma 3.2**.**

For any bounded function $v(x)$ defined on I, there exists a constant $c$ independent of $v$ such that

[TABLE]

Theorem 3.2**.**

Let $u(x)$ be the exact solution to (3.1), $u_{N}^{\lambda}(x)$ is the solution of discrete problem (3.5). If $0<\mu<1,-1<\alpha,\beta\leq-\frac{1}{2}$, $K(x,s)\in C^{m}(I,I)$, and $u(x^{\frac{1}{\lambda}})\in B^{m,1}_{\alpha,\beta}(I),m\geq 1$, then we have

[TABLE]

where $K^{*}$ is defined in (3.7).

Proof.

Similar to (3.13), we deduce from Lemma 3.1

[TABLE]

To bound the right hand side, we first use Proposition 2.4 to get

[TABLE]

Then applying Lemma 3.2 and (3.16) gives

[TABLE]

Finally, we notice $\|\cdot\|_{\omega^{\alpha,\beta,1}}$ can be bounded by $\|\cdot\|_{\infty}$, and

[TABLE]

where $\{z_{i}\}_{i=0}^{N}$ are zeros of shifted Jacobi polynomial $J^{\alpha,\beta,1}_{N+1}(z).$ Thus it follows from (2.42), Lemma 2.13, and Lemma 3.2

[TABLE]

Then we deduce from Theorem 3.1

[TABLE]

Combining (3.21), (3.22), and (3.24) proves the theorem. ∎

4. Numerical results

In this section, we present some numerical results to verify the error estimates obtained in the previous section for the proposed method. The main goal is to demonstrate that the new method has the capability to capture the typical solutions of VIEs with high accuracy. To this end a series of numerical tests are carried out in the following. All the errors to be presented below are in $L^{\infty}(I)-$ and $L^{2}_{\omega^{\alpha,\beta,\lambda}}(I)-$norms, where $\alpha,~{}\beta$ are related to the corresponding collocation points, i.e., the zeros of $J^{\alpha,\beta,\lambda}_{N+1}(x)$. In order to calculate the errors, for the test cases in which the exact solution is not available, a numerical solution obtained with a very large $N$ will be served as the exact solution.

Example 4.1**.**

Consider the Volterra integral equation (3.1) with $g(x)=1$ and $K(x,s)=\exp(x-s).$

It has been known from [9] that for this smooth $g$ and kernel function $K$, the solution of (3.1) can be expressed as

[TABLE]

where $G:=\{(j,k):j,k$ are non-negative integers}, $\gamma_{j,k}$ are constants, and $u_{r}$ is a smooth function.

In virtue of the solution structure, we want to choose the value of $\lambda$ such that $u(x^{1/\lambda})$ is smooth or as regular as possible. This can be easily done by taking $\lambda$ such that $(j+k\mu)/\lambda$ is integer or as large as possible. In this first example, we test for $\mu=0.2$ and $0.5$. In these cases there are many possible choices for $\lambda$. The obtained numerical results are shown in semi-log scale as a function of $N$ in Figure 1 with $\alpha=\beta=-\frac{1}{2}$. It is clear that the exponential convergence rate is attained with $\lambda=1/5,1/10,1/15$ for $\mu=0.2$, and $\lambda=1/2,1/4,1/6$ for $\mu=0.5$ since all these values make $(j+k\mu)/\lambda$ integer.

Example 4.2**.**

Consider the equation (3.1) with $K(x,s)=1,$ and limited regular source term $g(x)=x^{0.5}$.

In this case, there is not available result on the regularity of the exact solution. So it is interesting to see if high accuracy is attainable for the fractional spectral method by choosing suitable $\lambda$. Figure 2 shows the errors in semi-log scale as a function of $N$ for a number of values of $\mu$ and $\lambda$. In computing the errors we have assumed that the numerical solution obtained with $N=80$ is the “exact” solution since the real solution is unknown. Surprisingly, for all tested values of $\mu$ and $\lambda$, the error curves are almost linear versus the degrees of the $\lambda$-polynomial. This means that the convergence rate of the proposed method is exponential even if there is no information about the exact solution. A reasonable explanation for this excellent result is that the transformed solution, i.e., $u(x^{1/\lambda})$, becomes smooth or regular enough if a suitably small $\lambda$ is used in the approximation.

Example 4.3**.**

([21] Example 5.4) Consider the fabricated solution $u(x)=x^{-\mu}\sin(x)$, corresponding to $K(x,s)=1$ and

[TABLE]

where $\mathcal{B}(\cdot,\cdot)$ is the Bessel function, i.e.,

[TABLE]

Clearly, the exact solution of this example has singularity at the left end point, i.e., $u^{\prime}(x)\sim x^{-\mu}$ at $x=0$. In Figure 3, we plot the errors in the $L^{\infty}(I)$- and $L^{2}_{\omega^{-1/2,-1/2,\lambda}}$-norm in semi-log scale as a function of $N$ for $\mu=\frac{1}{2},\frac{2}{3}$. In the calculation we have used $\lambda=\frac{1}{2}$ and $\frac{1}{4}$ for $\mu=\frac{1}{2}$, and $\lambda=\frac{1}{3}$ and $\frac{1}{6}$ for $\mu=\frac{2}{3}$. We observe that the exponential convergence is obtained by the proposed fractional spectral method, which is an expected result since for the tested $\mu$, $u(x^{1/\lambda})$ is smooth for the $\lambda$ values used in the computation. It is worth to notice that this example has also been investigated in [21] by using a classical spectral method. The strategy used in that paper is a smoothing transformation approach. That is, the original integral equation is first transformed into a new one having smooth solution, then a pseudo spectral method based on the classical Chebyshev or Legendre polynomials is constructed to approximate the smooth solution. We emphasize that the current paper makes use of a new idea, which is completely different than the smoothing transformation approach used in [21]. The new idea allows full analysis for the proposed numerical schemes, and is applicable to a wider class of problems.

Example 4.4**.**

Still consider a fabricated solution of (3.1) with $K(x,s)=1$ and the exact solution as follows: $(i)$ $u(x)=x^{\gamma_{1}}+x^{\gamma_{2}}$; $(ii)$ $u(x)=\sin(x^{\gamma_{1}}+x^{\gamma_{2}}).$

This solution, which looks quite unrealistic, is served as a good example to test the capability of the proposed method in difficult situation. It is clear that for two general parameters $\gamma_{1}$ and $\gamma_{2}$, it is not always possible to make $u(x^{1/\lambda})$ smooth. However, a careful examination shows that $u(x^{1/\lambda})\in B^{2\gamma/\lambda+\beta+1-\epsilon}_{\omega^{\alpha,\beta,1}}(I)$ for any $\epsilon>0$, where

[TABLE]

The numerical errors in log-log scale for a number of different $\gamma_{1},\gamma_{2},\mu,\lambda$, and $\beta$ are presented in Figure 4. Since $K(x,s)$ is smooth, the second terms in the error estimates (3.6) and (3.19) are expected to be negligible as compared with the first terms in these estimates for large $N$. As a consequence, the result shown in the figure should reflect the impact of the regularity of $u(x^{1/\lambda})$ on the error behavior. To closely observe the error decay rates, the $N^{-(2\gamma/\lambda+\beta+1)}$ and $N^{-(2\gamma/\lambda+\beta+1)+1/2}$ decay rates are also plotted in the figure. We have the following observations from this test:

1. all the error curves are straight lines in this log-log representation, which indicates the algebraic convergence. This is consistent with the limited regularity of $u(x^{1/\lambda})$;
2. the convergence rate is in a good agreement with the theoretical prediction given in (3.6) and (3.19), i.e., order $N^{-(2\gamma/\lambda+\beta+1)+1/2}$ and $N^{-(2\gamma/\lambda+\beta+1)}$ respectively;
3. decreasing $\lambda$ increases the regularity of $u(x^{1/\lambda})$, thus results in improvement of the convergence.

5. Concluding remarks

In this work we investigated numerical solutions of the Volterra integral equations with the weakly singular kernel $(x-s)^{-\mu}$, $0<\mu<1$. The main difficulty in approximating this equation is that the derivative of the solution is singular. This has resulted in low order convergence for any methods using the traditional polynomials. To overcome this difficulty, we established a framework of approach based on the fractional Jacobi polynomials for the Volterra integral equations. Precisely, a fractional Jacobi spectral-collocation method is constructed and analyzed for the underlying equation. The significance of the approach is that it makes use of the fractional polynomials as the approximation space, which can well capture typical solution structures of the singular integral equations. It was known that a typical solution of the Volterra integral equations with weakly singular kernel behaves like a series of power functions $x^{\gamma}$. We have seen from the presented analysis that such functions can be much better approximated by the space span$\{1,x^{\lambda},x^{2\lambda},\dots,x^{N\lambda}\}$ than the traditional polynomial space span$\{1,x,x^{2},\dots,x^{N}\}$ as long as $\lambda$ is suitably chosen. In order to carry out a rigorous error analysis, we first established some approximation results for the weighted projection and interpolation operators. Then we derived the error estimates in the $L^{\infty}-$ and weighted $L^{2}-$norms for the proposed method. A series of numerical examples were carried out to verify the theoretical claims. The most remarkable property of the new method is its capability to achieve spectral convergence for solutions of limited regularity. It is worth mentioning that the choice of $\lambda$ is also of importance for the efficiency of the new method. Although there does not exist optimal choice for $\lambda$ for general problems, it can be made according to the following strategy: in case the regularity of the exact solution is unavailable, the parameter $\lambda$ can be taken like $1/q$ with moderately large integer $q$ so that $u(x^{q})$ is smooth enough. Our numerical experiments have shown that doing this can increase the convergence rate about $q$ times.