Singularity separation Chebyshev collocation method for weakly singular Volterra integral equations of the second kind

Abstract

Volterra integral equation of the second kind with weakly singular kernel usually exhibits singular behavior at the origin, which deteriorates the accuracy of standard numerical methods. This paper develops a singularity separation Chebyshev collocation method to solve this kind of Volterra integral equation by splitting the interval into a singular subinterval and a regular one. In the singular subinterval, the general psi-series expansion for the solution about the origin or its Padé approximation is used to approximate the solution. In the regular subinterval, the Chebyshev collocation method is used to discretize the equation. The details of the implementation are also discussed. Specifically, a stable and fast recurrence procedure is derived to evaluate the singular weight integrals involving Chebyshev polynomials analytically. The convergence of the method is proved. We further extend the method to the nonlinear Volterra integral equation by using the Newton method. Three numerical examples are provided to show that the singularity separation Chebyshev collocation method in this paper can effectively solve linear and nonlinear weakly singular Volterra integral equations with high precision.

Introduction

Volterra integral equations with weakly singular kernels arise in many applications in science and technology such as mathematical physics and chemical reactions, including chemical kinetics, heat transfer, heat conduction, shock wave, theory of superfluidity, and boundary layer [1–6].

In this paper, we mainly consider a general linear weakly singular Volterra integral equation (WSVIE) of the second kind

1

$$\begin{array}{cc}& u(x)=f(x)+{\int }_0^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)u(t)\text{d}t,0\le \alpha <1,\mu =0,1,\hfill \end{array}$$where the forcing term f(x) and the kernel $\overline{K}(x,t)$ are continuous functions, which may involve singularities about derivative at $x=0$ or $t=0$ (for $\overline{K}(x,t)$ only), and further $\overline{K}(t,t)\ne 0$. (1) includes two typical kinds of weakly singular algebraic and/or logarithmic kernels ($\mu =0$, $0<\alpha <1$ or $\mu =1$, $\alpha =0$). The theory about the existence, uniqueness, and regularity of the solution has been established [2]. A typical feature of WSVIE is that the derivative of the solution is usually singular at the origin [7, 8], which will deteriorate the convergence orders of standard numerical methods such as Galerkin, collocation, and product integration schemes.

Three popular techniques have been developed to construct efficient numerical methods for solving WSVIEs by considering the singular behavior of the solutions at the origin. The first one uses graded meshes [4, 9–12]. Brunner [9] designed a piecewise polynomial collocation method on graded mesh and developed the convergence theory. Liang and Brunner [11] further proved the convergence of collocation solution on graded mesh in continuous piecewise polynomial space for WSVIE. The second one uses nonpolynomial singular functions in the collocation method, which is suitable for piecewise collocation method [8, 13–16] and spectral collocation method [17–21]. The third one is performing variable transformation for the equation such that the solution is smooth enough. Hence, standard numerical methods, including piecewise polynomial collocation method [22, 23], spectral collocation [24–27] or Galerkin method [28], and product integration method [29–31], can be effectively used to solve WSVIEs. For all the three techniques, it is better to know the asymptotic expansion of the solution near zero. Here, we simply assume $u(x)\sim {\sum }_{k=1}^{\infty }{c}_k{x}^{{\alpha }_k}$, where $0<{\alpha }_1<1$ and ${\alpha }_1<{\alpha }_2<\cdots$. For the method on graded mesh, the nodes are defined by $x_i={(i/n)}^{r}T$, $i=0,1,\dots ,n$, where T is the terminal time. If we choose $r=m/{\alpha }_1$, the expected convergence order m of the piecewise collocation method can theoretically be recovered. But in practical computation, round-off error may deteriorate the accuracy when ${\alpha }_1$ is very small. For the nonpolynomial approximation, a finite-term truncation of the base functions $\{x^{{\alpha }_1},x^{{\alpha }_2},\dots \}$ should be used to construct collocation space, which will enlarge the dimension of the space when convergence order is given. For the smooth transformation, a simple choice is $x=y^r$ in (1) ($\mu =0$) for a suitable exponent $r>1$. A better choice for r is making the first few of $\{r{\alpha }_k\}$ be integers and r should be as small as possible. In addition, the transformation will make the equation more complicated than the original one.

Recently, for the following nonlinear WSVIE of the second kind

2

$$\begin{array}{cc}& u(x)=f(x)+{\int }_0^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}F(x,t,u(t))\text{d}t,0\le \alpha <1,\mu =0,1,\hfill \end{array}$$the first author et al. [32] derived a general asymptotic expansion for the solution about the origin. We summarize the main results in the following lemma.

Lemma 1

([32]) Suppose that F(x, t, u) is infinitely differentiable with respect to x, t, and u, or more generally, $\overline{F}(y,\tau ,u)=F(x,t,u)$ is infinitely differentiable with respect to y, $\tau$, and u, where $y=x^{\gamma }$, $\tau =t^{\omega }$ and $0<\gamma ,\omega \le 1$, and further f(x) possesses the following series expansion about $x=0^{+}$:$$\begin{array}{cc}\hfill f(x)& \sim \sum _{k=1}^{\infty }\sum _{j=1}^{p_k}{f}_{k,j}{x}^{{\beta }_k}{\left(log,x\right)}^{{\varsigma }_{k,j}},x\to 0^{+},\hfill \end{array}$$where ${\beta }_k$ are all reals satisfying $0\le {\beta }_1<{\beta }_2<\cdots \to \infty$, and $p_k$, ${\varsigma }_{k,j}$ are all nonnegative integers, then the solution u(x) of WSVIE (2) holds the finite-term asymptotic expansion about $x=0^{+}$:

3

$$\begin{array}{cc}& u(x)=\sum _{k=1}^N\sum _{j=1}^{q_k}{\eta }_{k,j}{x}^{{\alpha }_k}{(logx)}^{{\phi }_{k,j}}+r_u(x):=u_p(x)+r_u(x),\hfill \end{array}$$where N is a given integer, ${\alpha }_k$ are all reals satisfying $0\le {\alpha }_1<{\alpha }_2<\cdots <{\alpha }_N$, and $q_k$, ${\phi }_{k,j}$ are all nonnegative integers. The expansion $u_p(x)$ can be efficiently recovered by the Picard iteration algorithm, and the remainder $r_u(x)=o\left(x^{{\alpha }_N}\right)$ can be sufficiently smooth by choosing suitably large N.

We note that the expansion (3) can also be derived by Laplace transform [33, 34] for the linear case $\overline{K}(x,t)=1$ in (1). The expansion (3) was called Puiseux expansion in the references [32–34]. It was also referred to as psi-series in [35, 36], which can accurately characterize the singular behavior of the solution involving algebraic and logarithmic singularities at a point. Since the asymptotic expansion of the solution about $x=0$ is a good approximation to the solution when x is small, we can introduce a positive real $\delta$ and split the interval [0, T] as $[0,\delta ]\cup [\delta ,T]$. For $x\in [0,\delta ]$, the asymptotic expansion approximates the solution u(x). Hence, we only need to solve (1) on the remaining regular interval $[\delta ,T]$. Since the solution on this interval is very smooth, we can use standard numerical algorithms such as the piecewise collocation method [7] or spectral collocation method [37, 38] to solve the equation efficiently. In this paper, we discretize the equation using the Chebyshev collocation method and call the method as Singularity Separation Chebyshev Collocation Method (SSCCM).

We further list a few of the works relevant to this paper. Orsi [39] split the interval as a singular interval and a regular one and then solved the WSVIE using the spectral collocation method and product integration method on the singular and regular intervals, respectively. Allaei et al. [40] also split the interval into two subintervals. On the first subinterval, a series solution was obtained by Picard iteration, and on the second subinterval, two product integration methods based on rectangular and trapezoidal rules were analyzed. Huang and Stynes [41] showed that the convergence order is low when the Chebyshev collocation method is directly applied to the WSVIE with a weakly singular solution. Liu and Chen [25] first transformed the WSVIE (1) ($\alpha =1/2$, $\mu =0$) into an equivalent one possessing better regularity and then discretized the equation using the Chebyshev collocation method. Khater et al. [10] solved the WSVIE with only logarithmic kernel by piecewise Chebyshev polynomials on graded meshes. Compared with the above-mentioned works, our proposed method has the following distinct features:

• The kernel of (1) has algebraic and/or logarithmic singularities and $\overline{K}(x,t)$, f(x) may also involve singularities about derivatives at $x=0$ or $t=0$ (for $\overline{K}(x,t)$ only), which means that (1) is a more general WSVIE of the second kind.

• For this more general WSVIE, we use all the singular information to design the Chebyshev collocation method, and the above-mentioned three techniques to treat the singularities are not necessary at all.

• In the implementation of the Chebyshev collocation method on the regular interval, the weight integrals ${\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}{T}_{\rho }(s)\text{d}s$ ($\nu =0,1$, $\rho =0,1,\dots$) are evaluated analytically using a stable and fast recurrence procedure. We note that these singular integrals (without logarithmic singularity) are usually calculated by Jacobi-Gauss quadratures in literature [24, 25], which is unsuitable for evaluating the integrals with logarithmic singularity.

The outline of the paper is as follows. In Section 2, we design a Chebyshev collocation method based on the known psi-series expansion of the solution about the origin and discuss the implementation procedure by using Chebyshev interpolating bases in detail. In Section 3, we analyze the convergence of SSCCM and prove the method has optimal convergence order. In Section 4, we further extend the method to a nonlinear integral equation and show that the Newton iteration combined with SSCCM can efficiently solve nonlinear WSVIE of the second kind. In Section 5, three typical examples are provided to show that SSCCM is accurate and efficient for solving these kinds of linear and nonlinear WSVIEs. In Section 6, we end with a concise conclusion.

The singularity separation Chebyshev collocation method

The typical feature of the solution for WSVIE (1) is that the derivative may be singular at the origin. Lemma 1 shows that we can obtain the truncated psi-series for the solution about the origin by Picard iteration [32]. In this section, we develop a Chebyshev collocation method to solve (1) on a finite interval [0, T] using the psi-series expansion. Since $u_p(x)$ in (3) is exactly the truncated asymptotic expansion for u(x) about $x=0^{+}$, we can conclude that

4

$$\begin{array}{cc}& u(x)-u_p(x)=x^{{\alpha }_N}z(x),\hfill \end{array}$$where z(x) is continuous on the interval [0, T].

Remark 1

In Lemma 1, we have pointed out that the remainder has the form $r_u(x)=u(x)-u_p(x)=o\left(x^{{\alpha }_N}\right)=O\left(x^{{\alpha }_{N+1}}\right)$ (${\alpha }_N<{\alpha }_{N+1}$). Here, we write $r_u(x)=x^{{\alpha }_N}z(x)$ and $z(x)=x^{{\alpha }_{N+1}-{\alpha }_N}(z_0+z_{1}logx+z_2{(logx)}^2+\cdots )$ to guarantee that z(x) is a continuous function on the interval $[0,\delta ]$ ($0<\delta <T$). Then, we can estimate the upper bound of $|z(t)|$ on the interval $[0,\delta ]$, which will be used in Theorem 1.

We know from (4) that $u_p(x)$ can be accurate enough to u(x) near $x=0^{+}$. Hence, we introduce a small positive real $\delta$ and approximate u(t) by its expansion $u_p(t)$ on the subinterval $[0,\delta ]$. We further split (1) as

5

$$\begin{array}{cc}& u(x)=f(x)+G_{\delta }(x)+R_{\delta }(x)+{\int }_{\delta }^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)u(t)\text{d}t,x\in (\delta ,T],\hfill \end{array}$$where

6

$$\begin{array}{cc}& G_{\delta }(x)={\int }_0^{\delta }\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)u_p(t)\text{d}t,\hfill \end{array}$$

7

$$\begin{array}{cc}& R_{\delta }(x)={\int }_0^{\delta }\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)(u(t)-u_p(t))\text{d}t,x\in (\delta ,T].\hfill \end{array}$$We note that the integral in (6) is nearly singular at $t=\delta$ when x is close to $\delta$, which can be effectively evaluated using the following method. If $x-\delta \ge 0.2$, we evaluate it by the modified composite Gauss-Legendre rule [42, 43]. If $0<x-\delta <0.2$, we introduce a $\overline{\delta }=min\{\delta /4,0.2\}$ and split the integral (6) as

8

$$\begin{array}{cc}& G_{\delta }(x)={\int }_0^{\delta -\overline{\delta }}\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)u_p(t)\text{d}t+{\int }_{\delta -\overline{\delta }}^{\delta }\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)u_p(t)\text{d}t.\hfill \end{array}$$The first integral in (8) is evaluated by the Gauss-Legendre rule [42]. For the second integral in (8), a variable change $t=x-s^{\beta }$ yields

9

$$\begin{array}{cc}& {\int }_{\delta -\overline{\delta }}^{\delta }\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)u_p(t)\text{d}t\hfill \\ \hfill & ={\beta }^{1+\mu }{\int }_{{(x-\delta )}^{1/\beta }}^{{(x-\delta +\overline{\delta })}^{1/\beta }}{s}^{\beta (1-\alpha )-1}{log}^{\mu }s\overline{K}(x,x-s^{\beta })u_p(x-s^{\beta })\text{d}s.\hfill \end{array}$$If we choose $\beta =\frac{2}{1-\alpha }$, the integral on the right-hand side of (9) is no longer nearly singular, which can be evaluated by the standard Gauss-Legendre rule.

Theorem 1

Suppose that $\overline{K}(x,t)$ is continuously differentiable with respect to $x\in [\delta ,T]$ and continuous with respect to $t\in [0,T]$, and $u_p(x)$ defined in (3) is  the truncated psi-series for u(x) about $x=0$, then there exists a positive number $\sigma$ such that

10

$$\begin{array}{cc}& |R_{\delta }(x)|\le {\overline{K}}_{max}{z}_{max}\frac{1}{{(\sigma \text{e})}^{\mu }}\text{B}(1+{\alpha }_N,1-\alpha -\sigma \mu ){\delta }^{{\alpha }_N+1-\alpha -\sigma \mu },x\in (\delta ,T],\hfill \end{array}$$

11

$$\begin{array}{cc}& |R_{\delta }(x)-R_{\delta }(y)|\le \left\{\begin{array}{c}{C}_{\delta ,0}{|x-y|}^{1-\alpha }+C_{\delta ,0}^{\ast }|x-y|,x,y\in (\delta ,T],\mu =0,\hfill \\ C_{\delta ,1}{|x-y|}^{1-\alpha -\sigma }+C_{\delta ,1}^{\ast }|x-y|,x,y\in (\delta ,{\text{e}}^{1/\alpha }],\mu =1,\hfill \end{array}\right)\hfill \end{array}$$where ${\overline{K}}_{max}={max}_{x\in [\delta ,T],t\in [0,\delta ]}|\overline{K}(x,t)|$, $z_{max}={max}_{t\in [0,\delta ]}|z(t)|$, $\sigma$ satisfies $0<\alpha +\sigma \mu <1$, $\text{B}(\xi ,\eta )$ is the beta function, and$$\begin{array}{cc}& C_{\delta ,0}={\overline{K}}_{max}{z}_{max}{\delta }^{{\alpha }_N},C_{\delta ,1}={\overline{K}}_{max}{z}_{max}\frac{2{\delta }^{{\alpha }_N}}{\sigma \text{e}},\hfill \\ \hfill & C_{\delta ,\mu }^{\ast }={\overline{K}}_{max}^x\frac{z_{max}}{{(\sigma \text{e})}^{\mu }}\text{B}(1+{\alpha }_N,1-\alpha -\sigma \mu ){\delta }^{{\alpha }_N+1-\alpha -\sigma \mu },\hfill \\ \hfill & {\overline{K}}_{max}^x=\underset{x\in [\delta ,T],t\in [0,\delta ]}{max}\left|\frac{\partial \overline{K}(x,t)}{\partial x}\right|.\hfill \end{array}$$

Proof

Using (4) and (7), we have when $x\in (\delta ,T]$$$\begin{array}{cc}\hfill |R_{\delta }(x)|& \le {\overline{K}}_{max}{z}_{max}{\int }_0^{\delta }\frac{|{log}^{\mu }(x-t)|}{{(x-t)}^{\alpha }}{t}^{{\alpha }_N}\text{d}t.\hfill \end{array}$$Applying the inequality $|logx|\le x^{-\sigma }/(\sigma \text{e})$ ($x>0$, $\sigma >0$) to ${log}^{\mu }(x-t)$, we know there exists a positive number $\sigma$ satisfying $0<\alpha +\sigma \mu <1$, such that

12

$$\begin{array}{cc}\hfill |R_{\delta }(x)|& \le {\overline{K}}_{max}{z}_{max}\frac{1}{{(\sigma \text{e})}^{\mu }}{\int }_0^{\delta }\frac{t^{{\alpha }_N}}{{(x-t)}^{\alpha +\sigma \mu }}\text{d}t(\text{using}t=\delta s)\hfill \\ \hfill & ={\overline{K}}_{max}{z}_{max}\frac{{\delta }^{{\alpha }_N+1-\alpha -\sigma \mu }}{{(\sigma \text{e})}^{\mu }}{\int }_0^1{s}^{{\alpha }_N}{\left(\frac{x}{\delta },-,s\right)}^{-\alpha -\sigma \mu }\text{d}s\hfill \\ \hfill & \le {\overline{K}}_{max}{z}_{max}\frac{{\delta }^{{\alpha }_N+1-\alpha -\sigma \mu }}{{(\sigma \text{e})}^{\mu }}{\int }_0^1{s}^{{\alpha }_N}{\left(1,-,s\right)}^{-\alpha -\sigma \mu }\text{d}s,\hfill \end{array}$$which implies that the estimate (10) holds. We note that the last inequality of (12) holds since $x>\delta$. As for (11), we have$$\begin{array}{cc}\hfill R_{\delta }(x)-R_{\delta }(y)=& {\int }_0^{\delta }\left(\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }},-,\frac{{log}^{\mu }(y-t)}{{(y-t)}^{\alpha }}\right)\overline{K}(x,t)t^{{\alpha }_N}z(t)\text{d}t\hfill \\ \hfill & +{\int }_0^{\delta }\frac{{log}^{\mu }(y-t)}{{(y-t)}^{\alpha }}\left(\overline{K},(,x,,,t,),-,\overline{K},(,y,,,t,)\right)t^{{\alpha }_N}z(t)\text{d}t.\hfill \end{array}$$For $\mu =1$, since $\omega (y)=y^{-\alpha }logy$ is a strictly monotone increasing function on $(0,{\text{e}}^{1/\alpha }]$, we know when $0<\delta <x<y\le {\text{e}}^{1/\alpha }$

13

$$\begin{array}{cc}\hfill |R_{\delta }(x)-R_{\delta }(y)|& \le {\overline{K}}_{max}{z}_{max}{\delta }^{{\alpha }_N}\left|{\int }_0^{\delta },\left(\frac{log(y-t)}{{(y-t)}^{\alpha }},-,\frac{log(x-t)}{{(x-t)}^{\alpha }}\right),\text{d},t\right|\hfill \\ \hfill & +{\overline{K}}_{max}^x{z}_{max}|x-y|{\int }_0^{\delta }\left|\frac{log(y-t)}{{(y-t)}^{\alpha }}\right|t^{{\alpha }_N}\text{d}t.\hfill \end{array}$$Direct integrations show$$\begin{array}{cc}& {\int }_0^{\delta }\left(\frac{log(y-t)}{{(y-t)}^{\alpha }},-,\frac{log(x-t)}{{(x-t)}^{\alpha }}\right)\text{d}t={\int }_x^y\frac{logt}{t^{\alpha }}\text{d}t-{\int }_x^y\frac{log(t-\delta )}{{(t-\delta )}^{\alpha }}\text{d}t.\hfill \end{array}$$Thus,

14

$$\begin{array}{cc}\hfill \left|{\int }_0^{\delta },\left(\frac{log(y-t)}{{(y-t)}^{\alpha }},-,\frac{log(x-t)}{{(x-t)}^{\alpha }}\right),\text{d},t\right|& \le \frac{2}{\sigma \text{e}}{\int }_x^y\frac{1}{{(t-x)}^{\alpha +\sigma }}\text{d}t=\frac{2}{\sigma \text{e}}{(y-x)}^{1-\alpha -\sigma }.\hfill \end{array}$$Substituting (14) into (13), and noticing (10), we can conclude that (11) holds for the case $\mu =1$. The case $\mu =0$ can be proved analogously. The theorem is proved.$\square$

Theorem 1 tells us that the truncation error $R_{\delta }(x)$ can be as small as possible by choosing ${\alpha }_N$ suitably large and $\delta$ suitably small.

Remark 2

In practical computation, it is preferred that $u_p(t)$ in (3) is replaced by its Padé approximation $u_r(t)$ [44] to extend the convergence range, which means that we can choose a larger $\delta$. As for the selection of $\delta$, we define an error indicator

15

$$\begin{array}{cc}& e_r(x)=u_r(x)-f(x)-{\int }_0^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\overline{K}(x,t)u_r(t)\text{d}t.\hfill \end{array}$$Then, we choose $\delta >0$ as large as possible such that $|e_r(\delta )|\le \epsilon$, and $\epsilon$ is the prescribed precision.

Next, we discretize (5) on the finite interval $[\delta ,T]$ by the Chebyshev collocation method. By taking two variable transformations$$\begin{array}{cc}& t=\frac{T-\delta }{2}s+\frac{T+\delta }{2},x=\frac{T-\delta }{2}y+\frac{T+\delta }{2},\hfill \end{array}$$we have $t,x\in [\delta ,T]\to s,y\in [-1,1]$. Further by letting $v(y)=u(x)$, $g(y)=f(x)+G_{\delta }(x)$, $K(y,s)=\overline{K}(x,t)$, and $r_{\delta }(y)=R_{\delta }(x)$, we can transform (5) to

16

$$\begin{array}{cc}\hfill v(y)=& g(y)+r_{\delta }(y)+{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }\hfill \\ \hfill & \times {\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}K(y,s)v(s)\text{d}s,y\in [-1,1].\hfill \end{array}$$We choose the zeros of the Chebyshev polynomial $T_n(s)=cos(narccoss)$ as collocation points. They are$$\begin{array}{c}\hfill s_j=-cos{\theta }_j,\text{and}{\theta }_j=\frac{(2j-1)\pi }{2n},j=1,2,\dots ,n.\end{array}$$For the solution v(s), $s\in [-1,1]$, we construct a Lagrange interpolating polynomial ${\pi }_{n}v(s)$ of degree $n-1$ using the nodes $s_j$, $j=1,2,\dots ,n$, which is

17

$$\begin{array}{cc}& {\pi }_{n}v(s)=\sum _{j=1}^n{b}_j(s)v(s_j),\hfill \end{array}$$where the Lagrange basis functions $b_j(s)$ are$$\begin{array}{cc}& b_j(s)=\prod _{\rho =1,\rho \ne j}^n\frac{s-s_{\rho }}{s_j-s_{\rho }}=\frac{T_n(s)}{(s-s_j)T_n^{\prime }(s_j)},j=1,2,\dots ,n.\hfill \end{array}$$We can also reformulate the Lagrange interpolation (17) as a sum of Chebyshev polynomials in the form [45]

18

$$\begin{array}{cc}& {\pi }_{n}v(s)=\frac{c_0}{2}+\sum _{j=1}^{n-1}{c}_j{T}_j(s),\text{and}{c}_j=\frac{2}{n}\sum _{k=1}^{n}v(s_k)T_j(s_k).\hfill \end{array}$$Replacing v(s) in (16) by its interpolation (17) and dropping away the truncation error $r_{\delta }(y)$ yield

19

$$\begin{array}{cc}\hfill v(y)\approx & g(y)+{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }\sum _{j=1}^n{\overline{w}}_j(\nu ,y)v(s_j),\hfill \end{array}$$where

20

$$\begin{array}{cc}& {\overline{w}}_j(\nu ,y)={\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}K(y,s)b_j(s)\text{d}s,j=1,2,\dots ,n.\hfill \end{array}$$Collocating y in (19) with $s_i$, and denoting by $v_i$ the approximate value of $v(s_i)$, we obtain a Chebyshev collocation scheme

21

$$\begin{array}{cc}\hfill v_i=& g(s_i)+{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{j=1}^n\left[\sum _{\nu =0}^{\mu },\left(\begin{array}{c}\mu \\ \nu \end{array}\right),{\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu },{\overline{w}}_j,(\nu ,s_i)\right]v_j,\hfill \\ & i=1,2,\dots ,n.\hfill \end{array}$$Once $v_i$, $i=1,2,\dots ,n$ are determined, we obtain an approximate solution for (16)

22

$$\begin{array}{cc}& v_n(y)=\sum _{j=1}^n{b}_j(y)v_j.\hfill \end{array}$$The linear system (21) is derived by separating the singularity on a small subinterval and discretizing the equation on the remaining regular subinterval using the Chebyshev collocation method. As mentioned in the “Introduction”, this method is called the singularity separation Chebyshev collocation method (SSCCM) for WSVIE of the second kind.

The evaluation of the weight integrals ${\overline{w}}_j(\nu ,y)$ in (20) is not an easy task. Usually, we can approximate K(y, s) by its Chebyshev interpolation ${\pi }_{n,s}K(y,s)$ to yield

23

$$\begin{array}{cc}& {\overline{w}}_j(\nu ,y)\approx \sum _{k=1}^{n}K(y,s_k){\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}{b}_k(s)b_j(s)\text{d}s.\hfill \end{array}$$In our practical test, we find that the evaluation of the integrals in (23) is very time-consuming. In the following, we develop an efficient technique to implement SSCCM. Our start point is using the interpolation formula (18). By denoting

24

$$\begin{array}{cc}& w_j(\nu ,y)={\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}K(y,s)T_j(s)\text{d}s,j=0,1,\dots ,n-1,\hfill \end{array}$$and replacing v(s) in (16) by its interpolation (18), we obtain the following scheme

25

$$\begin{array}{cc}\hfill \frac{c_0}{2}& +\sum _{j=1}^{n-1}{c}_j{T}_j(s_i)=g(s_i)+{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }\hfill \\ & \times \left(w_0,(\nu ,s_i),\frac{c_0}{2},+,\sum _{j=1}^{n-1},w_j,(\nu ,s_i),c_j\right),i=1,2,\dots ,n,\hfill \end{array}$$which is a linear system of equations with unknowns $c_0,c_1,\dots ,c_{n-1}$. Once $c_j$, $j=0,1,\dots ,n-1$ are obtained, we can use the interpolation in (18) to compute $v_i$, the approximate values of v(s) at the nodes $s_i$, $i=1,2,\dots ,n$. Then, the approximate solution $v_n(y)$ is obtained via (22).

Next, we discuss the efficient evaluation of the weights $w_j(\nu ,y)$ in (24). These integrals are weakly singular and oscillatory, which are better to be evaluated by analytical methods. For this purpose, we apply the Chebyshev interpolation formula (18) to K(y, s) with respect to the argument s to yield

26

$$\begin{array}{cc}& K(y,s)\approx {\pi }_{n,s}K(y,s)=\frac{d_0(y)}{2}+\sum _{m=1}^{n-1}{d}_m(y)T_m(s),\hfill \end{array}$$where for $m=0,1,\dots ,n-1$,$$\begin{array}{cc}& d_m(y)=\frac{2}{n}\sum _{k=1}^{n}K(y,s_k)T_m(s_k),\text{and}{T}_m(s_k)={(-1)}^{m}cos(m{\theta }_k).\hfill \end{array}$$Substituting (26) into (24) yields

27

$$\begin{array}{cc}& w_j(\nu ,y)\approx {\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}\left(\frac{d_0(y)}{2},+,\sum _{m=1}^{n-1},d_m,(y),T_m,(s)\right)T_j(s)\text{d}s.\hfill \end{array}$$Using the formula [45]$$\begin{array}{cc}& T_m(s)T_j(s)=\frac{1}{2}\left(T_{m+j},(s),+,T_{|m-j|},(s)\right),\hfill \end{array}$$and introducing the notation

28

$$\begin{array}{cc}\hfill I_{\rho }& (\nu ,y)={\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}{T}_{\rho }(s)\text{d}s,\hfill \end{array}$$we can convert (27) to

29

$$\begin{array}{cc}& w_j(\nu ,y)\approx \frac{d_0(y)}{2}{I}_j(\nu ,y)+\frac{1}{2}\sum _{m=1}^{n-1}{d}_m(y)\left(I_{m+j},(\nu ,y),+,I_{|m-j|},(\nu ,y)\right).\hfill \end{array}$$

Theorem 2

For the integrals defined by (28), we have the following recurrences$$\begin{array}{cc}& I_0(\nu ,y)={(1+y)}^{1-\alpha }\sum _{l=0}^{\nu }\frac{{(-1)}^l{(\nu )}_l}{{(1-\alpha )}^{l+1}}{\left(log,(,1,+,y,)\right)}^{\nu -l},\hfill \\ \hfill & I_1(\nu ,y)=\frac{y}{2-\alpha }{I}_0(\nu ,y)-\frac{{(1+y)}^{1-\alpha }}{2-\alpha }{log}^{\nu }(1+y)-\frac{\nu }{2-\alpha }{I}_1(\nu -1,y),\hfill \\ \hfill & I_2(\nu ,y)=\frac{4y}{3-\alpha }{I}_1(\nu ,y)-I_0(\nu ,y)+\frac{2{(1+y)}^{1-\alpha }}{3-\alpha }{log}^{\nu }(1+y)\hfill \\ \hfill & -\frac{\nu }{3-\alpha }\left[I_2,(\nu -1,y),+,I_0,(\nu -1,y)\right],\hfill \end{array}$$and for $\rho \ge 2$$$\begin{array}{cc}\hfill I_{\rho +1}(\nu ,y)=& \frac{2(\rho +1)y}{2-\alpha +\rho }{I}_{\rho }(\nu ,y)+\frac{(\rho +1)(2-\alpha -\rho )}{(\rho -1)(2-\alpha +\rho )}{I}_{\rho -1}(\nu ,y)\hfill \\ \hfill & +{(-1)}^{\rho }\frac{2{(1+y)}^{1-\alpha }{log}^{\nu }(1+y)}{(\rho -1)(2-\alpha +\rho )}\hfill \\ \hfill & -\frac{\nu }{2-\alpha +\rho }\left(I_{\rho +1},(\nu -1,y),-,\frac{\rho +1}{\rho -1},I_{\rho -1},(\nu -1,y)\right),\hfill \end{array}$$where ${(\nu )}_l$ are defined by ${(\nu )}_0=1$, ${(\nu )}_l={(\nu )}_{l-1}(\nu -l+1)$, and $I_{\rho }(-1,y)=0$.

Proof

The simple case of $\nu =0$ in this theorem has been proved in [46]. Here, we prove the general case with a logarithmic factor. $I_0(\nu ,y)$ is obtained by integrating the integrand directly. For $I_1(\nu ,y)$, we have$$\begin{array}{cc}\hfill I_1(\nu ,y)=& -{\int }_{-1}^y{(y-s)}^{1-\alpha }{log}^{\nu }(y-s)\text{d}s+y{I}_0(\nu ,y)\hfill \\ \hfill =& -{(y-s)}^{1-\alpha }{log}^{\nu }(y-s)s{|}_{-1}^y-(1-\alpha )I_1(\nu ,y)\hfill \\ \hfill & -\nu I_1(\nu -1,y)+y{I}_0(\nu ,y),\hfill \end{array}$$from which we know the recurrence for $I_1(\nu ,y)$ holds. The recurrence for $I_2(\nu ,y)$ can be proved similarly to the general case. Now, we consider the general case $\rho \ge 2$. We need two identities for Chebyshev polynomials [45, 47]$$\begin{array}{cc}& T_{\rho +1}(s)=2s{T}_{\rho }(s)-T_{\rho -1}(s),\hfill \\ \hfill & \int T_{\rho }(s)\text{d}s=\frac{1}{2}\left(\frac{T_{\rho +1}(s)}{\rho +1},-,\frac{T_{\rho -1}(s)}{\rho -1}\right),\rho \ge 2.\hfill \end{array}$$Then,$$\begin{array}{cc}& I_{\rho +1}(\nu ,y)=2{\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}s{T}_{\rho }(s)\text{d}s-I_{\rho -1}(\nu ,y)\hfill \\ \hfill & =-2{\int }_{-1}^y{(y-s)}^{1-\alpha }{log}^{\nu }(y-s)T_{\rho }(s)\text{d}s+2y{I}_{\rho }(\nu ,y)-I_{\rho -1}(\nu ,y)\hfill \\ \hfill & ={(-1)}^{\rho }\frac{2{(1+y)}^{1-\alpha }}{{\rho }^2-1}{log}^{\nu }(1+y)-(1-\alpha ){\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}\left(\frac{T_{\rho +1}(s)}{\rho +1},,-,,\frac{T_{\rho -1}(s)}{\rho -1}\right)\text{d}s\hfill \\ \hfill & -\nu {\int }_{-1}^y\frac{{log}^{\nu -1}(y-s)}{{(y-s)}^{\alpha }}\left(\frac{T_{\rho +1}(s)}{\rho +1},-,\frac{T_{\rho -1}(s)}{\rho -1}\right)\text{d}s+2y{I}_{\rho }(\nu ,y)-I_{\rho -1}(\nu ,y)\hfill \\ \hfill & ={(-1)}^{\rho }\frac{2{(1+y)}^{1-\alpha }}{{\rho }^2-1}{log}^{\nu }(1+y)-\frac{1-\alpha }{\rho +1}{I}_{\rho +1}(\nu ,y)+\frac{2-\alpha -\rho }{\rho -1}{I}_{\rho -1}(\nu ,y)\hfill \\ \hfill & -\nu \left(\frac{1}{\rho +1},I_{\rho +1},(\nu -1,y),-,\frac{1}{\rho -1},I_{\rho -1},(\nu -1,y)\right)+2y{I}_{\rho }(\nu ,y),\hfill \end{array}$$from which we can prove the three-term recurrence for $I_{\rho }(\nu ,y)$, $\rho \ge 2$.

Remark 3

By test, we find the recurrence formulas in Theorem 2 are fast and stable to evaluate the integrals $I_{\rho }(\nu ,y)$. For instance, we obtain $I_{100}(0,0.99)=-0.41137631093062577$ with precision $3.4972\times {10}^{-15}$ when $\alpha =\frac{4}{5}$.

For the integrals $I_{\rho }(\nu ,y)$, $\nu =0,1$, we evaluate first $I_{\rho }(0,s_i)$, $\rho =0,1,\dots ,2(n-1)$, then $I_{\rho }(1,s_i)$, $\rho =0,1,\dots ,2(n-1)$ using the recurrences in Theorem 2 if necessary, where $i=1,2,\dots ,n$. Once all $I_{\rho }(\nu ,s_i)$ are evaluated, we can obtain the approximations of the weights $w_j(\nu ,s_i)$ via (29). We note that the evaluation of the weights $w_j(\nu ,s_i)$ can be simplified when K(y, s) is a constant, or more generally, it is a function for y only.

Convergence analysis

In this section, we discuss the convergence of SSCCM (21) combining with (23). Using (22), we rewrite the scheme (21) as

30

$$\begin{array}{cc}\hfill v_i=g(s_i)& +{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }\hfill \\ & \times {\int }_{-1}^{s_i}\frac{{log}^{\nu }(s_i-s)}{{(s_i-s)}^{\alpha }}{\pi }_{n,s}K(s_i,s)v_n(s)\text{d}s.\hfill \end{array}$$Letting $y=s_i$ in (16) and subtracting (30), we get

31

$$\begin{array}{cc}\hfill v(s_i)-v_i=& r_{\delta }(s_i)+{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }\hfill \\ & \times {\int }_{-1}^{s_i}\frac{{log}^{\nu }(s_i-s)}{{(s_i-s)}^{\alpha }}\left[K,(s_i,s),v,(s),-,{\pi }_{n,s},K,(s_i,s),v_n,(s)\right]\text{d}s.\hfill \end{array}$$Denote the error function by $e(y)=v(y)-v_n(y)$. Multiplying both sides of (31) by $b_i(y)$ and summing over i from 1 to n yield the error equation

32

$$\begin{array}{cc}& {\pi }_{n}v(y)-v_n(y)={\pi }_n{r}_{\delta }(y)+{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }\sum _{i=1}^n{b}_i(y)\hfill \\ & \times {\int }_{-1}^{s_i}\frac{{log}^{\nu }(s_i-s)}{{(s_i-s)}^{\alpha }}\left[(K(s_i,s)-{\pi }_{n,s}K(s_i,s)),v_n,(s),+,K,(s_i,s),e,(s)\right]\text{d}s.\hfill \end{array}$$By denoting

33

$$\begin{array}{cc}& {\widehat{\theta }}_{\nu }(y)={\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}(K(y,s)-{\pi }_{n,s}K(y,s))v_n(s)\text{d}s,\hfill \end{array}$$

34

$$\begin{array}{cc}& {\overline{\theta }}_{\nu }(y)={\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}K(y,s)e(s)\text{d}s,\hfill \end{array}$$we can further write the error (32) as

35

$$\begin{array}{cc}\hfill e(y)=& J_1(y)+J_2(y)+J_3(y)+J_4(y)+J_5(y),\hfill \end{array}$$where$$\begin{array}{cc}& J_1(y)=v(y)-{\pi }_{n}v(y),\hfill \\ & J_2(y)={\pi }_n{r}_{\delta }(y),\hfill \\ & J_3(y)={\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }{\pi }_n{\widehat{\theta }}_{\nu }(y),\hfill \\ & J_4(y)={\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }({\pi }_n{\overline{\theta }}_{\nu }(y)-{\overline{\theta }}_{\nu }(y)),\hfill \\ & J_5(y)={\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,\frac{T-\delta }{2}\right)}^{\mu -\nu }{\overline{\theta }}_{\nu }(y).\hfill \end{array}$$Using the technique in estimating (12), we can show that there exists a positive number $\sigma$ satisfying $0<\alpha +\sigma \nu <1$, such that$$\begin{array}{c}\hfill |{\overline{\theta }}_{\nu }(y)|\le \frac{1}{{(\sigma \text{e})}^{\nu }}\underset{y,s\in [-1,1]}{max}|K(y,s)|{\int }_{-1}^y\frac{|e(s)|}{{(y-s)}^{\alpha +\sigma \nu }}\text{d}s.\end{array}$$Then, the generalized Gronwall lemma ([7], P359, Theorem 6.1.17) implies

36

$$\begin{array}{cc}\hfill {\Vert e\Vert }_{\infty }\le & C(\Vert J_1{\Vert }_{\infty }+\Vert J_2{\Vert }_{\infty }+\Vert J_3{\Vert }_{\infty }+\Vert J_4{\Vert }_{\infty }),\hfill \end{array}$$where the maximum norm ${\Vert e\Vert }_{\infty }={max}_{y\in [-1,1]}|e(y)|$, and C is a positive constant independent of e(y).

Noting that (33) and (34) have the same structure, we discuss their properties similarly by denoting

37

$$\begin{array}{cc}& {\theta }_{\nu }(y)={\int }_{-1}^y\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }}Q(y,s)\text{d}s,\nu =0,1.\hfill \end{array}$$We point out that ${\theta }_{\nu }(y)\in C^{0,1-\alpha -\sigma \nu }[-1,1]$ ($0<1-\alpha -\sigma \nu <1$), which can be seen by taking $s=((y+1)\tau +(y-1))/2$ in (37). Here, $C^{0,\kappa }[-1,1]$ ($0<\kappa <1$) is the Hölder continuous function space with exponent $\kappa$. Its norm is usually defined by

38

$$\begin{array}{c}\hfill \Vert {\theta }_{\nu }{\Vert }_{0,\kappa }=\underset{y\in [-1,1]}{max}|{\theta }_{\nu }(y)|+\underset{y\ne z}{sup}\frac{|{\theta }_{\nu }(y)-{\theta }_{\nu }(z)|}{{|y-z|}^{\kappa }}.\end{array}$$

Lemma 2

For the integral defined by (37), there exist general positive constants $C_1,C_2$ independent of Q(y, s), such that

39

$$\begin{array}{cc}& |{\theta }_{\nu }(y)|\le \frac{1}{{(\sigma \text{e})}^{\nu }}\frac{2^{1-\alpha -\sigma \nu }}{1-\alpha -\sigma \nu }\underset{y,s\in [-1,1]}{max}|Q(y,s)|,\hfill \end{array}$$

40

$$\begin{array}{cc}& \frac{|{\theta }_{\nu }(y)-{\theta }_{\nu }(z)|}{{|y-z|}^{\kappa }}\le C_1\underset{y,s\in [-1,1]}{max}|Q(y,s)|+C_2\underset{y,s\in [-1,1]}{max}|Q_y^{\prime }(y,s)|,\hfill \end{array}$$where $\nu =0,1$, $\sigma >0$ satisfies $1-\alpha -\sigma \nu >0$, and $0<\kappa \le 1-\alpha -\sigma \nu$. Furthermore, we have

41

$$\begin{array}{c}\hfill \Vert {\theta }_{\nu }{\Vert }_{0,\kappa }\le C_1\underset{y,s\in [-1,1]}{max}|Q(y,s)|+C_2\underset{y,s\in [-1,1]}{max}|Q_y^{\prime }(y,s)|.\end{array}$$

Proof

The inequality (39) holds obviously. We only prove (40). For $\forall y,z\in [-1,1]$, without loss of generality, we assume $y<z$. Then,

42

$$\begin{array}{cc}& |{\theta }_{\nu }(y)-{\theta }_{\nu }(z)|\le \left|{\int }_{-1}^y,\left[\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }},-,\frac{{log}^{\nu }(z-s)}{{(z-s)}^{\alpha }}\right],Q,(y,s),\text{d},s\right|\hfill \\ \hfill & +\left|{\int }_{-1}^y,\frac{{log}^{\nu }(z-s)}{{(z-s)}^{\alpha }},[Q(y,s)-Q(z,s)],\text{d},s\right|+\left|{\int }_y^z,\frac{{log}^{\nu }(z-s)}{{(z-s)}^{\alpha }},Q,(z,s),\text{d},s\right|\hfill \\ \hfill & :=E_1+E_2+E_3.\hfill \end{array}$$For the term $E_1$, we have

43

$$\begin{array}{cc}& E_1\le \underset{y,s\in [-1,1]}{max}|Q(y,s)|{\int }_{-1}^y\left|\frac{{log}^{\nu }(y-s)}{{(y-s)}^{\alpha }},-,\frac{{log}^{\nu }(z-s)}{{(z-s)}^{\alpha }}\right|\text{d}s.\hfill \end{array}$$For the case $\nu =0$, since $\omega (y)=y^{-\alpha }$ is a strictly monotone decreasing function for $y>0$, we know$$\begin{array}{cc}\hfill {\int }_{-1}^y\left|\frac{1}{{(y-s)}^{\alpha }},-,\frac{1}{{(z-s)}^{\alpha }}\right|\text{d}s& =\frac{1}{1-\alpha }\left[{(z-y)}^{1-\alpha },-,{(z+1)}^{1-\alpha },+,{(y+1)}^{1-\alpha }\right]\hfill \\ \hfill & \le \frac{1}{1-\alpha }{(z-y)}^{1-\alpha }.\hfill \end{array}$$Hence,

44

$$\begin{array}{cc}& E_1\le \frac{1}{1-\alpha }\underset{y,s\in [-1,1]}{max}|Q(y,s)|{(z-y)}^{1-\alpha }\text{for}\nu =0.\hfill \end{array}$$For the case $\nu =1$, it is easy to show that $\omega (y)=y^{-\alpha }logy$ is a strictly monotone increasing function on (0, 2]. A direct integration shows

45

$$\begin{array}{cc}& {\int }_{-1}^y\left[\frac{log(z-s)}{{(z-s)}^{\alpha }},-,\frac{log(y-s)}{{(y-s)}^{\alpha }}\right]\text{d}s\hfill \\ \hfill & =\frac{{(z-y)}^{1-\alpha }}{{(1-\alpha )}^2}-\frac{{(z-y)}^{1-\alpha }}{1-\alpha }log(z-y)+{\int }_y^z\frac{log(s+1)}{{(s+1)}^{\alpha }}\text{d}s.\hfill \end{array}$$Applying $|log(s+1)|\le \frac{{(s+1)}^{-\sigma }}{\sigma \text{e}}$ and $s+1\ge s-y(y\ge -1)$, we obtain

46

$$\begin{array}{cc}\hfill \left|{\int }_y^z,\frac{log(s+1)}{{(s+1)}^{\alpha }},\text{d},s\right|& \le \frac{1}{\sigma \text{e}}{\int }_y^z\frac{\text{d}s}{{(s+1)}^{\alpha +\sigma }}\hfill \\ \hfill & \le \frac{1}{\sigma \text{e}}{\int }_y^z\frac{\text{d}s}{{(s-y)}^{\alpha +\sigma }}=\frac{{(z-y)}^{1-\alpha -\sigma }}{\sigma \text{e}(1-\alpha -\sigma )}.\hfill \end{array}$$Combining (43), (45), and (46), we know

47

$$\begin{array}{cc}& E_1\le C_1^{\ast }\underset{y,s\in [-1,1]}{max}|Q(y,s)|{(z-y)}^{1-\alpha -\sigma }\text{for}\nu =1,\hfill \end{array}$$where$$\begin{array}{cc}& C_1^{\ast }=\frac{2^{\sigma }}{{(1-\alpha )}^2}+\frac{1}{1-\alpha }max\left\{\frac{1}{\sigma \text{e}},,,2^{\sigma },log,2\right\}+\frac{1}{\sigma \text{e}(1-\alpha -\sigma )}.\hfill \end{array}$$Obviously, (44) and (47) can be uniformly written as

48

$$\begin{array}{cc}& E_1\le C_1\underset{y,s\in [-1,1]}{max}|Q(y,s)|{(z-y)}^{1-\alpha -\sigma \nu }.\hfill \end{array}$$Furthermore, it is easy to show that $E_2$ and $E_3$ in (42) can be estimated by

49

$$\begin{array}{cc}\hfill E_2& =\left|{\int }_{-1}^y,\frac{{log}^{\nu }(z-s)}{{(z-s)}^{\alpha }},[Q(y,s)-Q(z,s)],\text{d},s\right|\hfill \\ \hfill & \le \underset{y,s\in [-1,1]}{max}|Q_y^{\prime }(y,s)||z-y|\frac{{(z+1)}^{1-\alpha -\sigma \nu }-{(z-y)}^{1-\alpha -\sigma \nu }}{1-\alpha -\sigma \nu },\hfill \end{array}$$

50

$$\begin{array}{cc}\hfill E_3& =\left|{\int }_y^z,\frac{{log}^{\nu }(z-s)}{{(z-s)}^{\alpha }},Q,(z,s),\text{d},s\right|\le \underset{y,s\in [-1,1]}{max}|Q(y,s)|\frac{{(z-y)}^{1-\alpha -\sigma \nu }}{1-\alpha -\sigma \nu }.\hfill \end{array}$$Combining (42) with (48)–(50), we know there exist positive constants $C_1,C_2$, such that the inequality (40) holds for $\forall \kappa \in (0,1-\alpha -\sigma \nu ]$. Noticing the norm definition (38), we know that (41) also holds, but $C_1$ has been changed. The proof is complete.

We also need the following lemmas to complete the convergence analysis.

Lemma 3

([38], P192, Lemma 5.1) For $\forall \theta (y)\in C^{0,\kappa }[-1,1]$ with $0<\kappa <1$, there exist a polynomial $p_n(y)$ of degree n and a positive constant $C_{0,\kappa }$, such that

51

$$\begin{array}{cc}& \Vert \theta -p_n{\Vert }_{\infty }\le C_{0,\kappa }{n}^{-\kappa }{\Vert \theta \Vert }_{0,\kappa }.\hfill \end{array}$$

Lemma 4

([45], Theorem 6.13) For the Lagrange basis functions $b_j(s)$ in the Chebyshev interpolation (17), we have

52

$$\begin{array}{cc}& \underset{s\in [-1,1]}{sup}\sum _{j=1}^n|b_j(s)|=\frac{2}{\pi }log(n-1)+0.9625+\mathcal{O}\left(\frac{1}{n}\right),n\to \infty .\hfill \end{array}$$

Theorem 3

Suppose that the kernel function K(y, s) is sufficiently smooth with respect to $y,s\in [-1,1]$, $v(y)\in C^n[-1,1]$ is the exact solution of WSVIE (16), and $v_n(y)$ defined by (22) is the approximate solution of SSCCM, then when n is suitably large, the error $e(y)=v(y)-v_n(y)$ is estimated by

53

$$\begin{array}{cc}\hfill {\Vert e\Vert }_{\infty }& \le \frac{C}{n!2^{n-1}}||v^{(n)}{||}_{\infty }+C{(n-1)}^{-\kappa }log(n-1){\delta }^{{\alpha }_N}+C{\delta }^{{\alpha }_N+1-\alpha -\sigma \mu }\hfill \\ & +\frac{C}{n!2^{n-1}}{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left|log,\frac{T-\delta }{2}\right|}^{\mu -\nu }{\Vert v\Vert }_{\infty },\hfill \end{array}$$where $C>0$ is a general constant, $0<\kappa \le 1-\alpha -\sigma \mu$, and $\sigma >0$ satisfies $1-\alpha -\sigma \mu >0$.

Proof

We have shown that the error e(y) is estimated by (36). Now, we sequentially estimate the terms $J_{\rho }(y)$, $\rho =1,2,3,4$.

The term $J_1(y)$ is the Chebyshev interpolation error. Since $v(y)\in C^n[-1,1]$, we have the standard error estimate

54

$$\begin{array}{cc}& \Vert J_1{\Vert }_{\infty }\le \frac{1}{n!2^{n-1}}{\Vert v^{(n)}\Vert }_{\infty }.\hfill \end{array}$$We next estimate $J_4(y)$. Since ${\overline{\theta }}_{\nu }(y)\in C^{0,\kappa }[-1,1]$ with $0<\kappa \le 1-\alpha -\sigma \nu <1$, Lemma 3 implies that there exists a polynomial $p_{n-1}(y)$ of degree $n-1$, such that$$\begin{array}{cc}& \Vert {\overline{\theta }}_{\nu }-p_{n-1}{\Vert }_{\infty }\le C_{0,\kappa }{(n-1)}^{-\kappa }{\Vert {\overline{\theta }}_{\nu }\Vert }_{0,\kappa }.\hfill \end{array}$$Noting that ${\pi }_n{p}_{n-1}(y)=p_{n-1}(y)$, we further have by Lemma 4

55

$$\begin{array}{cc}\hfill \Vert {\pi }_n{\overline{\theta }}_{\nu }-{\overline{\theta }}_{\nu }{\Vert }_{\infty }& \le \Vert {\pi }_n{\overline{\theta }}_{\nu }-{\pi }_n{p}_{n-1}{\Vert }_{\infty }+{\Vert {\overline{\theta }}_{\nu }-p_{n-1}\Vert }_{\infty }\hfill \\ \hfill & \le (1+\Vert {\pi }_n{\Vert }_{\infty })\Vert {\overline{\theta }}_{\nu }-p_{n-1}{\Vert }_{\infty }\hfill \\ \hfill & \le C_{0,\kappa }{(n-1)}^{-\kappa }\left[\frac{2}{\pi },log,(n-1),+,1.9625,+,\mathcal{O},\left(\frac{1}{n}\right)\right]{\Vert {\overline{\theta }}_{\nu }\Vert }_{0,\kappa }.\hfill \end{array}$$Applying (41) to $Q(y,s)=K(y,s)e(s)$, we finally obtain

56

$$\begin{array}{cc}\hfill \Vert {\pi }_n{\overline{\theta }}_{\nu }-{\overline{\theta }}_{\nu }{\Vert }_{\infty }& \le C_{0,\kappa }{(n-1)}^{-\kappa }\left[\frac{2}{\pi },log,(n-1),+,1.9625,+,\mathcal{O},\left(\frac{1}{n}\right)\right]\hfill \\ \hfill & \times \left(C_1,\underset{y,s\in [-1,1]}{max},|K(y,s)|,+,C_2,\underset{y,s\in [-1,1]}{max},|K_y^{\prime }(y,s)|\right){\Vert e\Vert }_{\infty }.\hfill \end{array}$$Because ${(n-1)}^{-\kappa }log(n-1)\to 0$ as $n\to \infty$, we can conclude that there exists an integer $n_c$, such that when $n>n_c$

57

$$\begin{array}{cc}\hfill \Vert J_4{\Vert }_{\infty }& \le {\epsilon }_1{\Vert e\Vert }_{\infty },\hfill \end{array}$$where ${\epsilon }_1\le 1/4$ is a suitably small number.

For $J_2(y)$, we rewrite it as $J_2(y)=({\pi }_n{r}_{\delta }(y)-r_{\delta }(y))+r_{\delta }(y)$. The estimates (10) and (11) imply $\Vert r_{\delta }{\Vert }_{0,\kappa }\le {\overline{C}}_{\delta ,0}^{\ast }{\delta }^{{\alpha }_N}$ for $\forall \kappa \in (0,1-\alpha -\sigma \mu ]$. Further using the similar deduction for $J_4(y)$, we obtain

58

$$\begin{array}{cc}\hfill \Vert J_2{\Vert }_{\infty }& \le {\overline{C}}_{\delta ,0}{(n-1)}^{-\kappa }log(n-1){\delta }^{{\alpha }_N}+{\overline{C}}_{\delta ,1}{\delta }^{{\alpha }_N+1-\alpha -\sigma \mu }.\hfill \end{array}$$The estimate of $J_3(y)$ is also similar to $J_4(y)$. Noting that the Chebyshev interpolation of K(y, s) has the following error estimate$$\begin{array}{cc}& |K(y,s)-{\pi }_{n,s}K(y,s)|\le \frac{1}{n!2^{n-1}}\underset{s\in [-1,1]}{max}\left|\frac{{\partial }^{n}K(y,s)}{\partial s^n}\right|,\hfill \end{array}$$we have by (33), (39), (41), and (56)$$\begin{array}{cc}& \Vert {\widehat{\theta }}_{\nu }{\Vert }_{\infty }\le \frac{1}{n!2^{n-1}}\frac{2^{1-\alpha -\sigma \nu }}{{(\sigma \text{e})}^{\nu }(1-\alpha -\sigma \nu )}\underset{y,s\in [-1,1]}{max}\left|\frac{{\partial }^{n}K(y,s)}{\partial s^n}\right|{\Vert v_n\Vert }_{\infty },\hfill \\ \hfill & \Vert {\pi }_n{\widehat{\theta }}_{\nu }-{\widehat{\theta }}_{\nu }{\Vert }_{\infty }\le C_{0,\kappa }^{\ast }{(n-1)}^{-\kappa }log(n-1)\frac{1}{n!2^{n-1}}\hfill \\ \hfill & \times \left(C_1,\underset{y,s\in [-1,1]}{max},\left|\frac{{\partial }^{n}K(y,s)}{\partial s^n}\right|,+,C_2,\underset{y,s\in [-1,1]}{max},\left|\frac{{\partial }^{n+1}K(y,s)}{\partial y\partial s^n}\right|\right){\Vert v_n\Vert }_{\infty }.\hfill \end{array}$$Then, $J_3(y)$ can be estimated by

59

$$\begin{array}{cc}& \Vert J_3{\Vert }_{\infty }\le {\epsilon }_2{\Vert e\Vert }_{\infty }+{\left(\frac{T-\delta }{2}\right)}^{1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left|log,\frac{T-\delta }{2}\right|}^{\mu -\nu }\frac{C_{3,\nu }}{n!2^{n-1}}{\Vert v\Vert }_{\infty },\hfill \end{array}$$where ${\epsilon }_2\le 1/4$ is a suitably small number when $n>n_c$ and$$\begin{array}{cc}\hfill C_{3,\nu }=& \frac{2^{1-\alpha -\sigma \nu }}{{(\sigma \text{e})}^{\nu }(1-\alpha -\sigma \nu )}\underset{y,s\in [-1,1]}{max}\left|\frac{{\partial }^{n}K(y,s)}{\partial s^n}\right|+C_{0,\kappa }^{\ast }{(n-1)}^{-\kappa }log(n-1)\hfill \\ \hfill & \times \left(C_1,\underset{y,s\in [-1,1]}{max},\left|\frac{{\partial }^{n}K(y,s)}{\partial s^n}\right|,+,C_2,\underset{y,s\in [-1,1]}{max},\left|\frac{{\partial }^{n+1}K(y,s)}{\partial y\partial s^n}\right|\right).\hfill \end{array}$$Substituting (54), (57)–(59) into (36), and taking ${\epsilon }_1+{\epsilon }_2\le 1/2$, we can obtain (53). The proof is complete.

In Theorem 3, we only consider the transformed interval $[-1,1]$. By the variable change $x=\frac{T-\delta }{2}y+\frac{T+\delta }{2}$, we have the error estimate of SSCCM on $[\delta ,T]$.

Corollary 1

For the solution u(x) of (1), we have the error estimate when $x\in [\delta ,T]$

60

$$\begin{array}{cc}& \Vert u-u_n{\Vert }_{\infty }\le \frac{C}{n!2^{n-1}}{\left(\frac{T-\delta }{2}\right)}^n{\Vert u^{(n)}\Vert }_{\infty }+C{(n-1)}^{-\kappa }log(n-1){\delta }^{{\alpha }_N}\hfill \\ & +C{\delta }^{{\alpha }_N+1-\alpha -\sigma \mu }+\frac{C}{n!2^{n-1}}{\left(\frac{T-\delta }{2}\right)}^{n+1-\alpha }\sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left|log,\frac{T-\delta }{2}\right|}^{\mu -\nu }{\Vert u\Vert }_{\infty }.\hfill \end{array}$$

The error estimate (53) or (60) shows that the error of our SSCCM includes two main parts. One is generated by approximating the solution with its truncated psi-series $u_p(x)$ or its Padé approximation $u_r(x)$ on the singular interval $[0,\delta ]$, which is decreased by increasing ${\alpha }_N$ or decreasing $\delta$. The other is introduced by interpolating the solution using Chebyshev nodes on the regular interval $[\delta ,T]$. Noting that $u_r(x)$ is also an approximation to u(x) even when $x>\delta$, we can check the correctness of the collocation solution by comparing it with $u_r(x)$ at some points near $\delta$.

Remark 4

By Stirling formula $n!=\sqrt{2\pi n}{(n/\text{e})}^{n}exp({\theta }_n/(12n))$, ${\theta }_n\in (0,1)$, we can say that SSCCM has exponential convergence when $\delta$ is suitably chosen.

Extension to nonlinear equation

This section extends SSCCM to nonlinear WSVIE (2). Let $u_p(x)$ be the truncated psi-series for u(x) about $x=0$, which is accurate enough to u(x) on $[0,\delta ]$. As we have done in Section 2, we split (2) as

61

$$\begin{array}{cc}& u(x)=f(x)+G_{\delta }(x)+R_{\delta }(x)+{\int }_{\delta }^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}F(x,t,u(t))\text{d}t,x\in (\delta ,T],\hfill \end{array}$$where

62

$$\begin{array}{cc}& G_{\delta }(x)={\int }_0^{\delta }\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}F(x,t,u_p(t))\text{d}t,\hfill \end{array}$$

63

$$\begin{array}{cc}& R_{\delta }(x)={\int }_0^{\delta }\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\left[F,(x,t,u(t)),-,F,(x,t,u_p(t))\right]\text{d}t,x\in (\delta ,T].\hfill \end{array}$$Equation (61) can be directly discretized by interpolating F(x, t, u(t)) to t using the zeros of Chebyshev polynomials. Here, we use the Newton method to linearize (61) [48]. Drop away the error term $R_{\delta }(x)$ in (61), and let

64

$$\begin{array}{cc}& H(u)(x)=u(x)-f(x)-G_{\delta }(x)-{\int }_{\delta }^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}F(x,t,u(t))\text{d}t,\hfill \end{array}$$then the Newton method for $H(u)=0$ is

65

$$\begin{array}{cc}& {\overline{u}}_{l+1}={\overline{u}}_l-{[H^{\prime }({\overline{u}}_l)]}^{-1}H({\overline{u}}_l),l=0,1,\dots ,\hfill \end{array}$$where $H^{\prime }(u)$ is interpreted as Gâteaux derivative [48]. By denoting ${\varphi }_{l+1}={\overline{u}}_{l+1}-{\overline{u}}_l$, the iteration scheme (65) is equivalent to

66

$$\begin{array}{cc}& {\varphi }_{l+1}(x)-{\int }_{\delta }^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}\frac{\partial F(x,t,{\overline{u}}_l(t))}{\partial u}{\varphi }_{l+1}(t)\text{d}t=-H({\overline{u}}_l)(x).\hfill \end{array}$$The linear integral equation (66) can be effectively solved by SSCCM developed in Section 2. The initial guess ${\overline{u}}_0(x)$ can be chosen as the Padé approximation of the series expansion $u_p(x)$. The termination criterion for the iteration is $max|H({\overline{u}}_l)(x)|\le ϵ$, and $ϵ$ is the prescribed precision. Atkinson et al. [48] discussed the local convergence of the Newton method for nonlinear operator equations, but the given conditions are usually difficult to verify. Here, we choose the Padé approximation of $u_p(x)$ as a good initial function to guarantee that the Newton iteration converges as much as possible.

By our computational experience, the weakly singular integral in $H({\overline{u}}_l)(x)$ is better to be evaluated by the high precision integration method. Let $t=x-(x-\delta )s^{\beta }$, then the integral on the right-hand side of (64) is converted to

67

$$\begin{array}{cc}& {\int }_{\delta }^x\frac{{log}^{\mu }(x-t)}{{(x-t)}^{\alpha }}F(x,t,{\overline{u}}_l(t))\text{d}t={(x-\delta )}^{1-\alpha }\beta \sum _{\nu =0}^{\mu }\left(\begin{array}{c}\mu \\ \nu \end{array}\right){\left(log,(,x,-,\delta ,)\right)}^{\mu -\nu }{\beta }^{\nu }\hfill \\ & \times {\int }_0^1{s}^{\beta (1-\alpha )-1}{log}^{\nu }sG(x,s)\text{d}s,\hfill \end{array}$$where $G(x,s)=F(x,x-(x-\delta )s^{\beta },{\overline{u}}_l(x-(x-\delta )s^{\beta }))$. Here, we only consider the most used case that $\alpha =p/q$ ($p<q$) is a rational number, where p and q are coprime numbers. If we choose $\beta =q$, then the integral on the right-hand side of (67) is regular when $\nu =0$, which can be effectively evaluated by the standard Gauss-Legendre rule. As for the case $\nu =1$, we further perform a variable substitution $s={\text{e}}^{-y}$ to obtain

68

$$\begin{array}{cc}& {\int }_0^1{s}^{q-p-1}logsG(x,s)\text{d}s=-{\int }_0^{\infty }{\text{e}}^{-(q-p)y}yG(x,{\text{e}}^{-y})\text{d}y,\hfill \end{array}$$which can be effectively evaluated by the standard Gauss-Laguerre rule [47].

Numerical examples

In this section, we provide three examples, including linear and nonlinear equations, to illustrate the efficiency and accuracy of SSCCM developed in this paper. We note that all the experiments are implemented on a Laptop with Intel Core i5-6200U CPU (2.30 GHz) and 8 GB RAM by using Mathematica.

Example 1

Consider the following linear WSVIE of the second kind

69

$$\begin{array}{cc}& u(x)={\text{e}}^{-x}-{\int }_0^x\frac{{\text{e}}^{-2(x-t)}}{\sqrt{x-t}}u(t)\text{d}t,0\le x\le 10.\hfill \end{array}$$The exact solution is [49]$$\begin{array}{cc}\hfill u(x)& =\frac{1}{\pi -1}\left(\sqrt{\pi },{\text{e}}^{-x},\text{erf},\left(\sqrt{x}\right),+,\pi ,{\text{e}}^{(\pi -2)x},\text{erfc},\left(\sqrt{\pi x}\right),-,{\text{e}}^{-x}\right),\hfill \end{array}$$where $\text{erf}(x)$ and $\text{erfc}(x)$ are the error function and complementary error function, respectively.

For this example, we can obtain the truncated psi-series expansion for u(x) about $x=0$ by Picard iteration and its Padé approximation:$$\begin{array}{cc}\hfill u& {}_p(x)=1-2\sqrt{x}+2.14159x-1.52212x^{3/2}+0.722413x^2-0.428391x^{5/2}\hfill \\ \hfill & +0.441567x^3-0.292109x^{7/2}+0.0843558x^4-0.0190779x^{9/2}+\cdots \hfill \\ \hfill & +8.53383\times {10}^{-12}{x}^{15}-8.00654\times {10}^{-11}{x}^{31/2}+5.6109\times {10}^{-13}{x}^{16},\hfill \\ \hfill u& {}_r(x)=\frac{1-1.05861\sqrt{x}+1.13463x+\cdots +1.89725\times {10}^{-6}{x}^8}{1+0.941394\sqrt{x}+0.875828x+\cdots +2.00536\times {10}^{-9}{x}^8}.\hfill \end{array}$$Given $\epsilon ={10}^{-10}$ such that $|e_r(\delta )|\le \epsilon$ in (15), we obtain $\delta =2$, which is used to split the interval [0, 10]. We take $n=30$ to solve (69) by the Chebyshev collocation method on [0, 10] and [2, 10], respectively, and the absolute errors on the common logarithmic scale are plotted in Fig. 1, where E_CCM(n) ($\delta =0$) and E_SSCCM(n) ($\delta =2$) denote these errors using the partition number n. Here, we plot E_CCM(n) to compare the standard Chebyshev collocation method with our method. In this figure, we also plot the absolute errors of the truncated psi-series solution $u_p(x)$ and its Padé approximation $u_r(x)$, denoted by E_Psi and E_Pade, respectively.

Fig. 1 [Images not available. See PDF.]

The errors of approximate solutions on logarithmic scale in Example 1

It can be seen from Fig. 1 that the truncated psi-series $u_p(x)$ and its Padé approximation $u_r(x)$ are truly accurate enough to u(x) when x is small, but they become less accurate when x tends to be large. In addition, $u_r(x)$ is more accurate than $u_p(x)$. Hence, it is necessary to use SSCCM on the interval $[\delta ,T]$, and it is better to replace u(x) by $u_r(x)$ on $[0,\delta ]$. Figure 1 also shows that SSCCM gets a more accurate numerical solution than the standard Chebyshev collocation method (CCM, $\delta =0$), which means that it is necessary to separate the singularity for the equation with an insufficiently smooth solution.

We further list the maximal absolute errors ($E_{\delta }(n)$) and the implementation times in seconds (${\text{CPU}}_{\delta }(n)$) of the Chebyshev collocation methods for some n in Table 1, where in this table $2.60982\text{E}-5=2.60982\times {10}^{-5}$. As we expected, the SSCCM on the interval [2, 10] converges much faster than the standard CCM on [0, 10]. Table 1 shows that SSCCM and standard CCM cost less CPU time. Even for $n=30$, they cost less than one second. Certainly, SSCCM costs slightly much CPU time than the standard CCM since SSCCM needs to evaluate a singular integral on [0, 2]. By the way, we note that if we take the variable changes $y=\sqrt{x}$, $s=\sqrt{t}$ in (69), the solution is infinitely differentiable. Since the transformation leads $0\le x\le 10$ to $0\le y\le 100$, we must solve the transformed equation on a very large interval, which will cost much CPU time.

Table 1. The results ($E_{\delta }(n)$, ${\text{CPU}}_{\delta }(n)$) of CCMs in Example 1 ($\delta =0,2$)

n	20	22	25	28	30
$E_2(n)$	$2.60982\text{E}-5$	$8.52835\text{E}-7$	$3.52252\text{E}-9$	$1.18557\text{E}-10$	$4.48786\text{E}-11$
${\text{CPU}}_2(n)$	0.25	0.28125	0.328125	0.390625	0.46875
$E_0(n)$	0.231298	0.0187022	0.00289961	0.0024192	0.0021376
${\text{CPU}}_0(n)$	0.109375	0.140625	0.15625	0.234375	0.265625

Example 2

Consider the following linear WSVIE of the second kind with algebraic-logarithmic singular kernel

70

$$\begin{array}{cc}& u(x)={\text{e}}^{-\sqrt[3]{x}}+{\int }_0^x\frac{sin(x+t)log(x-t)}{\sqrt[3]{x-t}}u(t)\text{d}t,0\le x\le 2.\hfill \end{array}$$The exact solution is unknown.

The derivative of the forcing term in this example is weakly singular at $x=0$. For this example, we only list the truncated psi-series expansion for u(x) about $x=0$:$$\begin{array}{cc}& u_p(x)=1-x^{1/3}+\frac{1}{2}{x}^{2/3}-\frac{1}{6}x+\frac{1}{24}{x}^{4/3}+x^{5/3}(-4.14833+2.4logx)\hfill \\ \hfill & +x^2(3.91319-2.01533logx)+x^{7/3}(-1.86125+0.880024logx)+\cdots \hfill \\ \hfill & +x^{35/3}\left(-,1600.88,+,3782.5,log,x,-,3742.28,{log}^2,x,+,\cdots ,+,0.534098,{log}^7,x\right).\hfill \end{array}$$The Padé approximation $u_r(x)$ is omitted here. Since the analytical solution is unknown, we define the error functions of $u_p(x)$ and $u_r(x)$ by (15), as shown in Fig. 2 on common logarithmic scale. These errors increase quickly as x tends to be large. Hence, we only choose a small number as a splitting point. In the following computation, we take $\delta =0.1$ and denote the Chebyshev interpolation of the numerical solution by $u_n(x)$ on the interval $[\delta ,2]$. On this regular interval, we define the error function by$$\begin{array}{cc}\hfill e_n(x)=& u_n(x)-{\int }_0^{\delta }\frac{sin(x+t)log(x-t)}{\sqrt[3]{x-t}}{u}_r(t)\text{d}t\hfill \\ \hfill & -{\int }_{\delta }^x\frac{sin(x+t)log(x-t)}{\sqrt[3]{x-t}}{u}_n(t)\text{d}t.\hfill \end{array}$$We take $n=30$ in SSCCM, and the error is plotted in Fig. 2 (E_SSCCM(30)). As a comparison, we also plot the error of the standard CCM on the interval [0, 2] in this figure (E_CCM(30)). Obviously, SSCCM is much more accurate than the standard CCM.

Fig. 2 [Images not available. See PDF.]

The errors of approximate solutions on logarithmic scale in Example 2

Figure 2 shows that the results in this example are very similar to the ones in Example 1, which means that SSCCM is also successful in solving the problem with algebraic-logarithmic singular kernel and non-smooth forcing term.

Example 3

Solve the following nonlinear WSVIE of the second kind [6, 31, 39]

71

$$\begin{array}{cc}& u(x)=-\frac{1}{\sqrt{\pi }}{\int }_0^x\frac{1}{\sqrt{x-t}}{(u(t)-sint)}^3\text{d}t,0\le x\le 8,\hfill \end{array}$$which arises from a parabolic boundary value problem in the theory of superfluidity [6]. The exact solution is unknown, but we know $u(8)\approx 0.3236412904$ [50].

For this example, the truncated series expansion for u(x) about $x=0$ is$$\begin{array}{cc}\hfill u& {}_p(x)=0.51583x^{7/2}-0.208416x^{11/2}-0.61875x^6+0.0389044x^{15/2}\hfill \\ \hfill & +0.397182x^8+0.8972x^{17/2}-0.00434139x^{19/2}-0.123018x^{10}+\cdots \hfill \\ \hfill & +6.90231x^{37/2}-0.0824675x^{19}-0.934097x^{39/2}-5.00071x^{20}.\hfill \end{array}$$Its Padé approximation $u_r(x)$ has the following form$$\begin{array}{cc}\hfill u_r(x)& =\frac{0.51583x^{7/2}-0.34348x^4+0.0924853x^{9/2}+\cdots -0.116223x^{10}}{1-0.665878\sqrt{x}+0.179294x-0.394491x^{3/2}+\cdots -0.122341x^{10}}.\hfill \end{array}$$Noting that the coefficients of $u_p(x)$ are not decreasing as the exponents of x tend to be large, we can conclude that the series expansion for u(x) about $x=0$ is not convergent for large x. Even though, we can check that $u_r(0.4)$ has the precision of order ${10}^{-10}$ by (15). Hence, we take $\delta =0.4$ in this example. We use different n to solve this example by SSCCM on the interval [0.4, 8], and some errors of $E_{\delta }(n,8)=|u(8)-u_n(8)|$ are listed in Table 2. We also list the CPU time in seconds (${\text{CPU}}_{\delta }(n)$) in this table. The computational accuracy increases as n becomes large. For this nonlinear equation, the Newton method (65) is iterated about 8 times with stop criterion $ϵ={10}^{-10}$, and $u_r(x)$ is taken as the initial value. Considering the Newton iteration needs to solve a series of linear WSVIEs, it costs much CPU time compared with Table 1 for the linear equation.

Table 2. The results computed by SSCCM ($\delta =0.4$) in Example 3

n	20	30	40	45
$E_{\delta }(n,8)$	$6.45824\text{E}-4$	$7.22768\text{E}-5$	$7.12443\text{E}-6$	$3.22921\text{E}-7$
${\text{CPU}}_{\delta }(n)$	1.89063	4.04688	8.14063	10.9531

Conclusion

In this paper, we design a singularity separation Chebyshev collocation method to solve second kind linear and nonlinear Volterra integral equations with algebraic and/or logarithmic singular kernels and non-smooth forcing terms. Since the asymptotic expansion of the solution about the origin can be obtained by Picard iteration, which is a good approximation to the solution near the origin, we split the interval into a singular subinterval and a regular one. In the singular subinterval, the truncated asymptotic expansion or its Padé approximation is used to approximate the solution directly. In the remaining subinterval, the solution is regular even though the equation is still singular. The Chebyshev collocation method is used to discretize the equation. As for the singular kernel, we design a stable and fast recurrence algorithm to evaluate the singular weight integrals analytically. Convergence analysis shows that the error of our SSCCM includes two main parts, one is generated by approximating the solution with its truncated psi-series on the singular subinterval and the other is introduced by interpolating the solution using Chebyshev nodes on the regular subinterval. We can balance the two errors in practical computation by choosing an appropriate splitting point for the interval.

Acknowledgements

The authors are very grateful to Editors and Referees for the valuable comments, which improve the quality of the paper significantly.

Author Contributions

All authors contributed to the study conception and design. The manuscript was written by Tongke Wang. All authors read and approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China under Grant No. 11971241 and the Program for Innovative Research Team in Universities of Tianjin under Grant No. TD13-5078.

Data availability

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Declarations

Conflict of interest

The authors declare no competing interests.

Ethical approval and consent to participate

Not applicable.

Consent for publication

All of the material is owned by the authors, and no permissions are required.

Human and animal ethics

Not applicable.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Gorenflo, R; Vessella, S. Abel Integral Equations; 1991; Berlin, Springer: [DOI: https://dx.doi.org/10.1007/BFb0084665]

2. Brunner, H. Volterra Integral Equations: An Introduction to Theory and Applications; 2017; Cambridge, Cambridge University Press: [DOI: https://dx.doi.org/10.1017/9781316162491]

3. Lighthill, MJ. Contributions to the theory of heat transfer through a laminar boundary layer. Proc. R. Soc. Lond. A; 1950; 202, pp. 359-377.1950RSPSA.202.359LMathSciNet ID: 37147[COI: 1:CAS:528:DyaG3MXis1Wmsg%3D%3D] [DOI: https://dx.doi.org/10.1098/rspa.1950.0106]

4. Diogo, T., Ma, J.T., Rebelo, M.: Fully discretized collocation methods for nonlinear singular Volterra integral equations. J. Comput. Appl. Math. 247,84–101 (2013). https://doi.org/10.1016/j.cam.2013.01.002

5. Huang, C., Wang, L.L.: An accurate spectral method for the transverse magnetic mode of Maxwell equations in Cole-Cole dispersive media. Adv. Comput. Math. 45,707–734 (2019). https://doi.org/10.1007/s10444-018-9636-2

6. Levinson, N. A nonlinear Volterra equation arising in the theory of superfluidity. J. Math. Anal. Appl.; 1960; 1, pp. 1-11. [DOI: https://dx.doi.org/10.1016/0022-247X(60)90028-7]

7. Brunner, H. Collocation Methods for Volterra Integral and Related Functional Equations; 2004; Cambridge, Cambridge University Press: [DOI: https://dx.doi.org/10.1017/CBO9780511543234]

8. Cao, Y.Z., Herdman, T., Xu, Y.S.: A hybrid collocation method for Volterra integral equations with weakly singular kernels. SIAM J. Numer. Anal. 41, 364–381 (2004). https://doi.org/10.1137/S0036142901385593

9. Brunner, H. The numerical solution of weakly singular Volterra integral equations by collocation on graded meshes. Math. Comp.; 1985; 45, pp. 417-437.MathSciNet ID: 804933[DOI: https://dx.doi.org/10.1090/S0025-5718-1985-0804933-3]

10. Khater, AH; Shamardan, AB; Callebaut, DK; Sakran, MRA. Solving integral equations with logarithmic kernels by Chebyshev polynomials. Numer. Algorithms; 2008; 47, pp. 81-93.2008NuAlg.47..81KMathSciNet ID: 2372893[DOI: https://dx.doi.org/10.1007/s11075-007-9148-5]

11. Liang, H; Brunner, H. The convergence of collocation solutions in continuous piecewise polynomial spaces for weakly singular Volterra integral equations. SIAM J. Numer. Anal.; 2019; 57, pp. 1875-1896.MathSciNet ID: 3986560[DOI: https://dx.doi.org/10.1137/19M1245062]

12. Kant, K; Nelakanti, G. Approximation methods for second kind weakly singular Volterra integral equations. J. Comput. Appl. Math.; 2020; 368,MathSciNet ID: 4020569[DOI: https://dx.doi.org/10.1016/j.cam.2019.112531]

13. Brunner, H. Nonpolynomial spline collocation for Volterra equations with weakly singular kernels. SIAM J. Numer. Anal.; 1983; 20, pp. 1106-1119.1983SJNA..20.1106BMathSciNet ID: 723827[DOI: https://dx.doi.org/10.1137/0720080]

14. Rebelo, M; Diogo, T. A hybrid collocation method for a nonlinear Volterra integral equation with weakly singular kernel. J. Comput. Appl. Math.; 2010; 234, pp. 2859-2869.MathSciNet ID: 2652131[DOI: https://dx.doi.org/10.1016/j.cam.2010.01.034]

15. Ma, JJ; Liu, HL. Fractional collocation boundary value methods for the second kind Volterra equations with weakly singular kernels. Numer. Algorithms; 2020; 84, pp. 743-760.MathSciNet ID: 4098865[COI: 1:CAS:528:DC%2BB3cXovFGr] [DOI: https://dx.doi.org/10.1007/s11075-019-00777-9]

16. Huang, QM; Wang, M. Superconvergence of interpolated collocation solutions for weakly singular Volterra integral equations of the second kind. Comput. Appl. Math.; 2021; 40, 71.MathSciNet ID: 4226463[DOI: https://dx.doi.org/10.1007/s40314-021-01435-4]

17. Eshaghi, J; Adibi, H; Kazem, S. Solution of nonlinear weakly singular Volterra integral equations using the fractional-order Legendre functions and pseudospectral method. Math. Meth. Appl. Sci.; 2016; 39, pp. 3411-3425.MathSciNet ID: 3521264[DOI: https://dx.doi.org/10.1002/mma.3788]

18. Cai, HT; Chen, YP. A fractional order collocation method for second kind Volterra integral equations with weakly singular kernels. J. Sci. Comput.; 2018; 75, pp. 970-992.MathSciNet ID: 3780795[DOI: https://dx.doi.org/10.1007/s10915-017-0568-7]

19. Cai, HT. A fractional spectral collocation for solving second kind nonlinear Volterra integral equations with weakly singular kernels. J. Sci. Comput.; 2019; 80, pp. 1529-1548.MathSciNet ID: 3995976[DOI: https://dx.doi.org/10.1007/s10915-019-00987-2]

20. Hou, DM; Lin, YM; Azaiez, M; Xu, CJ. A Müntz-collocation spectral method for weakly singular Volterra integral equations. J. Sci. Comput.; 2019; 81, pp. 2162-2187.MathSciNet ID: 4036983[DOI: https://dx.doi.org/10.1007/s10915-019-01078-y]

21. Talaei, Y. Chelyshkov collocation approach for solving linear weakly singular Volterra integral equations. J. Appl. Math. Comput.; 2019; 60, pp. 201-222.MathSciNet ID: 3969080[DOI: https://dx.doi.org/10.1007/s12190-018-1209-5]

22. Pedas, A; Vainikko, G. Smoothing transformation and piecewise polynomial collocation for weakly singular Volterra integral equations. Computing; 2004; 73, pp. 271-293.MathSciNet ID: 2106251[DOI: https://dx.doi.org/10.1007/s00607-004-0088-9]

23. Zhao, JJ; Long, T; Xu, Y. Super implicit multistep collocation methods for weakly singular Volterra integral equations. Numer. Math. Theor. Meth. Appl.; 2019; 12, pp. 1039-1065.MathSciNet ID: 3978738[DOI: https://dx.doi.org/10.4208/nmtma.OA-2018-0084]

24. Chen, YP; Tang, T. Convergence analysis of the Jacobi spectral-collocation methods for Volterra integral equations with a weakly singular kernel. Math. Comp.; 2010; 79, pp. 147-167.2010MaCom.79.147CMathSciNet ID: 2552221[DOI: https://dx.doi.org/10.1090/S0025-5718-09-02269-8]

25. Liu, X; Chen, YP. Convergence analysis for the Chebyshev collocation methods to Volterra integral equations with a weakly singular kernel. Adv. Appl. Math. Mech.; 2017; 9, pp. 1506-1524.MathSciNet ID: 3748995[DOI: https://dx.doi.org/10.4208/aamm.OA-2016-0049]

26. Allaei, SS; Diogo, T; Rebelo, M. The Jacobi collocation method for a class of nonlinear Volterra integral equations with weakly singular kernel. J. Sci. Comput.; 2016; 69, pp. 673-695.MathSciNet ID: 3551341[DOI: https://dx.doi.org/10.1007/s10915-016-0213-x]

27. Sohrabi, S; Ranjbar, H; Saei, M. Convergence analysis of the Jacobi-collocation method for nonlinear weakly singular Volterra integral equations. Appl. Math. Comput.; 2017; 299, pp. 141-152.MathSciNet ID: 3588581[DOI: https://dx.doi.org/10.1016/j.amc.2016.11.022]

28. Li, XJ; Tang, T; Xu, CJ. Numerical solutions for weakly singular Volterra integral equations using Chebyshev and Legendre pseudo-spectral Galerkin methods. J. Sci. Comput.; 2016; 67, pp. 43-64.MathSciNet ID: 3473692[DOI: https://dx.doi.org/10.1007/s10915-015-0069-5]

29. Baratella, P; Orsi, AP. A new approach to the numerical solution of weakly singular Volterra integral equations. J. Comput. Appl. Math.; 2004; 163, pp. 401-418.2004JCoAm.163.401BMathSciNet ID: 2042821[DOI: https://dx.doi.org/10.1016/j.cam.2003.08.047]

30. Baratella, P. A Nyström interpolant for some weakly singular nonlinear Volterra integral equations. J. Comput. Appl. Math.; 2013; 237, pp. 542-555.MathSciNet ID: 2966927[DOI: https://dx.doi.org/10.1016/j.cam.2012.06.024]

31. Lü, T; Huang, Y. Extrapolation method for solving weakly singular nonlinear Volterra integral equations of the second kind. J. Math. Anal. Appl.; 2006; 324, pp. 225-237.MathSciNet ID: 2262467[DOI: https://dx.doi.org/10.1016/j.jmaa.2005.12.013]

32. Wang, TK; Qin, M; Zhang, ZY. The Puiseux expansion and numerical integration to nonlinear weakly singular Volterra integral equations of the second kind. J. Sci. Comput.; 2020; 82, 64.MathSciNet ID: 4067294[COI: 1:CAS:528:DC%2BB3MXpvFyhs7o%3D] [DOI: https://dx.doi.org/10.1007/s10915-020-01167-3]

33. Wang, TK; Gu, YS; Zhang, ZY. An algorithm for the inversion of Laplace transforms using Puiseux expansions. Numer. Algorithms; 2018; 78, pp. 107-132.MathSciNet ID: 3788266[DOI: https://dx.doi.org/10.1007/s11075-017-0369-y]

34. Wang, TK; Qin, M; Lian, H. The asymptotic approximations to linear weakly singular Volterra integral equations via Laplace transform. Numer. Algorithms; 2020; 85, pp. 683-711.MathSciNet ID: 4146530[DOI: https://dx.doi.org/10.1007/s11075-019-00832-5]

35. Hemmi, MA; Melkonian, S. Convergence of psi-series solutions of nonlinear ordinary differential equations. Canad. Appl. Math. Q.; 1995; 3, pp. 43-88.MathSciNet ID: 1341920

36. Wang, TK; Liu, ZF; Kong, YT. The series expansion and Chebyshev collocation method for nonlinear singular two-point boundary value problems. J. Eng. Math.; 2021; 126, 5.MathSciNet ID: 4206633[DOI: https://dx.doi.org/10.1007/s10665-020-10077-0]

37. Tang, T; Xu, X; Cheng, J. On spectral methods for Volterra integral equations and the convergence analysis. J. Comput. Math.; 2008; 26, pp. 825-837.MathSciNet ID: 2464738

38. Shen, J; Tang, T; Wang, LL. Spectral Methods: Algorithms; 2011; New York, Analysis and Applications. Springer: [DOI: https://dx.doi.org/10.1007/978-3-540-71041-7]

39. Orsi, AP. Product integration for Volterra integral equations of the second kind with weakly singular kernels. Math. Comp.; 1996; 65, pp. 1201-1212.1996MaCom.65.1201PMathSciNet ID: 1344619[DOI: https://dx.doi.org/10.1090/S0025-5718-96-00736-3]

40. Allaei, SS; Diogo, T; Rebelo, M. Analytical and computational methods for a class of nonlinear singular integral equations. Appl. Numer. Math.; 2017; 114, pp. 2-17.MathSciNet ID: 3601626[DOI: https://dx.doi.org/10.1016/j.apnum.2016.06.001]

41. Huang, C; Stynes, M. A spectral collocation method for a weakly singular Volterra integral equation of the second kind. Adv. Comput. Math.; 2016; 42, pp. 1015-1030.MathSciNet ID: 3554362[DOI: https://dx.doi.org/10.1007/s10444-016-9451-6]

42. Wang, TK; Liu, ZF; Zhang, ZY. The modified composite Gauss type rules for singular integrals using Puiseux expansions. Math. Comp.; 2017; 86, pp. 345-373.MathSciNet ID: 3557802[DOI: https://dx.doi.org/10.1090/mcom/3105]

43. Wang, TK; Zhang, ZY; Liu, ZF. The practical Gauss type rules for Hadamard finite-part integrals using Puiseux expansions. Adv. Comput. Math.; 2017; 43, pp. 319-350.MathSciNet ID: 3627469[DOI: https://dx.doi.org/10.1007/s10444-016-9487-7]

44. Brezinski, C; Van Iseghem, J. A taste of Padé approximation. Acta Numerica; 1995; 4, pp. 53-103.1995AcNum..4..53B[DOI: https://dx.doi.org/10.1017/S096249290000252X]

45. Mason, JC; Handscomb, D. Chebyshev Polynomials; 2003; Boca Raton, Chapman & Hall/CRC:

46. Li, CP; Zeng, FR; Liu, FW. Spectral approximations to the fractional integral and derivative. Fract. Calc. Appl. Anal.; 2012; 15, pp. 383-406.MathSciNet ID: 2944106[DOI: https://dx.doi.org/10.2478/s13540-012-0028-x]

47. Davis, PJ; Rabinowitz, P. Methods of Numerical Integration; 1984; 2 New York, Academic Press:

48. Atkinson, K; Han, WM. Theoretical Numerical Analysis: A Functional Analysis Framework; 2009; New York, Springer:

49. Polyanin, AD; Manzhirov, AV. Handbook of Integral Equations; 2008; 2 London, Chapman & Hall/CRC: [DOI: https://dx.doi.org/10.1201/9781420010558]

50. Lubich, C. Fractional linear multistep methods for Abel-Volterra integral equations of the second kind. Math. Comp.; 1985; 45, pp. 463-469.MathSciNet ID: 804935[DOI: https://dx.doi.org/10.1090/S0025-5718-1985-0804935-7]