Piecewise Fractional Jacobi Polynomial Approximations for Volterra Integro-Differential Equations with Weakly Singular Kernels

Abstract

This paper is concerned with numerical solutions to Volterra integro-differential equations with weakly singular kernels. Making use of the transformed fractional Jacobi polynomials, we develop a class of piecewise fractional Galerkin methods for solving this kind of Volterra equation. Then, we study the existence, uniqueness and convergence properties of Galerkin solutions by exploiting the decaying rate of the coefficients of the transformed fractional Jacobi series. Finally, numerical experiments are carried out to illustrate the performance of the piecewise Galerkin solution.

1. Introduction

In recent years, considerable interest in the field of numerical solutions to Volterra integro-differential equations (VIDEs) with singular kernels has been stimulated by their wide application [1] and close relation to fractional differential equations. Consider the following VIDE with a weakly singular kernel:

$$\begin{array}{c}\hfill y^{\prime }\left(t\right)=p\left(t\right)y\left(t\right)+q\left(t\right)+{\int }_0^t{(t-s)}^{\mu }K(t,s)y\left(s\right)ds,t\in [0,T],\end{array}$$

(1)

where $y\left(0\right)=y_0$ denotes the initial value, $p\left(t\right)$ and $q\left(t\right)$ are continuous on $[0,T],$ $K(t,s)$ is continuous in the domain $D:=\left\{\right(t,s):0\le s\le t\le T\}$ and $\mu \in (-1,0).$ Without loss of generality, we assume $T=1.$ According to [2] (Theorem 7.1.1), the VIDE in Equation (1) possesses a unique solution $y\left(t\right)\in C\left(\right[0,1\left]\right).$ Aside from that, the regularity of the derivative of $y\left(t\right)$ is given by

$$|y^{″}\left(t\right)|\le B{t}^{\mu },t\in (0,1],$$

which implies the second derivative of $y\left(t\right)$ is unbounded at $t=0.$ In the remaining parts, we denote various constants to be B for simplicity. The singularity of $y^{″}\left(t\right)$ presents difficulties in the construction of high-order numerical methods. In fact, standard algorithms with uniform grids usually result in a low convergence rate $O\left(h^{\mu +1}\right)$, with h denoting the stepsize ([2], Theorem 7.2.3).

To address this thorny difficulty, graded meshes are employed by several authors. In [3], Brunner, Pedas and Vainikko transformed the original VIDE (Equation (1)) into the second type of Volterra integral equation by employing integration and Dirichlet’s formula. By applying collocation methods to the transformed Volterra integral equations at the graded meshes, efficient numerical solutions to the weakly singular problem were constructed in Equation (1). Furthermore, a comprehensive theory of the optimal convergence estimates was established. Let $[0,1]$ be partitioned by grid points ${\displaystyle t_j={\left(\frac{j}{N}\right)}^r}$ with ${\displaystyle r>\frac{m}{1+\mu },}$ where m denotes the number of collocation points in each subinterval. Then, the convergence order of the collocation solution to the VIDE in Equation (1) could attain $O\left(h^m\right).$ Similarly, by utilizing geometrically graded meshes, Brunner and Schötzau analyzed the discontinuous Galerkin method for the VIDE in Equation (1) in [4]. By imposing particular bounds for $p\left(t\right)$ and the initial data, the Galerkin approximation guaranteed an exponential convergence. Parabolic Volterra integro-differential equations with weakly singular kernels were considered in [5], where error estimates with respect to all parameters were developed. By using non-uniformly refined time steps, Mustapha devised superconvergent discontinuous Galerkin approximations in [6]. Moreover, continuous Galerkin methods were discussed by Yi and Guo in [7], where algebraic convergence rates were achieved with quasi-uniform meshes. Its extension to Volterra delay-integro-differential equations was considered in [8]. Based on the reducible quadrature developed from the boundary value technique for solving the differential equation, the block boundary value method was investigated in [9].

Although numerical algorithms at graded meshes are able to increase the convergence rate, clustering of the collocation points near $t=0$ makes themselves suffer round-off errors if the graded parameter r is large. An applicable alternative approach is to employ fractional polynomials or variable transformation, which has gained more and more attention in last several years. In [10], Diogo et al. studied the method of smoothing the solutions of VIDEs with weakly singular kernels. The solutions were regularized by the transformation such that $\phi \left(t\right)=t^d,$ and piecewise collocation methods with a graded mesh were employed. Let r be the graded parameter. Then, the magnitude of the transformed collocation error $e_N$ enjoyed the following estimates (see [10], Theorem 4.2):

• If $d(1+\mu )<m<1+d(1+\mu ),$ and $r\ge 1,$ then

  $$e_N=O\left(N^{-m}\right),N\to \infty .$$
• If $m=1+d(1+\mu ),$ then

  $$e_N=\left\{\begin{array}{cc}O\left(N^{-m}(1+logN)\right),\hfill & r=1,\hfill \\ O\left(N^{-m}\right),\hfill & r>1.\hfill \end{array}\right.$$
• If $m>1+d(1+\mu ),$ then

  $$e_N=\left\{\begin{array}{cc}O\left(N^{-r(1+d(1+\mu \left)\right)}\right),\hfill & 1\le r<m/(1+d(1+\mu \left)\right),\hfill \\ O\left(N^{-m}(1+logN)\right),\hfill & r=m/(1+d(1+\mu \left)\right),\hfill \\ O\left(N^{-m}\right),\hfill & r>m/(1+d(1+\mu \left)\right).\hfill \end{array}\right.$$

A similar technique was extended to the numerical solution of the weakly singular Volterra integral equations in two dimensions (see [11]). It can be seen from the above theoretical results that the regularized algorithm results in pretty fast convergence rates, even in the case where the weak singularity is not known. With the help of fractional Jacobi polynomials, a class of fractional spectral Galerkin methods for VIDEs was discussed by Hou and Xu in [12]. Under restrictions on the coefficient of $y\left(t\right)$ and the kernel function $K(t,s),$ the existence and uniqueness of the Galerkin solution and corresponding error estimates were established. The fractional polynomial also plays an important role in solving other weakly singular problems. In [13], Cai and Chen proposed a class of spectral collocation methods with the help of fractional Lagrange interpolation. Aside from that, they considered the conditioning number of the discretized linear system and gave the error estimates with respect to special Chebyshev-type weight functions. Its further extension to nonlinear problems was studied in [14]. By utilizing the zeros of fractional Jacobi polynomials, Hou, et al. discussed the fractional spectral method in [15], where detailed convergence properties with regard to the $L^{\infty }$- and weighted $L^2$-norms were proposed. Through transforming the initial value problem into the boundary value problem with approximated end values, the authors developed a kind of fractional collocation boundary value method for solving the second type of Volterra integral equations with weakly singular kernels in [16]. Fractional Jacobi polynomials also found application in the calculation of highly oscillatory integrals. In [17], a special kind of Petrov–Galerkin method for solving Levin’s equation was discussed, which led to an efficient sparse fractional Jacobi–Galerkin–Levin quadrature rule. The third kind of Volterra integral equation was studied with the help of fractional interpolation in [18], where numerical studies indicated that the fractional collocation method provided more accurate approximations than the graded mesh did.

This paper is devoted to studying the piecewise fractional Galerkin approximation to the solutions of VIDEs (Equation (1)). The remaining parts of this paper are organized as follows. In Section 2, we give the description of piecewise fractional Galerkin methods (PFG). Then, the solvability and convergence property of PFG are discussed in Section 3. Numerical experiments are conducted in Section 4 to illustrate the performance of PFG. In the final section, we conclude with some remarks.

2. Formulation of PFG

In this section we will introduce transformed fractional Jacobi polynomials and construct PFG for a VIDE (Equation (1)).

By letting ${\omega }^{\alpha ,\beta }\left(x\right)={(1-x)}^{\alpha }{(1+x)}^{\beta }$ denote the weight function, we obtain Jacobi polynomials ${\displaystyle {\left\{P_n^{\alpha ,\beta }\left(x\right)\right\}}_{n=0}^{\infty }}$, which are orthogonal with respect to the weight function ${\omega }^{\alpha ,\beta }\left(x\right)$ on the interval $[-1,1]$; that is, we have

$$\begin{array}{c}\hfill {\int }_{-1}^1{\omega }^{\alpha ,\beta }\left(x\right)P_j^{\alpha ,\beta }\left(x\right)P_k^{\alpha ,\beta }\left(x\right)dx={\gamma }_j^{\alpha ,\beta }{\delta }_{j,k},\end{array}$$

(2)

where (see [19], p. 73)

$$\begin{array}{c}\hfill {\delta }_{j,k}=\left\{\begin{array}{c}{\displaystyle 1,j=k,}\hfill \\ {\displaystyle 0,j\ne k,}\hfill \end{array}\right.{\gamma }_j^{\alpha ,\beta }=\frac{2^{\alpha +\beta +1}\mathsf{\Gamma }(j+\beta +1)\mathsf{\Gamma }(j+\alpha +1)}{(2j+\alpha +\beta +1)\mathsf{\Gamma }(j+1)\mathsf{\Gamma }(j+\alpha +\beta +1)}.\end{array}$$

To handle weakly singular problems, it is convenient to resort to a class of high-order Jacobi approximations by extending the definition of Müntz–Legendre polynomials ([20,21]) (i.e., the transformed fractional Jacobi polynomials):

Definition 1.

The transformed fractional Jacobi polynomial of a degree n over the interval $[a,b]$ with $a\ge 0$ is defined as

$$\begin{array}{c}\hfill P_{n,a,b}^{\lambda ,\alpha ,\beta }\left(x\right):=P_n^{\alpha ,\beta }\left(2\frac{x^{1/\lambda }-a^{1/\lambda }}{b^{1/\lambda }-a^{1/\lambda }}-1\right),\forall x\in [a,b],\end{array}$$

(3)

where $\lambda \in [1,\infty )$ denotes the regularization parameter and $\alpha ,\beta >-1.$ In the case of $\lambda =1,$ $P_{n,a,b}^{1,\alpha ,\beta }\left(x\right)$ reduces to the classical transformed Jacobi polynomial.

Noting that

$$\begin{array}{cc}\hfill & {\int }_a^b{\lambda }^{-1}{(b^{1/\lambda }-x^{1/\lambda })}^{\alpha }{(x^{1/\lambda }-a^{1/\lambda })}^{\beta }{x}^{1/\lambda -1}{P}_{j,a,b}^{\lambda ,\alpha ,\beta }\left(x\right)P_{k,a,b}^{\lambda ,\alpha ,\beta }\left(x\right)dx\hfill \\ \hfill =& {\int }_{a^{1/\lambda }}^{b^{1/\lambda }}{(b^{1/\lambda }-t)}^{\alpha }{(t-a^{1/\lambda })}^{\beta }{P}_{j,a,b}^{\lambda ,\alpha ,\beta }\left(t^{\lambda }\right)P_{k,a,b}^{\lambda ,\alpha ,\beta }\left(t^{\lambda }\right)dt\hfill \\ \hfill =& {\left(\frac{b^{1/\lambda }-a^{1/\lambda }}{2}\right)}^{1+\alpha +\beta }{\int }_{-1}^1{(1-s)}^{\alpha }{(1+s)}^{\beta }{P}_j^{\alpha ,\beta }\left(s\right)P_k^{\alpha ,\beta }\left(s\right)ds,\hfill \end{array}$$

where variable transformations

$${\displaystyle x=t^{\lambda }\mathrm{and}t=\frac{b^{1/\lambda }-a^{1/\lambda }}{2}s+\frac{a^{1/\lambda }+b^{1/\lambda }}{2}}$$

are employed, we arrive at the following lemma with the help of Equation (2):

Lemma 1.

The fractional Jacobi polynomials ${\left\{P_{n,a,b}^{\lambda ,\alpha ,\beta }\left(x\right)\right\}}_{n=0}^{\infty }$ are orthogonal with respect to the weight function ${\omega }_{a,b}^{\lambda ,\alpha ,\beta }\left(x\right)={\lambda }^{-1}{(b^{1/\lambda }-x^{1/\lambda })}^{\alpha }{(x^{1/\lambda }-a^{1/\lambda })}^{\beta }{x}^{1/\lambda -1}$ over the interval $[a,b]$, where $a\ge 0.$

Let $X_h:=\{t_0,\cdots ,t_N:0=t_0<t_1<\cdots <t_N=1\}$ be a uniform grid with nodes $t_n=nh$, where $n=0,1,\cdots ,N$ and $h=1/N.$ By employing the graded parameter $\lambda ,$ we can divide the interval $[0,1]$ into N subintervals; that is, ${\displaystyle [t_0^{\lambda },t_1^{\lambda }],\cdots ,[t_{N-1}^{\lambda },t_N^{\lambda }],}$ where ${\displaystyle t_n^{\lambda }={\left(\frac{n}{N}\right)}^{\lambda }.}$ Furthermore, we define the piecewise fractional projection operator ${\mathcal{P}}^{\lambda ,m,N}$ with respect to the grid $\{t_0^{\lambda },t_1^{\lambda },\cdots ,t_N^{\lambda }\}$ as

$$\left({\mathcal{P}}^{\lambda ,m,N}f\right)\left(t\right)=({\mathcal{P}}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{m,\lambda }f)\left(t\right),\forall t\in (t_n^{\lambda },t_{n+1}^{\lambda }],$$

where ${\mathcal{P}}_{a,b}^{m,\lambda }$ denotes the projection operator, satisfying

$$<f,{\psi }_{k,a,b}^{\lambda ,\alpha ,\beta }{>}_{a,b}^{\lambda ,\alpha ,\beta }=<{\mathcal{P}}_{a,b}^{m,\lambda }f,{\psi }_{k,a,b}^{\lambda ,\alpha ,\beta }{>}_{a,b}^{\lambda ,\alpha ,\beta },k=0,1,\cdots ,m,$$

over the interval $(a,b]$ and where ${\mathcal{P}}_{a,b}^{m,\lambda }f\equiv 0$ on $[0,1]\setminus (a,b]$. Here, the test function ${\psi }_{k,a,b}^{\lambda ,\alpha ,\beta }\left(t\right)$ is defined by

$$\begin{array}{c}\hfill {\psi }_{k,a,b}^{\lambda ,\alpha ,\beta }\left(t\right)=\left\{\begin{array}{cc}{P}_{k,a,b}^{\lambda ,\alpha ,\beta }\left(t\right),\hfill & t\in (a,b],\hfill \\ 0,\hfill & \mathrm{others},\hfill \end{array}\right.\end{array}$$

In addition, the inner product $<·,·{>}_{a,b}^{\lambda ,\alpha ,\beta }$ is defined by

$$<f,g{>}_{a,b}^{\lambda ,\alpha ,\beta }={\int }_a^b{\lambda }^{-1}{(b^{1/\lambda }-t^{1/\lambda })}^{\alpha }{(t^{1/\lambda }-a^{1/\lambda })}^{\beta }{t}^{1/\lambda -1}f\left(t\right)g\left(t\right)dt.$$

By employing the transformation (see [22])

$$y\left(t\right)=y_0+\widehat{y}\left(t\right){\mathrm{e}}^{{\int }_0^{t}p\left(s\right)ds},$$

Then, we have

$$\begin{array}{c}\hfill {\widehat{y}}^{\prime }\left(t\right)=\widehat{q}\left(t\right)+{\int }_0^t\widehat{K}(t,s)\widehat{y}\left(s\right)ds\end{array}$$

(4)

with

$$\begin{array}{cc}\hfill \widehat{y}\left(0\right)=& 0,\widehat{K}(t,s)={(t-s)}^{\mu }K(t,s){\mathrm{e}}^{-{\int }_s^{t}p\left(\tau \right)d\tau },\hfill \\ \hfill \widehat{q}\left(t\right)=& \left(y_{0}p\left(t\right)+q\left(t\right)+y_0+{\int }_0^t{(t-s)}^{\mu }K(t,s)ds\right){\mathrm{e}}^{-{\int }_0^{t}p\left(s\right)ds}.\hfill \end{array}$$

Integration of both sides of Equation (4) results in

$${\int }_0^t{y}^{\prime }\left(s\right)ds={\int }_0^t\widehat{q}\left(s\right)ds+{\int }_0^t{\int }_0^s\widehat{K}(s,v)y\left(v\right)dvds.$$

With the help of Dirichlet’s formula

$${\int }_0^t{\int }_0^s\mathsf{\Phi }(s,v)dvds={\int }_0^t{\int }_v^t\mathsf{\Phi }(s,v)dvds,0\le v\le s\le t,$$

Then, we obtain

$$\begin{array}{cc}\hfill \widehat{y}\left(t\right)=& {\int }_0^t\widehat{q}\left(s\right)ds+{\int }_0^t\left({\int }_s^t\widehat{K}(v,s)dv\right)\widehat{y}\left(s\right)ds\hfill \\ \hfill =& \tilde{q}\left(t\right)+\left(\tilde{\mathcal{K}}\widehat{y}\right)\left(t\right),\hfill \end{array}$$

(5)

or

$$\begin{array}{c}\hfill y\left(t\right)=\tilde{q}\left(t\right)+\left(\tilde{\mathcal{K}}y\right)\left(t\right),\end{array}$$

(6)

for simplicity, where

$$\tilde{q}\left(t\right):={\int }_0^{t}q\left(s\right)ds\mathrm{and}\left(\tilde{\mathcal{K}}y\right)\left(t\right):={\int }_0^t\tilde{K}(t,s)y\left(s\right)ds={\int }_0^t\left({\int }_s^t\widehat{K}(v,s)dv\right)y\left(s\right)ds.$$

Recalling the operator ${\mathcal{P}}^{\lambda ,m,N},$ we obtain the piecewise fractional Galerkin solution $y_N\left(t\right)$ with

$$\begin{array}{c}\hfill {\mathcal{P}}^{\lambda ,m,N}{y}_N={\mathcal{P}}^{\lambda ,m,N}\tilde{q}+{\mathcal{P}}^{\lambda ,m,N}\tilde{\mathcal{K}}{y}_N,\end{array}$$

(7)

or equivalently

$$\begin{array}{c}\hfill <y_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }=<\tilde{q},{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }+<\tilde{\mathcal{K}}{y}_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }.\end{array}$$

(8)

In fact, we can represent the approximate solution $y_N\left(t\right)$ as

$$\begin{array}{c}\hfill y_N\left(t\right)=\sum _{n=0}^{N-1}\sum _{j=0}^m{c}_j^n{\varphi }_{j,n}^{\alpha ,\beta ,\lambda }\left(t\right),t\in (0,1],\end{array}$$

(9)

with the local basis function

$$\begin{array}{c}\hfill {\varphi }_{j,n}^{\alpha ,\beta ,\lambda }\left(t\right)=\left\{\begin{array}{cc}{\displaystyle P_{j,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right),}\hfill & t\in (t_n^{\lambda },t_{n+1}^{\lambda }],\hfill \\ 0,\hfill & \mathrm{others}.\hfill \end{array}\right.\end{array}$$

We define the notation

$$\begin{array}{c}\hfill {\phi }_{n,\lambda }\left(x\right)={\left(t_n+\frac{t_{n+1}-t_n}{2}(1+x)\right)}^{\lambda }.\end{array}$$

By employing variable transformations

$$t=s^{\lambda }\mathrm{and}s=t_n+\frac{t_{n+1}-t_n}{2}(1+x),$$

Then, we can compute

$$\begin{array}{cc}\hfill & <y_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\lambda }^{-1}{t}^{1/\lambda -1}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{y}_N\left(t\right)P_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)dt\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\lambda }^{-1}{t}^{1/\lambda -1}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)\sum _{j=0}^m{c}_j^n{P}_{j,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)dt\hfill \\ \hfill =& {\int }_{t_n}^{t_{n+1}}{(t_{n+1}-s)}^{\alpha }{(s-t_n)}^{\beta }{\lambda }^{-1}{s}^{1-\lambda }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(s^{\lambda }\right)\sum _{j=0}^m{c}_j^n{P}_{j,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(s^{\lambda }\right)\lambda s^{\lambda -1}ds\hfill \\ \hfill =& \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right)\sum _{j=0}^m{c}_j^n{P}_j^{\alpha ,\beta }\left(x\right)dx,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & <\tilde{q},{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\lambda }^{-1}{t}^{1/\lambda -1}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)\tilde{q}\left(t\right)dt\hfill \\ \hfill =& {\int }_{t_n}^{t_{n+1}}{(t_{n+1}-s)}^{\alpha }{(s-t_n)}^{\beta }{\lambda }^{-1}{s}^{1-\lambda }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(s^{\lambda }\right)\tilde{q}\left(s^{\lambda }\right)\lambda s^{\lambda -1}ds\hfill \\ \hfill =& \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right)\tilde{q}\left({\phi }_{n,\lambda }\left(x\right)\right)dx,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & <\tilde{\mathcal{K}}{y}_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\lambda }^{-1}{t}^{1/\lambda -1}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right){\int }_0^t\tilde{K}(t,v)y_N\left(v\right)dvdt\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\lambda }^{-1}{t}^{1/\lambda -1}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right){\int }_{t_n^{\lambda }}^t\tilde{K}(t,v)y_N\left(v\right)dvdt\hfill \\ \hfill & +{\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\lambda }^{-1}{t}^{1/\lambda -1}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)\sum _{l=0}^{n-1}{\int }_{t_l^{\lambda }}^{t_{l+1}^{\lambda }}\tilde{K}(t,v)y_N\left(v\right)dvdt\hfill \\ \hfill =& {\int }_{t_n}^{t_{n+1}}{\left(t_{n+1}-s\right)}^{\alpha }{\left(s-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(s^{\lambda }\right){\int }_{t_n^{\lambda }}^{s^{\lambda }}\tilde{K}(s^{\lambda },v)y_N\left(v\right)dvds\hfill \\ \hfill & +{\int }_{t_n}^{t_{n+1}}{\left(t_{n+1}-s\right)}^{\alpha }{\left(s-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(s^{\lambda }\right)\sum _{l=0}^{n-1}{\int }_{t_l^{\lambda }}^{t_{l+1}^{\lambda }}\tilde{K}(s^{\lambda },v)y_N\left(v\right)dvds\hfill \\ \hfill =& \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{\left(1-x\right)}^{\alpha }{\left(x+1\right)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right){\int }_{t_n^{\lambda }}^{{\phi }_{n,\lambda }\left(x\right)}\tilde{K}({\phi }_{n,\lambda }\left(x\right),v)y_N\left(v\right)dvdx\hfill \\ \hfill & +\frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{\left(1-x\right)}^{\alpha }{\left(x+1\right)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right)\sum _{l=0}^{n-1}{\int }_{t_l^{\lambda }}^{t_{l+1}^{\lambda }}\tilde{K}({\phi }_{n,\lambda }\left(x\right),v)y_N\left(v\right)dvdx.\hfill \end{array}$$

Hence, the moment integral $<y_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }$ can be computed in a closed form, and other moments, such as ${\displaystyle <\tilde{q},{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }}$ and $<\tilde{\mathcal{K}}{y}_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta },$ can be efficiently evaluated by Gauss-type quadrature rules in Chebfun (see [23]). Once the linear system in Equation (8) is solved step by step, we will obtain the piecewise Galerkin solution $y_N\left(t\right)$ immediately.

3. Solvability and Convergence Property of PFG

In this section, we will study the solvability and convergence property of the Galerkin approximation $y_N\left(t\right)$, defined in Equation (8). We begin with the approximation results of truncated spectral expansions. Let us restrict the considered function in the current paper to the following definition:

Definition 2.

Any function $f\left(x\right)$ is said to be in $\mathcal{A}\left(\right[a,b],\gamma )$ with $\gamma >-1$ if it satisfies the following conditions:

• $f\left(x\right)$ has the form of ${(x-a)}^{\gamma }{g}_1\left(x\right)+g_2\left(x\right)$, with γ being a positive real number;
• Both $g_1\left(x\right)$ and $g_2\left(x\right)$ are analytic in a sufficiently large domain containing the interval $[a,b].$

If we assume $f\left(x\right)$ belongs to $\mathcal{A}\left(\right[-1,1],\gamma )$, then its orthogonal polynomial expansion with respect to the weight function $\omega \left(x\right)={(1-x)}^{\alpha }{(x+1)}^{\beta }$ is defined as

$$f\left(x\right)=\sum _{n=0}^{\infty }{c}_n(\alpha ,\beta )P_n^{\alpha ,\beta }\left(x\right),$$

where

$$c_n(\alpha ,\beta )=\frac{1}{{\sigma }_n}{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }f\left(x\right)P_n^{\alpha ,\beta }\left(x\right)dx,$$

and

$${\sigma }_n={\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_n^{\alpha ,\beta }\left(x\right)P_n^{\alpha ,\beta }\left(x\right)dx.$$

Analysis of the decay rate of such coefficients helps with studying the error bound derived from the truncation of the Jacobi series (see [24,25]). In fact, according to Rodrigues’ formula ([19], p. 72), we have

$$\begin{array}{c}\hfill {(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_n^{\alpha ,\beta }\left(x\right)=\frac{{(-1)}^k}{2^k}\prod _{j=0}^{k-1}\frac{1}{n-j}\frac{d^k}{d{x}^k}\left({(1-x)}^{\alpha +k}{(1+x)}^{\beta +k}{P}_{n-k}^{\alpha +k,\beta +k}\left(x\right)\right).\end{array}$$

A direct integration by parts results in

$$\begin{array}{c}\hfill c_n(\alpha ,\beta )=\frac{{\displaystyle {\int }_{-1}^1{(1-x)}^{\alpha +k}{(1+x)}^{\beta +k}{f}^{\left(k\right)}\left(x\right)P_{n-k}^{\alpha +k,\beta +k}\left(x\right)dx}}{{\displaystyle {\sigma }_n{2}^k\prod _{j=0}^{k-1}(n-j)}}.\end{array}$$

(10)

Now, let us consider the Jacobi expansion on an arbitrary interval $[a,b]$ with $a\ge 0.$ Suppose that $f\left(x\right)={(x-a)}^{\gamma }{g}_1\left(x\right)+g_2\left(x\right)$ is defined on $[a,b],$ where $g_1\left(x\right)$ and $g_2\left(x\right)$ are analytic in a sufficiently large domain containing the interval $[a,b]$. Then, the coefficients of the Jacobi expansion of $f\left(x\right)$ can be computed by the following for $k=0,1,\cdots ,n$:

$$\begin{array}{cc}\hfill c_n(\alpha ,\beta )=& \frac{1}{{\sigma }_n^{a,b}}{\left(\frac{b-a}{2}\right)}^{1+\alpha +\beta }{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }f\left(\frac{b-a}{2}x+\frac{a+b}{2}\right)P_n^{\alpha ,\beta }\left(x\right)dx\hfill \\ \hfill =& {\left(\frac{b-a}{2}\right)}^{k+1+\alpha +\beta }\frac{{\displaystyle {\int }_{-1}^1{(1-x)}^{\alpha +k}{(1+x)}^{\beta +k}{f}^{\left(k\right)}\left(\frac{b-a}{2}x+\frac{a+b}{2}\right)P_{n-k}^{\alpha +k,\beta +k}\left(x\right)dx}}{{\displaystyle {\sigma }_n^{a,b}{2}^k\prod _{j=0}^{k-1}(n-j)}},\hfill \end{array}$$

(11)

where

$${\sigma }_n^{a,b}={\int }_a^b{(b-x)}^{\alpha }{(x-a)}^{\beta }{P}_{n,a,b}^{1,\alpha ,\beta }\left(x\right)P_{n,a,b}^{1,\alpha ,\beta }\left(x\right)dx$$

Furthermore, by letting $b-a\to 0,$ we have the following by noting for Equation (11) that $n=0,1,2,\cdots$:

$$\begin{array}{c}\hfill c_n(\alpha ,\beta )=\left\{\begin{array}{cc}O\left({(b-a)}^{\gamma }\right),\hfill & \gamma <n,\hfill \\ O\left({(b-a)}^n\right),\hfill & \gamma \ge n.\hfill \end{array}\right.\end{array}$$

(12)

Let ${\mathcal{P}}_{a,b}^n$ denote the projection operator which satisfies

$$<y,P_{k,a,b}^{1,\alpha ,\beta }{>}_{a,b}^{\alpha ,\beta }=<{\mathcal{P}}_{a,b}^{n}y,P_{k,a,b}^{1,\alpha ,\beta }{>}_{a,b}^{\alpha ,\beta },k=0,1,\cdots ,n.$$

Here, the inner product $<·,·{>}_{a,b}^{\alpha ,\beta }$ is defined by

$$<f,g{>}_{a,b}^{\alpha ,\beta }={\int }_a^b{(b-x)}^{\alpha }{(x-a)}^{\beta }f\left(x\right)g\left(x\right)dx.$$

Since

$$(I-{\mathcal{P}}_{a,b}^n)f=\sum _{j=n+1}^{\infty }{c}_j(\alpha ,\beta )P_{j,a,b}^{1,\alpha ,\beta }\left(x\right),$$

Then, we obtain the estimate for the truncated Jacobi approximation; that is, we have the following:

Lemma 2.

Suppose $f\left(x\right)$ belongs to $\mathcal{A}\left(\right[a,b],\gamma )$. Then, for the Jacobi expansion, it follows that as $b-a\to 0$, the following holds:

$$\parallel (I-{\mathcal{P}}_{a,b}^n){f\parallel }_{\infty }=\left\{\begin{array}{cc}O\left({(b-a)}^{\gamma }\right),\hfill & \gamma <n+1,\hfill \\ O\left({(b-a)}^{n+1}\right),\hfill & \gamma \ge n+1.\hfill \end{array}\right.$$

The above classical Jacobi polynomial theory provides us powerful tools for analyzing the fractional Jacobi expansions. Suppose that

$$e\left(x\right):=\left((I-{\mathcal{P}}_{a,b}^{n,\lambda })f\right)\left(x\right)=\sum _{j=n+1}^{\infty }{c}_j\left(\lambda \right)P_{j,a,b}^{\lambda ,\alpha ,\beta }\left(x\right).$$

For $j\ge n+1$, we have

$$\begin{array}{cc}\hfill c_j\left(\lambda \right)=& \frac{{\displaystyle {\int }_a^b{\lambda }^{-1}{x}^{1/\lambda -1}{\left(b^{1/\lambda }-x^{1/\lambda }\right)}^{\alpha }{\left(x^{1/\lambda }-a^{1/\lambda }\right)}^{\beta }e\left(x\right)P_{j,a,b}^{\lambda ,\alpha ,\beta }\left(x\right)dx}}{{\displaystyle {\int }_a^b{\lambda }^{-1}{x}^{1/\lambda -1}{\left(b^{1/\lambda }-x^{1/\lambda }\right)}^{\alpha }{\left(x^{1/\lambda }-a^{1/\lambda }\right)}^{\beta }{P}_{j,a,b}^{\lambda ,\alpha ,\beta }\left(x\right)P_{j,a,b}^{\lambda ,\alpha ,\beta }\left(x\right)dx}}\hfill \\ \hfill =& \frac{{\displaystyle {\int }_{a^{1/\lambda }}^{b^{1/\lambda }}{\left(b^{1/\lambda }-t\right)}^{\alpha }{\left(t-a^{1/\lambda }\right)}^{\beta }e\left(t^{\lambda }\right)P_{j,a,b}^{\lambda ,\alpha ,\beta }\left(t^{\lambda }\right)dt}}{{\displaystyle {\int }_{a^{1/\lambda }}^{b^{1/\lambda }}{\left(b^{1/\lambda }-t\right)}^{\alpha }{\left(t-a^{1/\lambda }\right)}^{\beta }{P}_{j,a,b}^{\lambda ,\alpha ,\beta }\left(t^{\lambda }\right)P_{j,a,b}^{\lambda ,\alpha ,\beta }\left(t^{\lambda }\right)dt}}\hfill \\ \hfill =& \frac{{\displaystyle {\int }_{-1}^1{\left(1-s\right)}^{\alpha }{\left(s-1\right)}^{\beta }e\left({\left(\frac{b^{1/\lambda }-a^{1/\lambda }}{2}s+\frac{b^{1/\lambda }+a^{1/\lambda }}{2}\right)}^{\lambda }\right)P_j^{\alpha ,\beta }\left(s\right)ds}}{{\displaystyle {\int }_{-1}^1{\left(1-s\right)}^{\alpha }{\left(s-1\right)}^{\beta }{P}_j^{\alpha ,\beta }\left(s\right)P_j^{\alpha ,\beta }\left(s\right)ds}}\hfill \\ \hfill =& \frac{1}{{\sigma }_j}{\int }_{-1}^1{\left(1-s\right)}^{\alpha }{\left(s-1\right)}^{\beta }f\left({\left(\frac{b^{1/\lambda }-a^{1/\lambda }}{2}s+\frac{b^{1/\lambda }+a^{1/\lambda }}{2}\right)}^{\lambda }\right)P_j^{\alpha ,\beta }\left(s\right)ds.\hfill \end{array}$$

Here, the variable transformations

$$x=t^{\lambda }\mathrm{and}t=\frac{b^{1/\lambda }-a^{1/\lambda }}{2}s+\frac{b^{1/\lambda }+a^{1/\lambda }}{2}$$

are employed. If $f\left(x\right)\in \mathcal{A}\left(\right[a,b],\gamma ),$ then we have ${\displaystyle f\left({\left(\frac{b^{1/\lambda }-a^{1/\lambda }}{2}s+\frac{b^{1/\lambda }+a^{1/\lambda }}{2}\right)}^{\lambda }\right)}$ belonging to $\mathcal{A}([a^{1/\lambda },b^{1/\lambda }],\lambda \gamma ),$ which implies

$$\begin{array}{c}\hfill c_n^{\lambda }(\alpha ,\beta )=\left\{\begin{array}{cc}O\left({(b^{1/\lambda }-a^{1/\lambda })}^{\lambda \gamma }\right),\hfill & \lambda \gamma <n+1,\hfill \\ O\left({(b^{1/\lambda }-a^{1/\lambda })}^{n+1}\right),\hfill & \lambda \gamma \ge n+1.\hfill \end{array}\right.\end{array}$$

(13)

With Lemma 2 in mind, we arrive at the following result:

Corollary 1.

Suppose $f\left(x\right)$ belongs to $\mathcal{A}\left(\right[a,b],\gamma )$, where $a\ge 0.$ Then, for the fractional Jacobi expansion, it follows that as $b^{1/\lambda }-a^{1/\lambda }\to 0$, we have

$$\parallel (I-{\mathcal{P}}_{a,b}^{n,\lambda }){f\parallel }_{\infty }=\left\{\begin{array}{cc}O\left({(b^{1/\lambda }-a^{1/\lambda })}^{\lambda \gamma }\right),\hfill & \lambda \gamma <n+1,\hfill \\ O\left({(b^{1/\lambda }-a^{1/\lambda })}^{n+1}\right),\hfill & \lambda \gamma \ge n+1.\hfill \end{array}\right.$$

The solution to the VIDE in Equation (1) can be analyzed with the help of the basic Volterra theory of the second kind. For sufficiently smooth $p\left(t\right),q\left(t\right)$ and $K(t,s),$ its asymptotic expansion can be obtained by utilizing Picard’s iteration ([2], Theorem 7.1.4). We summarize these theoretical results in the following lemma:

Lemma 3.

Assume $q\in \mathcal{A}\left(\right[0,1],\mu +1)$ and $K\in C^m\left(D\right)$, where $K(t,t)\ne 0$, $p\in C^m\left([0,1]\right)$ and $\mu \in (-1,0).$ Then, the solution y of Equation (1) can be written in the following form:

$$\begin{array}{c}\hfill y\left(t\right)=\sum _{{(j,k)}_{\mu }}{\gamma }_{j,k}\left(\mu \right)t^{j+k(2+\mu )}+Y_{m+1}(t,\mu ),\end{array}$$

(14)

where

$${(j,k)}_{\mu }:=\{(j,k):j,k\in N,j+k(2+\mu )<m+1\}.$$

Moreover, $Y_{m+1}(t,\mu )\in C^{m+1}\left(I\right).$

It is noted that the exact solution $y\left(t\right)$ to the VIDE in Equation (1) belongs to $\mathcal{A}\left(\right[0,1],2+\mu )$ under the assumption of Lemma 3. Now, we arrive at the main theoretical result of the current paper:

Theorem 1.

Assume the following:

• The regularization parameter λ determining the mesh $X_N^{\lambda }$ is given;
• The given functions in Equation (1) satisfy $K(t,s)\in C^{m+1}\left(D\right)n$ $p\in C^m\left([0,1]\right)$ and $q\left(t\right)\in \mathcal{A}\left(\right[0,1],1+\mu ).$

Let $y_N$ denote the piecewise fractional Galerkin solution computed by Equation (8), and denote the error function as $e_N\left(t\right):=y\left(t\right)-y_N\left(t\right).$ Then, the maximum of $e_N\left(t\right)$ over the interval $[0,1]$ satisfies

$$\begin{array}{c}\hfill \parallel e_N{\left(t\right)\parallel }_{\infty }\le B\left\{\begin{array}{c}{N}^{-\lambda (2+\mu )},\lambda (2+\mu )<m+1,\hfill \\ N^{-m-1},\lambda (2+\mu )\ge m+1.\hfill \end{array}\right.\end{array}$$

(15)

Here, the constant B does not depend on m and $N.$

Proof. 

First, we study the existence and uniqueness of the Galerkin approximation $y_N\left(t\right).$ This is equivalent to proving the fact that $y_N\left(t\right)=0$ is the unique solution of the homogenous equation. In fact, the Galerkin solution for the corresponding homogenous version of Equation (1) satisfies

$$\begin{array}{c}\hfill {\mathcal{P}}^{\lambda ,m,N}{y}_N={\mathcal{P}}^{\lambda ,m,N}\tilde{\mathcal{K}}{y}_N.\end{array}$$

(16)

Noting the localization of the operator ${\mathcal{P}}^{\lambda ,m,N},$ we have

$$\begin{array}{c}\hfill <y_N,{\psi }_{k,0,t_1^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{0,t_1^{\lambda }}^{\lambda ,\alpha ,\beta }=<\tilde{\mathcal{K}}{y}_N,{\psi }_{k,0,t_1^{\lambda }}^{\lambda ,\alpha ,\beta }{>}_{0,t_1^{\lambda }}^{\lambda ,\alpha ,\beta }.\end{array}$$

(17)

The coefficient matrix $A_{left}$ on the left-hand side is diagonal and invertible due to the orthogonality of $P_n^{\alpha ,\beta }\left(t\right).$ Furthermore, by letting $h\to 0,$ we have $\parallel A_{left}{\parallel }_{\infty }=O\left(h^{1+\alpha +\beta }\right).$ For the right-hand side, it follows by a direct calculation that

$$\begin{array}{cc}\hfill <\tilde{\mathcal{K}}{y}_N,y_N{>}_{0,t_1^{\lambda }}^{\lambda ,\alpha ,\beta }=& {\int }_0^{t_1^{\lambda }}{(t_1-t^{1/\lambda })}^{\alpha }{t}^{\beta /\lambda }{y}_N\left(t\right){\int }_0^t\tilde{K}(t,v)y_N\left(v\right)dvdt\hfill \\ \hfill \le & \frac{\overline{K}}{\mu +1}{\parallel y_N\parallel }_{\infty }^2{\int }_0^{t_1^{\lambda }}{(t_1-t^{1/\lambda })}^{\alpha }{t}^{\beta /\lambda }{\int }_0^{t}dvdt\hfill \\ \hfill \le & B{h}^{\lambda +1+\alpha +\beta }{\parallel y_N\parallel }_{\infty }^2,\hfill \end{array}$$

where $\overline{K}$ denotes the maximum of $K(t,s)$ in the domain $\left\{\right(t,s):0\le s\le t\le 1\}.$ This implies that the maximum of the coefficient matrix on the right-hand side is $O\left(h^{\lambda +1+\alpha +\beta }\right).$ Hence, it follows that $y_N\left(t\right)\equiv 0$ for any $t\in [0,t_1^{\lambda }].$ Noting that

$$\begin{array}{c}\hfill {\mathcal{P}}^{\lambda ,m,N}{y}_N={\mathcal{P}}^{\lambda ,m,N}\left(\sum _{j=0}^{n-1}{\int }_{t_j^{\lambda }}^{t_{j+1}^{\lambda }}\tilde{K}(t,s)y_N\left(s\right)ds\right)+{\mathcal{P}}^{\lambda ,m,N}\left({\int }_{t_n^{\lambda }}^t\tilde{K}(t,s)y_N\left(s\right)ds\right)\end{array}$$

Then, we can continue $y_N\left(t\right)\equiv 0$ to $[0,1]$ by the induction, which implies that Equation (8) can be uniquely solved.

Secondly, let us study the convergence property of the fractional Galerkin solution $y_N\left(t\right)$. Noting that

$$\begin{array}{c}\hfill y_N={\mathcal{P}}^{\lambda ,m,N}\tilde{q}+{\mathcal{P}}^{\lambda ,m,N}\tilde{\mathcal{K}}{y}_N,y=\tilde{q}+\tilde{\mathcal{K}}y,\end{array}$$

Then, we can compute

$$\begin{array}{cc}\hfill e_N=& ({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{q}+{\mathcal{P}}^{\lambda ,m,N}\tilde{\mathcal{K}}{y}_N-\tilde{\mathcal{K}}y\hfill \\ \hfill =& ({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{q}+{\mathcal{P}}^{\lambda ,m,N}\tilde{\mathcal{K}}{e}_N+({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{\mathcal{K}}y.\hfill \end{array}$$

Let ${ϵ}_n$ denote the maximum of $e_N\left(t\right)$ on $(t_n^{\lambda },t_{n+1}^{\lambda }].$ Then, for any $t\in (t_n^{\lambda },t_{n+1}^{\lambda }],$ it follows that

$$\begin{array}{c}\hfill {ϵ}_n\le \parallel ({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{q}{\parallel }_{\infty }+\parallel {\mathcal{P}}^{\lambda ,m,N}\tilde{\mathcal{K}}{e}_N{\parallel }_{\infty }+{\parallel ({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{\mathcal{K}}y\parallel }_{\infty }.\end{array}$$

(18)

The remaining work is to estimate the above equation term by term. In fact, a direct calculation leads to

$$\begin{array}{cc}\hfill \parallel {\mathcal{P}}^{\lambda ,m,N}\tilde{\mathcal{K}}{e}_N{\parallel }_{\infty }\le & \overline{K}{\int }_0^t\left({\int }_s^t{(t-\tau )}^{\mu }d\tau \right)\left|e_N\left(s\right)\right|ds\hfill \\ \hfill =& \frac{\overline{K}}{\mu +1}{\int }_0^t{(t-s)}^{\mu +1}\left|e_N\left(s\right)\right|ds\hfill \\ \hfill =& \frac{\overline{K}}{\mu +1}\left(\sum _{j=0}^{n-1}{\int }_{t_j^{\lambda }}^{t_{j+1}^{\lambda }}{(t-s)}^{\mu +1}|e_N\left(s\right)|ds+{\int }_{t_n^{\lambda }}^t{(t-s)}^{\mu +1}\left|e_N\left(s\right)\right|ds\right)\hfill \\ \hfill \le & \frac{\overline{K}}{\mu +1}\sum _{j=0}^n{l}_j{ϵ}_j,\hfill \end{array}$$

where

$$l_j=\left\{\begin{array}{cc}{\displaystyle {\int }_{t_j^{\lambda }}^{t_{j+1}^{\lambda }}{(t-s)}^{\mu +1}ds,}\hfill & j=0,1,\cdots ,n-1,\hfill \\ {\displaystyle {\int }_{t_n^{\lambda }}^t{(t-s)}^{\mu +1}ds,}\hfill & j=n.\hfill \end{array}\right.$$

Then, we estimate $l_j$ in two cases: $j=n$ (Case I) and $0\le j\le n-1$ (Case II).

For Case I ($j=n$), we obtain

$$\begin{array}{cc}\hfill {\int }_{t_n^{\lambda }}^t{(t-s)}^{\mu +1}ds\le & {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{(t_{n+1}^{\lambda }-s)}^{\mu +1}ds\hfill \\ \hfill =& \frac{1}{2+\mu }{(t_{n+1}^{\lambda }-t_n^{\lambda })}^{2+\mu }\hfill \\ \hfill \le & \frac{1}{2+\mu }{N}^{-\lambda (2+\mu )}{(n+{\eta }_1)}^{(\lambda -1)(2+\mu )}\hfill \\ \hfill \le & \frac{1}{2+\mu }{N}^{-\lambda (2+\mu )}{(n+1)}^{(\lambda -1)(2+\mu )},\hfill \end{array}$$

(19)

with ${\eta }_1\in (0,1).$

For Case II ($0\le j\le n-1$), a direct calculation leads to

$$\begin{array}{cc}\hfill {\int }_{t_j^{\lambda }}^{t_{j+1}^{\lambda }}{(t-s)}^{\mu +1}ds\le & {\int }_{t_j^{\lambda }}^{t_{j+1}^{\lambda }}{(t_{n+1}^{\lambda }-s)}^{\mu +1}ds\hfill \\ \hfill \le & \left({\left(\frac{j+1}{N}\right)}^{\lambda }-{\left(\frac{j}{N}\right)}^{\lambda }\right){\left({\left(\frac{n+1}{N}\right)}^{\lambda }-{\left(\frac{j+1}{N}\right)}^{\lambda }\right)}^{\mu +1}\hfill \\ \hfill =& \lambda {\left(\frac{j+{\eta }_2}{N}\right)}^{\lambda -1}{\lambda }^{\mu +1}{\left(\frac{n-j}{N}\right)}^{\mu +1}{\left(\frac{j+1+{\eta }_3}{N}\right)}^{(\lambda -1)(\mu +1)}\hfill \\ \hfill \le & {\lambda }^{\mu +2}{N}^{-\lambda (\mu +2)}{(j+1)}^{\lambda -1}{j}^{(\lambda -1)(\mu +1)}{(n-j)}^{\mu +1}\hfill \\ \hfill \le & {\lambda }^{\mu +2}{N}^{-\lambda (\mu +2)}{(n-j)}^{\mu +1}{(j+1)}^{(\lambda -1)(\mu +2)},\hfill \end{array}$$

(20)

with ${\eta }_2,{\eta }_3\in (0,1).$

To sum up, it follows that

$$\begin{array}{cc}\hfill \sum _{j=0}^n{l}_j\le & B{N}^{-\gamma (\mu +2)}(n^{\alpha }+{(n-1)}^{\mu +1}{2}^{(\gamma -1)(\mu +2)}+\cdots \hfill \\ \hfill & +2^{\mu +1}{(n-1)}^{(\gamma -1)(\mu +2)}+1^{\mu +1}{n}^{(\gamma -1)(\mu +2)}+{(n+1)}^{(\gamma -1)(\mu +2)})\hfill \\ \hfill \le & B{(n+1)}^{-\gamma (\mu +2)}(n^{\mu +1}+{(n-1)}^{\mu +1}{2}^{(\gamma -1)(\mu +2)}+\cdots \hfill \\ \hfill & +2^{\mu +1}{(n-1)}^{(\gamma -1)(\mu +2)}+1^{\mu +1}{n}^{(\gamma -1)(\mu +2)}+{(n+1)}^{(\gamma -1)(\mu +2)})\hfill \\ \hfill =& B\left({(n+1)}^{-(\mu +2)}+\frac{1}{n+1}\sum _{j=1}^n{\left(1-\frac{j}{n+1}\right)}^{\mu +1}{\left(\frac{j}{n+1}\right)}^{(\gamma -1)(\mu +2)}\right)\hfill \\ \hfill \le & B\left(1+{\int }_0^1{(1-x)}^{\mu +1}{x}^{(\gamma -1)(\mu +2)}dx\right)\hfill \\ \hfill =& B\left(1+{\mathrm{B}}_{\mathrm{ETA}}\left(\right(\gamma -1\left)\right(\mu +2)+1,\mu +2)\right).\hfill \end{array}$$

(21)

Here, ${\mathrm{B}}_{\mathrm{ETA}}(·,·)$ denotes the beta function ([26], p. 142).

Noting that $({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{q}$ and $({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{\mathcal{K}}y$ belong to $\mathcal{A}\left(\right[0,1],1+\mu ),$ we obtain their error bounds by Lemma 2:

$$\begin{array}{c}\hfill \parallel ({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{q}{\parallel }_{\infty }:={ϵ}_q\le B\left\{\begin{array}{c}{N}^{-\lambda (2+\mu )},\lambda (2+\mu )<m+1,\hfill \\ N^{-m-1},\lambda (2+\mu )\ge m+1.\hfill \end{array}\right.\end{array}$$

(22)

$$\begin{array}{c}\hfill \parallel ({\mathcal{P}}^{\lambda ,m,N}-I)\tilde{\mathcal{K}}{y\parallel }_{\infty }:={ϵ}_{Ky}\le B\left\{\begin{array}{c}{N}^{-\lambda (2+\mu )},\lambda (2+\mu )<m+1,\hfill \\ N^{-m-1},\lambda (2+\mu )\ge m+1.\hfill \end{array}\right.\end{array}$$

(23)

Here, the constant B does not depend on m and $N.$

On the other hand, since the inequality

$$\begin{array}{cc}\hfill {ϵ}_n\le & {ϵ}_q+{ϵ}_{Ky}+\frac{\overline{K}}{\mu +1}\sum _{j=0}^n{l}_j{ϵ}_j\hfill \\ \hfill \le & {ϵ}_q+{ϵ}_{Ky}+\frac{\overline{K}{l}_nϵ}{\mu +1}+\frac{\overline{K}}{\mu +1}\sum _{j=0}^{n-1}{l}_j{ϵ}_j\hfill \end{array}$$

holds, we have the following with the help of Grönwall’s inequality

$$\begin{array}{c}\hfill {ϵ}_n\le \left({ϵ}_q+{ϵ}_{Ky}+\frac{\overline{K}{l}_{n-1}{ϵ}_n}{\mu +1}\right){\mathrm{e}}^{\frac{\overline{K}}{\mu +1}{\sum }_{j=0}^{n-2}{l}_j}\end{array}$$

(24)

Combining Equations (21)–(24) gives

$$\begin{array}{c}\hfill ϵ\le B\left\{\begin{array}{c}{N}^{-\lambda (2+\mu )},\lambda (2+\mu )<m+1,\hfill \\ N^{-m-1},\lambda (2+\mu )\ge m+1.\hfill \end{array}\right.\end{array}$$

(25)

This completes the proof. □

4. Numerical Experiments

In this section, we apply PFG to the weakly singular problem of the form in Equation (1) to demonstrate the effectiveness of the piecewise discretization. The singular parameter $\mu$ is chosen randomly to verify the estimate in Theorem 1.

Consider the weakly singular VIDE

$$\begin{array}{c}\hfill y^{\prime }\left(t\right)=q_{\mu }\left(t\right)+{\int }_0^t{(t-s)}^{\mu }y\left(s\right)ds,t\in [0,1].\end{array}$$

(26)

Here, $q_{\mu }\left(t\right)=(2+\mu )t^{1+\mu }-{\mathrm{B}}_{\mathrm{ETA}}(1+\mu ,3+\mu )t^{2\mu +3}$, and the exact solution is $y\left(t\right)=t^{2+\mu }.$

By letting

$${\mathbf{e}}_N:={\left(e_N\left(t_0^{\lambda }\right),e_N\left(t_0^{\lambda }+\frac{h_0}{10}\right),\cdots ,e_N\left(t_1^{\lambda }\right),\cdots ,e_N\left(t_N^{\lambda }\right)\right)}^T,h_j=t_{j+1}^{\lambda }-t_j^{\lambda },$$

Then, we compute the maximum of ${\mathbf{e}}_N$ and show the numerical results in

Table 1. Absolute errors and convergence rates of PFG for Equation (26) when $\alpha =\beta =0,$ $\mu =-0.2$ and $\lambda =1$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
$N=4$	$2.79\times {10}^{-4}$	–	$1.78\times {10}^{-5}$	–
$N=8$	$8.00\times {10}^{-5}$	1.8023	$5.10\times {10}^{-6}$	1.8004
$N=12$	$3.86\times {10}^{-5}$	1.8008	$2.46\times {10}^{-6}$	1.8002
$N=16$	$2.30\times {10}^{-5}$	1.8004	$1.47\times {10}^{-6}$	1.8001
$N=20$	$1.54\times {10}^{-5}$	1.8003	$9.81\times {10}^{-7}$	1.8001
Refer Order		1.8000		1.8000

,

Table 2. Absolute errors and convergence rates of PFG for Equation (26) when $\alpha =\beta =0,$ $\mu =-0.2$ and $\lambda =2$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
$N=4$	$8.00\times {10}^{-4}$	–	$9.46\times {10}^{-7}$	–
$N=8$	$1.04\times {10}^{-4}$	2.9485	$7.80\times {10}^{-8}$	3.5999
$N=12$	$3.15\times {10}^{-5}$	2.9081	$1.81\times {10}^{-8}$	3.6000
$N=16$	$1.37\times {10}^{-5}$	2.9383	$6.44\times {10}^{-9}$	3.6000
$N=20$	$7.08\times {10}^{-6}$	2.9534	$2.88\times {10}^{-9}$	3.6000
Refer Order		3.0000		3.6000

and

Table 3. Absolute errors and convergence rates of PFG for Equation (26) when $\alpha =\beta =0,$ $\mu =-0.2$ and $\lambda =3$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
$N=4$	$3.48\times {10}^{-3}$	–	$3.47\times {10}^{-6}$	–
$N=8$	$5.03\times {10}^{-4}$	2.7899	$1.11\times {10}^{-7}$	4.9645
$N=12$	$1.55\times {10}^{-4}$	2.9047	$1.45\times {10}^{-8}$	5.0125
$N=16$	$6.73\times {10}^{-5}$	2.8949	$3.44\times {10}^{-9}$	5.0173
$N=20$	$3.55\times {10}^{-5}$	2.8677	$1.12\times {10}^{-9}$	5.0169
Refer Order		3.0000		5.0000

. Here, we choose $\alpha =\beta =0.$ The “error” in the tables represents the maximum of ${\mathbf{e}}_N,$ and the “order” is computed by

$$\mathrm{order}=\frac{log{\mathrm{error}}_{\mathrm{current}}-log{\mathrm{error}}_{\mathrm{previous}}}{log{N}_{\mathrm{current}}-log{N}_{\mathrm{previous}}}.$$

That aside, “Refer Order” denotes the theoretical estimates given in Theorem 1. The results shown in Table 1, Table 2 and Table 3 illustrate that the theoretical convergence rate coincided with the numerical experiments well. In addition, we found that increasing the regularization parameter $\lambda$ promoted the accuracy of the numerical solutions from the comparison between the computed results in Table 1, Table 2 and Table 3.

Next, we compare the proposed PFG with the classical piecewise collocation methods with Gauss–Legendre points in the case of $\mu =\alpha =\beta =-0.5.$ The graded parameters for the collocation methods are denoted as $r.$ In

Table 4. Absolute errors and convergence rates of PFG for Equation (26) when $\alpha =\beta =-0.5,$ $\mu =-0.5$ and $\lambda =1$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
$N=4$	$1.11\times {10}^{-3}$	–	$1.52\times {10}^{-4}$	–
$N=8$	$3.89\times {10}^{-4}$	1.5116	$5.36\times {10}^{-5}$	1.5020
$N=12$	$2.11\times {10}^{-4}$	1.5049	$2.92\times {10}^{-5}$	1.5008
$N=16$	$1.37\times {10}^{-4}$	1.5029	$1.89\times {10}^{-5}$	1.5005
$N=20$	$9.81\times {10}^{-5}$	1.5020	$1.36\times {10}^{-5}$	1.5003
Refer Order		1.5000		1.5000

,

Table 5. Absolute errors and convergence rates of collocation methods on graded meshes for Equation (26) when $\mu =-0.5$ and $r=1$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
4	$4.34\times {10}^{-3}$	–	$6.58\times {10}^{-4}$	–
8	$1.49\times {10}^{-3}$	$1.5472$	$2.31\times {10}^{-4}$	$1.5106$
12	$8.00\times {10}^{-4}$	$1.5260$	$1.25\times {10}^{-4}$	$1.5060$
16	$5.17\times {10}^{-4}$	$1.5181$	$8.13\times {10}^{-5}$	$1.5042$
20	$3.69\times {10}^{-4}$	$1.5139$	$5.82\times {10}^{-5}$	$1.5033$
Refer Order		1.5000		1.5000

,

Table 6. Absolute errors and convergence rates of PFG for Equation (26) when $\alpha =\beta =-0.5,$ $\mu =-0.5$ and $\lambda =2$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
$N=4$	$4.89\times {10}^{-4}$	–	$6.07\times {10}^{-14}$	–
$N=8$	$6.10\times {10}^{-5}$	3.0022	$6.31\times {10}^{-14}$	–
$N=12$	$1.81\times {10}^{-5}$	3.0004	$6.39\times {10}^{-14}$	–
$N=16$	$7.63\times {10}^{-6}$	3.0001	$6.93\times {10}^{-14}$	–
$N=20$	$3.91\times {10}^{-6}$	3.0001	$7.57\times {10}^{-14}$	–
Refer Order		3.0000		3.0000

,

Table 7. Absolute errors and convergence rates of collocation methods on graded meshes for Equation (26) when $\mu =-0.5$ and $r=2$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
4	$8.86\times {10}^{-4}$	–	$8.22\times {10}^{-5}$	–
8	$1.13\times {10}^{-4}$	$2.9670$	$1.02\times {10}^{-5}$	$3.0042$
12	$3.38\times {10}^{-5}$	$2.9844$	$3.03\times {10}^{-6}$	$3.0013$
16	$1.43\times {10}^{-5}$	$2.9893$	$1.28\times {10}^{-6}$	$3.0007$
20	$7.34\times {10}^{-6}$	$2.9916$	$6.55\times {10}^{-7}$	$3.0004$
Refer Order		2.5000		3.0000

,

Table 8. Absolute errors and convergence rates of PFG for Equation (26) when $\alpha =\beta =-0.5,$ $\mu =-0.5$ and $\lambda =3$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
$N=4$	$2.63\times {10}^{-3}$	–	$1.33\times {10}^{-6}$	–
$N=8$	$3.68\times {10}^{-4}$	2.8360	$5.88\times {10}^{-8}$	4.5001
$N=12$	$1.12\times {10}^{-4}$	2.9285	$9.50\times {10}^{-9}$	4.5000
$N=16$	$4.80\times {10}^{-5}$	2.9550	$2.60\times {10}^{-9}$	4.5000
$N=20$	$2.48\times {10}^{-5}$	2.9675	$9.53\times {10}^{-10}$	4.5000
Refer Order		3.0000		4.5000

and

Table 9. Absolute errors and convergence rates of collocation methods on graded meshes for Equation (26) when $\mu =-0.5$ and $r=3$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
4	$1.01\times {10}^{-3}$	–	$1.56\times {10}^{-5}$	–
8	$8.07\times {10}^{-5}$	$3.6519$	$6.90\times {10}^{-7}$	$4.5014$
12	$1.74\times {10}^{-5}$	$3.7848$	$1.11\times {10}^{-7}$	$4.5008$
16	$5.76\times {10}^{-6}$	$3.8413$	$3.05\times {10}^{-8}$	$4.5002$
20	$2.43\times {10}^{-6}$	$3.8750$	$1.12\times {10}^{-8}$	$4.5000$
Refer Order		2.5000		4.5000

, we display the computed results. It can be seen that the theoretical estimate for the convergence rates of PFG was optimal in this case. In the case of $\lambda =2,$ PFG with $m=4$ resulted in numerical solutions which were able to achieve machine precision in Matlab (see the right hand side of Table 6), while the collocation method could not work out similar numerical approximations. We also found that although the theoretical convergence rates of the PFG and collocation methods were the same in the case of $m=4,$ the PFG provided more accurate approximations than the collocation method.

5. Final Remark

We studied the piecewise fractional Galerkin method for the weakly singular VIDE in Equation (1). The existence, uniqueness and convergence property of the piecewise fractional Galerkin solution were analyzed in detail. The theoretical and numerical results showed that the new Galerkin method is efficient at solving weakly singular problems. This is partly due to the approximate property of the transformed fractional Jacobi polynomial. Hence, it is expected that the transformed fractional Jacobi polynomial will be a competitive tool for solving many classes of singular problems.

It is noted that the proposed Galerkin method depends on the transformation of Volterra integral equations, which may encounter difficulties in some practical calculation and increases the computation cost due to the extra computation of integrals. An alternative approach is to apply the Gakerkin method directly to the original problem (Equation (1)). In fact, a direct Galerkin method can be constructed by carefully selecting the local basis functions and weight functions. For example, the approximate solution $y_N\left(t\right)$ can be represented by

$$\begin{array}{c}\hfill y_N\left(t\right)=\sum _{n=0}^{N-1}\left(2\frac{t^{1/\lambda }-t^n}{t_{n+1}-t_n}\sum _{j=0}^m{c}_j^n{\varphi }_{j,n}^{\alpha ,\beta ,\lambda }\left(t\right)+y_N\left(t_n^{\lambda }\right){\mathsf{\Phi }}_{t_n^{\lambda },t_{n+1}^{\lambda }}\left(t\right)\right),t\in (0,T],\end{array}$$

(27)

where ${\mathsf{\Phi }}_{a,b}\left(t\right)$ denotes the characteristic function of the interval $(a,b],$ expressed as

$${\mathsf{\Phi }}_{a,b}\left(t\right)=\left\{\begin{array}{cc}1,\hfill & t\in (a,b],\hfill \\ 0,\hfill & \mathrm{others},\hfill \end{array}\right.$$

In addition, the local basis function ${\varphi }_{j,n}^{\alpha ,\beta ,\lambda }\left(t\right)$ is defined in Section 2. Furthermore, we define the piecewise fractional projection operator ${\mathcal{Q}}^{\lambda ,m,N}$ with respect to the grid $\{t_0^{\lambda },t_1^{\lambda },\cdots ,t_N^{\lambda }\}$ by

$$\left({\mathcal{Q}}^{\lambda ,m,N}f\right)\left(t\right)={\overline{f}}_n\left(t\right),\forall t\in (t_n^{\lambda },t_{n+1}^{\lambda }],$$

where ${\overline{f}}_n\left(t\right)$ satisfies

$${(f,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }={({\overline{f}}_n,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta },k=0,1,\cdots ,m,$$

over the interval $(t_n^{\lambda },t_{n+1}^{\lambda }]$ and ${\overline{f}}_n\left(t\right)\equiv 0$ on $[0,1]\setminus (t_n^{\lambda },t_{n+1}^{\lambda }]$. Additionally, the modified inner product ${(·,·)}_{a,b}^{\lambda ,\alpha ,\beta }$ is defined by

$${(f,g)}_{a,b}^{\lambda ,\alpha ,\beta }={\int }_a^b{(b^{1/\lambda }-t^{1/\lambda })}^{\alpha }{(t^{1/\lambda }-a^{1/\lambda })}^{\beta }f\left(t\right)g\left(t\right)dt.$$

With the help of the operator ${\mathcal{Q}}^{\lambda ,m,N},$ we obtain the direct piecewise fractional Galerkin solution $y_N\left(t\right)$ by

$$\begin{array}{c}\hfill {\mathcal{Q}}^{\lambda ,m,N}{y}_N^{\prime }={\mathcal{Q}}^{\lambda ,m,N}p{y}_N+{\mathcal{Q}}^{\lambda ,m,N}q+{\mathcal{Q}}^{\lambda ,m,N}\mathcal{K}{y}_N,\end{array}$$

(28)

or the equivalent

$$\begin{array}{c}\hfill {(y_N^{\prime },{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }={(p{y}_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }+{(q,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }+{(\mathcal{K}{y}_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }.\end{array}$$

(29)

A direct calculation leads to

$$\begin{array}{cc}\hfill & {(y_N^{\prime },{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{y}_N^{\prime }\left(t\right)P_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)dt\hfill \\ \hfill =& {\left(\frac{h}{2}\right)}^{\alpha +\beta }{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right)\left(\sum _{j=0}^{m}2c_j^n{P}_j^{\alpha ,\beta }\left(x\right)\right)dx\hfill \\ \hfill & +{\left(\frac{h}{2}\right)}^{\alpha +\beta }{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right)\left(\sum _{j=1}^m{c}_j^n(j+\alpha +\beta +1)P_{j-1}^{\alpha +1,\beta +1}\left(x\right)\right)dx,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {(p{y}_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)p\left(t\right)y_N\left(t\right)dt\hfill \\ \hfill =& \lambda \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_n^{\alpha ,\beta }\left(x\right)p\left({\phi }_{n,\lambda }\left(x\right)\right)\left(\sum _{j=0}^m{c}_j^n{P}_j^{\alpha ,\beta }\left(x\right)\right){\left({\phi }_{n,\lambda }\left(x\right)\right)}^{1-1/\lambda }dx\hfill \\ \hfill & +y_N\left(t_n^{\lambda }\right)\lambda \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_n^{\alpha ,\beta }\left(x\right)p\left({\phi }_{n,\lambda }\left(x\right)\right){\left({\phi }_{n,\lambda }\left(x\right)\right)}^{1-1/\lambda }dx,\hfill \end{array}$$

$$\begin{array}{cc}\hfill & {(q,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right)q\left(t\right)dt\hfill \\ \hfill =& \lambda \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{(1-x)}^{\alpha }{(x+1)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right)q\left({\phi }_{n,\lambda }\left(x\right)\right){\left({\phi }_{n,\lambda }\left(x\right)\right)}^{1-1/\lambda }dx,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & {(\mathcal{K}{y}_N,{\psi }_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta })}_{t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\hfill \\ \hfill =& {\int }_{t_n^{\lambda }}^{t_{n+1}^{\lambda }}{\left(t_{n+1}-t^{1/\lambda }\right)}^{\alpha }{\left(t^{1/\lambda }-t_n\right)}^{\beta }{P}_{k,t_n^{\lambda },t_{n+1}^{\lambda }}^{\lambda ,\alpha ,\beta }\left(t\right){\int }_0^t{(t-v)}^{\mu }K(t,v)y_N\left(v\right)dvdt\hfill \\ \hfill =& \lambda \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{\left(1-x\right)}^{\alpha }{\left(x+1\right)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right){\int }_{t_n^{\lambda }}^{{\phi }_{n,\lambda }\left(x\right)}{({\phi }_{n,\lambda }\left(x\right)-v)}^{\mu }K({\phi }_{n,\lambda }\left(x\right),v)y_N\left(v\right)dv{\left({\phi }_{n,\lambda }\left(x\right)\right)}^{1-1/\lambda }dx\hfill \\ \hfill & +\lambda \frac{h^{1+\alpha +\beta }}{2}{\int }_{-1}^1{\left(1-x\right)}^{\alpha }{\left(x+1\right)}^{\beta }{P}_k^{\alpha ,\beta }\left(x\right)\sum _{l=0}^{n-1}{\int }_{t_l^{\lambda }}^{t_{l+1}^{\lambda }}{({\phi }_{n,\lambda }\left(x\right)-v)}^{\mu }K({\phi }_{n,\lambda }\left(x\right),v)y_N\left(v\right)dv{\left({\phi }_{n,\lambda }\left(x\right)\right)}^{1-1/\lambda }dx.\hfill \end{array}$$

With the help of Gauss-type quadrature rules, we can implement the above direct piecewise fractional Galerkin method (DPFG) easily. In

Table 10. Absolute errors and convergence rates of DPFG for Equation (26) when $\mu =-0.2$ and $\lambda =2$.

	$\mathit{m}=2$		$\mathit{m}=4$
	Error	Order	Error	Order
$N=4$	$1.80\times {10}^{-3}$	–	$3.49\times {10}^{-6}$	–
$N=8$	$2.10\times {10}^{-4}$	3.1022	$2.88\times {10}^{-7}$	3.6000
$N=12$	$6.17\times {10}^{-5}$	3.0189	$6.69\times {10}^{-8}$	3.6000
$N=16$	$2.63\times {10}^{-5}$	2.9653	$2.37\times {10}^{-8}$	3.6000
$N=20$	$1.35\times {10}^{-5}$	2.9740	$1.06\times {10}^{-8}$	3.6000
Refer Order		3.0000		3.6000

, we apply DPFG to Equation (26) again. The parameters employed in this numerical example are $\mu =-0.2,\lambda =2$ and $\alpha =\beta =0$, where m and N are both variables. From the computed results, we can conjecture that DPFG shares the same convergence rate with the PFG. However, the convergence analysis remains open.