Sinc Collocation Method to Simulate the Fractional Partial Integro-Differential Equation with a Weakly Singular Kernel

Abstract

In this article, we develop an efficient numerical scheme for dealing with fractional partial integro-differential equations (FPIEs) with a weakly singular kernel. The weight and shift Grünwald difference (WSGD) operator is adopted to approximate a time fractional derivative and the Sinc collocation method is applied for discretizing the spatial derivative.The exponential convergence of our proposed method is demonstrated in detail. Finally, numerical evidence is employed to verify the theoretical results and confirm the expected convergence rate.

1. Introduction

In recent years, fractional calculus has played an increasingly important role in various fields and has attracted much interest from scholars due to its extensive applications in modeling many complex problems [1,2,3,4,5,6,7,8]. The fractional integro-differential equation is one of the most active fields in fractional calculus [9,10,11,12,13,14,15], which can be seen as the extension of classical integral equations by replacing integer-order derivatives with fractional derivatives. Consequently, there is a growing need to explore solution techniques to study these equations. Although there are some ways to get exact solutions of these equations, the exact solutions of these equations are very difficult to find in most cases. For this reason, there has been much research on the effective numerical methods of fractional integral differential equations (FIDEs). Here, we list only a few of them. In [16,17], the homotopy analysis method is used to find the approximate solution of FIDEs. In [18], the spectral Jacobi-collocation method is presented by Ma et al. to solve the solution of general linear FIDEs. In [19], the compact finite difference scheme is constructed to approximate the solution of FIDEs with a weakly singular kernel. In [20], the alternating direction implicit difference scheme combined with a fractional trapezoidal rule is developed to solve two dimensional FIDEs. In [21], the Legendre wavelet collocation method based on the Gauss–Jacobi quadrature is introduced to solve the fractional delay-type integro-differential equations. In [22], the collocation method combined with fractional Genocchi functions is used for the solution of variable-order FIDEs. In [23], the meshless method based on the Laplace transform is constructed for approximating the solution of the two-dimensional multi-term FIDEs. In [24], the finite element method is proposed to solve the two-dimensional weakly singular FIDEs. In [25], the spectral Galerkin method based on Legendre polynomials is presented to solve the one and two-dimensional fourth-order FIDEs. In [26], the Adomian decomposition method and homotopy perturbation method are given to approximate the solution of the time FPIEs. Liu et al. [27,28,29,30,31] used the multigrid method and homotopy method to solve practical problems in the fractional flow formulation of the two-phase porous media flow equations and Biot elastic models.

In this article, we consider the fractional partial integro-differential equations (FPIEs) with a weakly singular kernel as follows:

$${}_0{D}_t^{\alpha }v(x, t)={\displaystyle \frac{{\partial }^{2}v(x, t)}{\partial x^2}}+{\int }_0^t{(t-s)}^{-1/2}{\displaystyle \frac{{\partial }^{2}v(x, t)}{\partial x^2}}ds+f(x, t),x\in [a, b],t\in [0, T],$$

(1)

with the initial condition:

$$\begin{array}{c}v(x, 0)=0,x\in [a, b],\hfill \end{array}$$

(2)

and boundary conditions

$$\begin{array}{c}v(a, t)=v(b, t)=0,t\in [0, T],\hfill \end{array}$$

(3)

where $0<\alpha <1$ and $f(x, t)$ is smooth enough. The time fractional derivative ${}_0{D}_t^{\alpha }v(x, t)$ is defined in Riemann–Liouville sense as

$${}_0{D}_t^{\alpha }v(x, t)=\frac{1}{\Gamma (1-\alpha )}\frac{d}{dt}{\int }_0^t\frac{v{(x, s)}^{\alpha }}{(t-s)}ds,$$

where $\Gamma (·)$ is the Gamma function.

The partial integro-differential equation of integer order has proven to describe some phenomena such as viscoelasticity, population dynamics and heat conduction in materials with memory [32,33,34]. Over the past three decades, various numerical methods based on the Sinc approximation have been presented, which have the advantages of a very fast convergence of exponential order and handling singularities effectively. The Sinc method proposed by Frank Stenger [35,36,37] has been increasingly applied to solve a variety of linear and nonlinear models that arise in scientific and engineering applications such as two-point boundary value problems [38], the Blasius equation [39], oceanographic problems with boundary layers [40], fourth-order partial integro-differential equation [41], the Volterra integro-differential equation [42,43], optimal control, heat distribution and astrophysics equations. According to the definition, it can be seen that fractional derivatives and integrals always deal with weak singularities. Therefore, in the past few years, the Sinc method has been widely extended to get the numerical solution of the fractional differential equations [44,45,46,47]. The main objective of this work is to provide a new attempt to develop a numerical solution via the use of the Sinc collocation method to solve the fractional partial integro-differential equation with a weakly singular kernel.

The remainder of this article is organized as follows. In Section 2, we introduce some basic formulation and theoretic results of Sinc functions which are required for our subsequent development. In Section 3, we propose a time discrete scheme based on the weight and shift Grünwald difference operator and a space discrete scheme by applying a collocation scheme based on the Sinc functions. In Section 4, the convergence analysis of our scheme is proved and in the meantime the exponential convergence is obtained. Some numerical results are described in Section 5 to illustrate the performance of our method. Finally, we give our conclusion in Section 6.

2. Definitions and Preliminaries

In this section, we describe some main notations and definitions of the Sinc function and review some known results that will be used in the following sections. The reader interested in learning more about the detailed properties of the Sinc function can investigate [35].

The Sinc function is basically defined on the whole real line $-\infty <x<\infty$ by

$$Sinc\left(x\right)=\left\{\begin{array}{cc}\frac{sin\left(\pi x\right)}{\pi x},\hfill & \hfill x\ne 0,\\ 1,\hfill & \hfill x=0.\end{array}\right.$$

For any mesh size $h>0$ and $k=0,\pm 1,\pm 2,\dots$, the Sinc basis functions with evenly spaced nodes given on $\mathbb{R}$ by

$$S(k,h)\left(x\right)=Sinc\left(\frac{y-kh}{h}\right).$$

Let $f\left(x\right)$ be a function defined on the real line, then for $h>0$ the series

$$C(f,h)\left(x\right)=\sum _{k=-\infty }^{\infty }f\left(kh\right)S(k,h)\left(x\right),$$

is called the Whittaker cardinal expansion of f, whenever this series converges. The properties of Whittaker cardinal expansions which are derived on the infinite strip $Q_s$ of the complex plane have been described and proved in detail in the literature [36].

$$Q_s=\left\{w=t+is:\left|s\right|<d\le \frac{\pi }{2}\right\}.$$

In order to construct approximations on the interval (a, b), we consider the conformal map

$$\varphi \left(z\right)=log\left(\frac{z-a}{b-z}\right),$$

which carries the eye-shaped domain of complex plan

$$Q_E=\left\{z=x+iy:\left|arg\left(\frac{z-a}{b-z}\right)\right|<d\le \frac{\pi }{2}\right\},$$

onto the infinite strip domain of complex plan $Q_s$.

Let $\psi$ be the inverse map of $w=\varphi \left(z\right)$, and we define the range of ϕ⁻¹ on $\mathbb{R}$ as

$$(a,b)=\{\psi \left(u\right)={\varphi }^{-1}\left(u\right)\in Q_E:-\infty <u<\infty \}.$$

For the uniform grid ${\left\{jh\right\}}_{j=-\infty }^{\infty }$ on $\mathbb{R}$, the Sinc points which correspond to these nodes are denoted by

$$x_j=\psi \left(jh\right)=\frac{a+b{e}^{jh}}{1+e^{jh}},j=0,\pm 1,\pm 2,\dots .$$

(4)

The basis functions on (a, b) for $z\in Q_E$ are taken to be the composite translated Sinc functions as

$$S_k\left(z\right)=S(k,h)\circ \varphi \left(z\right)=Sinc\left(\frac{\varphi \left(z\right)-kh}{h}\right),k=0,\pm 1,\pm 2,\dots .$$

(5)

Definition 1.

Let $B\left(Q_E\right)$ be the class of functions F which are analytic in $Q_E$ and satisfy

$${\int }_{\psi (t+\Sigma )}\left|F\left(z\right)\right|\mathrm{d}z\to 0,ast\to \pm \infty ,$$

where $\Sigma =\left\{\mathrm{i}\eta :\left|\eta \right|<d\le {\displaystyle \frac{\pi }{2}}\right\}$ and satisfy

$$N\left(f\right)={\int }_{\partial Q_E}\left|F\left(z\right)\mathrm{d}z\right|<\infty .$$

where $\partial Q_E$ represents the boundary of $Q_E$.

Lemma 1

([37]). If $F\left(x\right),{\varphi }^{\prime }\left(x\right)\in B\left(Q_E\right)$ and $h>0.$ Let ϕ be a one-to-one comformal map. For all $x\in (a,b)$

$$|F\left(x\right)-\sum _{j=-N}^{N}F\left(x_j\right)S_j\left(x\right)|\le {\displaystyle \frac{2N\left(f{\varphi }^{\prime }\right)}{\pi d}}{\mathrm{e}}^{-\pi d/h}.$$

(6)

Moreover, if $\left|F\left(x\right)\right|\le c{e}^{-\beta \left|\varphi \right(x\left)\right|}$, $x\in (a,b)$, for some positive constants c, l, N and β, and if the selection $h=\sqrt{\pi d/\beta N}$, then

$$\underset{x\in (0,1)}{sup}|F\left(x\right)-{\left({\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\right)}^l\sum _{j=-N}^{N}f\left(x_j\right)S_j\left(x\right)|\le C_1{N}^{(l+1)/2}{\mathrm{e}}^{-\sqrt{\pi d\beta N}},$$

(7)

where $C_1$ depends only on F, d and β.

Lemma 1 shows the Sinc interpolation on $B\left(Q_E\right)$ is exponentially convergent. We also need the derivative of the complex Sinc function to be evaluated at the nodes. So we introduce the lemma as follows:

Lemma 2

([48]). For the step size h and the Sinc points $x_j$ determined by (4), suppose that ϕ is the conformal one-to-one mapping of the simply connected domain $Q_E$ onto $Q_d$, then we have

$${\delta }_{kj}^{\left(0\right)}=[S(k,h)\circ \varphi \left(x\right)]{|}_{x=x_j}=\left\{\begin{array}{c}1,j=k,\hfill \\ 0,j\ne k,\hfill \end{array}\right.$$

(8)

$${\delta }_{kj}^{\left(1\right)}=h\frac{d}{d\varphi }[S(k,h)\circ \varphi \left(x\right)]{|}_{x=x_j}=\left\{\begin{array}{c}0,j=k,\hfill \\ \frac{{(-1)}^{j-k}}{j-k},j\ne k,\hfill \end{array}\right.$$

(9)

$${\delta }_{kj}^{\left(2\right)}=h^2\frac{d^2}{d{\varphi }^2}[S(k,h)\circ \varphi \left(x\right)]{|}_{x=x_j}=\left\{\begin{array}{c}\frac{-{\pi }^2}{3},j=k,\hfill \\ \frac{-2{(-1)}^{j-k}}{{(j-k)}^2},j\ne k.\hfill \end{array}\right.$$

(10)

To facilitate the representation of discrete systems, we give the definition of the following matrix as follows:

$$I^{\left(l\right)}=\left[{\delta }_{kj}^{\left(l\right)}\right],l=0,1,2,$$

(11)

where ${\delta }_{kj}^{\left(l\right)}$ is the (k, j) the element of the matrix $I^{\left(l\right)}$. The matrix $I^{\left(0\right)}$, $I^{\left(1\right)}$ and $I^{\left(2\right)}$ represents the identity matrix, the skew symmetric Toeplitz matrix and the symmetric Toeplitz matrix, respectively.

3. Derivation of the Numerical Scheme

3.1. The Time Semi-Discretization

For positive integer number N, let $\tau =\frac{T}{N}$ be the time mesh size, $t_n=n\tau$, $n=0,1,\cdots ,N,$ be the mesh points. Denote $v^n=v(x,t_n)$ and $f^n=f(x,t_n)$. Firstly, in order to apply the WSGD operator to discrete the time fractionl derivatives, the approximation order must be used. Therefore, we review the following lemma.

Lemma 3

([49]). Suppose that $\phi \left(t\right)\in L_1\left(\mathbb{R}\right)$ and $\phi \left(t\right)\in C^{\alpha +1}\left(\mathbb{R}\right),$ and define the shift Grünwald difference operator by

$$A_{\tau ,p}^{\alpha }\phi \left(t\right)=\frac{1}{{\tau }^{\alpha }}\sum _{k=0}^{\infty }{g}_k^{\left(\alpha \right)}\phi (t-(k-p)\tau ),$$

(12)

where p is an integer and the sequences $g_k^{\left(\alpha \right)}$ are the coefficients of the power series expansion of the function ${(1-z)}^{\alpha }$, i.e, $g_0^{\left(\alpha \right)}=1$, $g_k^{\left(\alpha \right)}={(-1)}^k\left(\genfrac{}{}{0pt}{}{\alpha }{k}\right)$, $k=1,2,\cdots$. Then

$$A_{\tau ,p}^{\alpha }\phi \left(t\right)={{}_{-\infty }D}_t^{\alpha }\phi \left(t\right)+O\left(\tau \right),$$

(13)

uniformly for $t\in \mathbb{R}$ as $\tau \to 0.$

Lemma 4

([50]). Let $\phi \left(t\right)\in L_1\left(\mathbb{R}\right)$, ${}_{-\infty }{D}_t^{\alpha +2}\phi \left(t\right)$ and its Fourier transform belong to $L_1\left(\mathbb{R}\right)$. Define the weighted and shifted Grünwald difference operator by

$$\begin{array}{c}{D}_{\tau ,p,q}^{\alpha }\phi \left(t\right)={\displaystyle \frac{2q-\alpha }{2(q-p)}}{A}_{\tau ,p}^{\alpha }\phi \left(t\right)+{\displaystyle \frac{2p-\alpha }{2(p-q)}}{A}_{\tau ,q}^{\alpha }\phi \left(t\right),\hfill \end{array}$$

(14)

where p and q are integers and $p\ne q$. Then,

$$D_{\tau ,p,q}^{\alpha }\phi \left(t\right)={{}_{-\infty }D}_t^{\alpha }\phi \left(t\right)+O\left({\tau }^2\right),$$

uniformly for $t\in \mathbb{R}$ as $\tau \to 0.$

Let p and q be equal to 0 and 1 in Lemma 4, we can get

$$\begin{array}{c}{}_0{D}_t^{\alpha }v(x_i,t_n)={\tau }^{-\alpha }\left({\displaystyle \frac{2+\alpha }{2}}{\displaystyle \sum _{k=0}^n}{g}_k^{\left(\alpha \right)}{v}_i^{n-k}-{\displaystyle \frac{\alpha }{2}}{\displaystyle \sum _{k=0}^{n-1}}{g}_k^{\left(\alpha \right)}{v}_i^{n-1-k}\right)+O\left({\tau }^2\right)\hfill \\ ={\tau }^{-\alpha }{\displaystyle \sum _{k=0}^n}{\lambda }_k{v}_i^{n-k}+R_{n+1,1},\hfill \end{array}$$

(15)

where $R_{n+1,1}=O\left({\tau }^2\right),{\lambda }_0={\displaystyle \frac{2+\alpha }{2}}{g}_0^{\left(\alpha \right)},{\lambda }_k={\displaystyle \frac{2+\alpha }{2}}{g}_k^{\left(\alpha \right)}-{\displaystyle \frac{1-\alpha }{2}}{g}_{k-1}^{\left(\alpha \right)},k\ge 1.$

By using unusual quadrature approximation, the integral term of (1) can be approximated as follows:

$$\begin{array}{c}{\displaystyle {\int }_0^{t_{n+1}}{(t_{n+1}-s)}^{-1/2}{\displaystyle \frac{{\partial }^{2}v(x,s)}{\partial x^2}}\mathrm{d}s}\hfill \\ {\displaystyle ={\displaystyle \sum _{l=0}^n}{\int }_{t_l}^{t_{l+1}}{(t_{n+1}-s)}^{-1/2}{\displaystyle \frac{{\partial }^{2}v(x,s)}{\partial x^2}}\mathrm{d}s}\hfill \\ {\displaystyle \approx {\displaystyle \sum _{l=0}^n}{\int }_{t_l}^{t_{l+1}}{(t_{n+1}-s)}^{-1/2}\left({\displaystyle \frac{t_{l+1}-s}{\tau }}{v}_{xx}^l\left(x\right)+{\displaystyle \frac{s-t_l}{\tau }}{v}_{xx}^{l+1}\left(x\right)\right)\mathrm{d}s}\hfill \\ {\displaystyle \approx {\displaystyle \frac{1}{\tau }}{\displaystyle \sum _{l=0}^n}\left(A_{n,l}{v}_{xx}^l\left(x\right)+B_{n,l}{v}_{xx}^{l+1}\left(x\right)\right)+R_{n+1,2},}\hfill \end{array}$$

(16)

where

$$\begin{array}{c}{R}_{n+1,2}=O(▵t^{3/2}),\hfill \\ {\displaystyle A_{n,l}={\int }_{t_l}^{t_{l+1}}{(t_{n+1}-s)}^{-1/2}(t_{l+1}-s)\mathrm{d}s,}\hfill \\ {\displaystyle B_{n,l}={\int }_{t_l}^{t_{l+1}}{(t_{n+1}-s)}^{-1/2}(s-t_l)\mathrm{d}s,}\hfill \end{array}$$

(17)

Substituting (15) and (16) into (1), we have

$$\begin{array}{c}{\lambda }_0{v}^{n+1}\left(x\right)-\left({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n}\right)v_{xx}^{n+1}\left(x\right)\hfill \\ {\displaystyle ={\tau }^{\alpha -1}{\displaystyle \sum _{l=0}^n}{\rho }_{n,l}{v}_{xx}^l\left(x\right)-{\displaystyle \sum _{k=0}^n}{\lambda }_k{v}^{n+1-k}\left(x\right)+{\tau }^{\alpha }{f}^{n+1}\left(x\right)+R_{n+1},}\hfill \end{array}$$

(18)

where

$$\begin{array}{c}|R_{n+1}|\le min\left\{|R_{n+1,1}|,|R_{n+1,2}|\right\},\hfill \\ {\rho }_{n,0}=A_{n,0},\hfill \\ {\rho }_{n,l}=A_{n,l}+B_{n,l-1}.\hfill \end{array}$$

Omitting the truncation error term $R_{n+1}$ from Equation (18), we get the following semidiscrete scheme of Equation (1):

$$\begin{array}{c}{\lambda }_0{v}^{n+1}\left(x\right)-\left({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n}\right)v_{xx}^{n+1}\left(x\right)\hfill \\ {\displaystyle ={\tau }^{\alpha -1}\sum _{l=0}^n{\rho }_{n,l}{v}_{xx}^l\left(x\right)-\sum _{k=0}^n{\lambda }_k{v}^{n+1-k}\left(x\right)+{\tau }^{\alpha }{f}^{n+1}\left(x\right),}\hfill \end{array}$$

(19)

and using the initial and boundary conditions (2), we have

$$\begin{array}{c}{v}^0\left(x\right)=g_0\left(x\right),\hfill \\ v^{n+1}\left(a\right)=0,v^{n+1}\left(b\right)=0.\hfill \end{array}$$

3.2. The Sinc Collocation Method for Spatial Discretization

Now we construct the Sinc collocation method to discrete the semidiscrete scheme (19). The approximation solution $v^n\left(x\right)$ of the semidiscrete scheme (19) can be approximated by

$$v^n\left(x\right)\approx V_m^n\left(x\right)=\sum _{k=-N}^N{c}_j^{n}S(k,h)\circ \varphi \left(x\right),m=2N+1,$$

(20)

where $c_j^n$ is the undetermined coefficient in (20).

$$\begin{array}{c}{\displaystyle {\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}{V}_m^n\left(x\right)=\sum _{j=-N}^N{c}_j^n{\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}[S(j,h)\circ \varphi \left(x\right)]}\hfill \\ {\displaystyle =\sum _{j=-N}^N{c}_j^n[{\varphi }^{″}\left(x\right)S_j^{\left(1\right)}\left(x\right)+{\left({\varphi }^{\prime }\left(x\right)\right)}^2{S}_j^{\left(2\right)}\left(x\right)],}\hfill \end{array}$$

where

$$S_j^{\left(l\right)}={\displaystyle \frac{{\mathrm{d}}^{\left(l\right)}}{\mathrm{d}{\varphi }^{\left(l\right)}}}[S(j,h)\circ \varphi \left(x\right)],l=1,2.$$

It then follows from Lemma 2 that

$${\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}{V}_m^n\left(x_i\right)=\sum _{j=-N}^N{c}_j^n\left[{\varphi }^{″}\left(x_i\right){\displaystyle \frac{{\delta }_{ji}^{\left(1\right)}}{h}}+{\left({\varphi }^{\prime }\left(x_i\right)\right)}^2{\displaystyle \frac{{\delta }_{ji}^{\left(2\right)}}{h^2}}\right].$$

(21)

Substituting (20) and (21) into (19), we have

$$\begin{array}{c}{\displaystyle {\lambda }_0\sum _{j=-N}^{j=N}{c}_j^{n+1}{\delta }_{ji}^{\left(0\right)}-({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n})\sum _{j=-N}^{j=N}{c}_j^{n+1}\left[{\varphi }^{″}\left(x\right){\displaystyle \frac{{\delta }_{ji}^{\left(1\right)}}{h}}+{\left({\varphi }^{\prime }\left(x\right)\right)}^2{\displaystyle \frac{{\delta }_{ji}^{\left(2\right)}}{h^2}}\right]}\hfill \\ {\displaystyle ={\tau }^{\alpha -1}\sum _{l=0}^n\sum _{j=-N}^{j=N}{\rho }_{n,l}{c}_j^{n+1}\left[{\varphi }^{″}\left(x\right){\displaystyle \frac{{\delta }_{ji}^{\left(1\right)}}{h}}+{\left({\varphi }^{\prime }\left(x\right)\right)}^2{\displaystyle \frac{{\delta }_{ji}^{\left(2\right)}}{h^2}}\right]}\hfill \\ {\displaystyle -\sum _{k=0}^n\sum _{j=-N}^{j=N}{\lambda }_k{c}_j^{n+1-k}{\delta }_{ji}^{\left(0\right)}+{\tau }^{\alpha }{f}_i^{n+1}.}\hfill \end{array}$$

(22)

A diagonal matrix of order $2N+1$ is defined as follows

$$D{\left(g\left(x\right)\right)}_{ij}=\left\{\begin{array}{cc}g\left(x_i\right),\hfill & \hfill i=j\\ 0,\hfill & \hfill i\ne j.\end{array}\right.$$

(23)

Multiplying both sides of the above equation by ${\displaystyle \frac{1}{{\left({\varphi }^{\prime }\left(x\right)\right)}^2}}$, we get

$$\begin{array}{c}{\displaystyle {\displaystyle \frac{1}{{\left({\varphi }^{\prime }\left(x_i\right)\right)}^2}}{c}_i^{n+1}-\left({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n}\right)\sum _{j=-N}^N{c}_j^{n+1}\left[{\displaystyle \frac{-{\varphi }^{″}\left(x_i\right)}{{\left({\varphi }^{\prime }\left(x_i\right)\right)}^2}}{\displaystyle \frac{{\delta }_{ij}^{\left(1\right)}}{h}}+{\displaystyle \frac{{\delta }_{ij}^{\left(2\right)}}{h^2}}\right]}\hfill \\ {\displaystyle ={\tau }^{\alpha -1}\sum _{l=0}^n\sum _{j=-N}^N{\rho }_{n,l}{c}_j^l\left[{\displaystyle \frac{-{\varphi }^{″}\left(x_i\right)}{{\left({\varphi }^{\prime }\left(x_i\right)\right)}^2}}{\displaystyle \frac{{\delta }_{ij}^{\left(1\right)}}{h}}+{\displaystyle \frac{{\delta }_{ij}^{\left(2\right)}}{h^2}}\right]+{\displaystyle \frac{1}{{\left({\varphi }^{\prime }\left(x_i\right)\right)}^2}}\sum _{k=0}^n{c}_i^k+{\displaystyle \frac{{\tau }^{\alpha }}{{\left({\varphi }^{\prime }\left(x_i\right)\right)}^2}}{f}_i^{n+1}.}\hfill \end{array}$$

(24)

Writing the above Equation (24) in matrix form as

$$\begin{array}{c}D\left[{\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^2\right]C^{n+1}-\left({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n}\right)\left[{\displaystyle \frac{1}{h}}D{\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^{\prime }{I}^{\left(1\right)}+{\displaystyle \frac{1}{h^2}}{I}^{\left(2\right)}\right]C^{n+1}\hfill \\ ={\tau }^{\alpha }D\left[{\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^2\right]F^{n+1}+D\left[{\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^2\right](C^n+C^{n-1}+\dots +C^1+C^0)\hfill \\ {\displaystyle +{\tau }^{\alpha -1}\sum _{l=0}^n{\rho }_{n,l}\left[{\displaystyle \frac{1}{h}}D{\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^{\prime }{I}^{\left(1\right)}+{\displaystyle \frac{1}{h^2}}{I}^{\left(2\right)}\right]C^l,}\hfill \end{array}$$

(25)

or in a compact form as

$$P{C}^{n+1}=R\left({\tau }^{\alpha }{F}^{n+1}+\sum _{m=0}^n{C}^m\right)+{\tau }^{\alpha -1}\sum _{l=0}^n{\rho }_{n,l}Q{C}^l,$$

(26)

where

$$\begin{array}{c}R=D\left[{\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^2\right],\hfill \\ Q={\displaystyle \frac{1}{h}}D{\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^{\prime }{I}^{\left(1\right)}+{\displaystyle \frac{1}{h^2}}{I}^{\left(2\right)},\hfill \\ C^{n+1}={(c_{-N}^{n+1},c_{-N+1}^{n+1},\dots ,c_N^{n+1})}^{\mathrm{T}},\hfill \\ F^{n+1}={(f_{-N}^{n+1},f_{-N+1}^{n+1},\dots ,f_N^{n+1})}^{\mathrm{T}},\hfill \\ P=R-B_{n,n}Q.\hfill \end{array}$$

(27)

If we set

$$G^{n+1}=R\left({\tau }^{\alpha }{F}^{n+1}+\sum _{m=0}^n{C}^m\right)+{\tau }^{\alpha -1}\sum _{l=0}^n{\rho }_{n,l}Q{C}^l,$$

then Equation (26) can be written as follows:

$$P{C}^{n+1}=G^{n+1},$$

(28)

with the additional initial condition

$$C^0={(V_0\left(x_{-N}\right),V_0\left(x_{-N+1}\right),\cdots ,V_0\left(x_N\right))}^{\mathrm{T}}.$$

For each n, Formula (28) is a system of $2N+1$ order linear equations including $2N+1$ equations. By solving this system of linear equations, the coefficients of the numerical solutions (20) can be obtained.

4. Convergence Analysis

In this section, we aim to analyze the convergence of the semidiscrete Equation (19) for the FPIEs (1)–(3).

For the sake of convenience, the semidiscrete Equation (19) can be rewritten as

$${\lambda }_0{v}^{n+1}\left(x\right)-\left({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n}\right)v_{xx}^{n+1}\left(x\right)=g\left(x\right),$$

(29)

where

$$g\left(x\right)={\tau }^{\alpha -1}\sum _{l=0}^n{\rho }_{n,l}{v}_{xx}^l\left(x\right)-\sum _{k=0}^n{\lambda }_k{v}^{n+1-k}\left(x\right)+{\tau }^{\alpha }{f}^{n+1}\left(x\right).$$

The numerical solution $W_m^{n+1}\left(x\right)$ of Equation (29) at the point $x_j$ can be obtained by

$$W_m^{n+1}\left(x\right)=\sum _{j=-N}^N{v}^{n+1}\left(x_j\right)S_j\left(x\right),$$

(30)

To obtain the bound of $|V^{n+1}\left(x\right)-V_m^{n+1}\left(x\right)|$, we can start by estimating the boundary of $|V_m^{n+1}\left(x\right)-W_m^{n+1}\left(x\right)|$.

Lemma 5

([51]). For $x\in {\varphi }^{-1}$ and the matrix P defined by Equation (27), we have

$${\displaystyle \frac{P+P^{*}}{2}}=H-{\displaystyle \frac{B_{n,n}}{h^2}}{I}^{\left(2\right)},$$

where $P^{*}$ is the conjugate transpose of P and

$$H=D\left[\mathrm{Re}\left({\left({\displaystyle \frac{1}{{\varphi }^{\prime }}}\right)}^2\right)\right]-{\displaystyle \frac{B_{n,n}}{2h}}\left\{D\left[{\left(\frac{1}{{\varphi }^{\prime }}\right)}^{\prime }\right]I^{\left(1\right)}-I^{\left(1\right)}D\left[{\left(\overline{\frac{1}{{\varphi }^{\prime }}}\right)}^{\prime }\right]\right\}.$$

If the eigenvalues of the matrix H are non-negative, then there exists a constant $c_0$ that doesn’t depend on N, such that

$$\parallel P^{-1}{\parallel }_2\le {\displaystyle \frac{4dN}{\beta \pi B_{n,n}}}\left(1+\frac{C_0}{N}\right),$$

for a sufficiently large $N.$

Theorem 1.

Suppose $V_m^{n+1}\left(x\right)$ is an approximate solution of Equation (19), $W_m^{n+1}$ is an approximate solution of Equation (1). Then, there exists a constant $C_4$ that doesn’t depend on N, such that

$$\underset{x\in [a,b]}{sup}|V_m^{n+1}\left(x\right)-W_m^{n+1}\left(x\right)|\le C_4{N}^3{\mathrm{e}}^{-\sqrt{\pi d\beta N}}.$$

Proof. 

By Equations (20) and (30) and the Cauchy–Schwarz inequality, we gain

$$\begin{array}{c}{\displaystyle |V_m^{n+1}\left(x\right)-W_m^{n+1}\left(x\right)|=|\sum _{j=-N}^N{c}_j^{n+1}{S}_j\left(x\right)-\sum _{j=-N}^N{v}^{n+1}\left(x_j\right)S_j\left(x\right)|}\hfill \\ {\displaystyle \le {\left(\sum _{j=-N}^N{|c_j^{n+1}-v^{n+1}\left(x_j\right)|}^2\right)}^{\frac{1}{2}}{\left(\sum _{j=-N}^N{\left|S_j\left(x\right)\right|}^2\right)}^{\frac{1}{2}}.}\hfill \end{array}$$

(31)

Since ${\displaystyle {\left(\sum _{j=-N}^N{\left|S_j\left(x\right)\right|}^2\right)}^{\frac{1}{2}}\le C_1,}$ where $C_1$ is a constant independent of N, we obtain

$$|V_m^{n+1}\left(x\right)-W_m^{n+1}\left(x\right)|\le C_1{\parallel C^{n+1}-U^{n+1}\parallel }_2,$$

(32)

where $C^{n+1}$ is given by (27) and denoting the vector $V^{n+1}$ by

$$U^{n+1}={\left(v^{n+1}\left(x_{-N}\right),v^{n+1}\left(x_{-N+1}\right),\cdots ,v^{n+1}\left(x_N\right)\right)}^{\mathrm{T}}.$$

Based on (19) and (28), we have

$$\parallel C^{n+1}-U^{n+1}{\parallel }_2=\parallel P^{-1}(P{C}^{n+1}-P{U}^{n+1}){\parallel }_2\le \parallel P^{-1}{\parallel }_2{\parallel P{U}^{n+1}-G^{n+1}\parallel }_2.$$

(33)

For simplicity, we denote

$$r_k={\left(P{U}^{n+1}-G^{n+1}\right)}_k,k=-N,\cdots ,N,$$

and using Equation (29), we obtain

$$\begin{array}{c}|r_k|=|g\left(x_k\right)-g_m\left(x_k\right)|\hfill \\ =|v^{n+1}\left(x_k\right)-\left({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n}\right){\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}{v}^{n+1}\left(x_k\right)-V_m^{n+1}\left(x_k\right)+\left({\tau }^{\alpha }+{\tau }^{\alpha -1}{B}_{n,n}\right){\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}{V}_m^{n+1}\left(x_k\right)|\hfill \\ \le |v^{n+1}\left(x_k\right)-V_m^{n+1}\left(x_k\right)|+B_{n,n}|{\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}{v}^{n+1}\left(x_k\right)-{\displaystyle \frac{{\mathrm{d}}^2}{\mathrm{d}{x}^2}}{V}_m^{n+1}\left(x_k\right)|.\hfill \end{array}$$

(34)

Now, using Theorem 1, we have

$$\begin{array}{c}\parallel r_k\parallel \le C_2{N}^{\frac{1}{2}}{\mathrm{e}}^{-\sqrt{\pi d\beta N}}+B_{n,n}{C}_3{N}^{\frac{3}{2}}{\mathrm{e}}^{-\sqrt{\pi d\beta N}}\hfill \\ \le {\mathrm{e}}^{-\sqrt{\pi d\beta N}}\left(C_2{N}^{\frac{3}{2}}+B_{n,n}{C}_3{N}^{\frac{3}{2}}\right)\hfill \\ =K{N}^{\frac{3}{2}}{\mathrm{e}}^{-\sqrt{\pi d\beta N}},\hfill \end{array}$$

(35)

where $C_2$ and $C_3$ are constants independent of N and $K=C_2+B_{n,n}{C}_3$.

$$\parallel P{U}^{n+1}-G^{n+1}{\parallel }_2\le \sqrt{2N+1}{\parallel P{U}^{n+1}-G^{n+1}\parallel }_{\infty },$$

and using inequality (35), we get

$$\parallel P{U}^{n+1}-G^{n+1}{\parallel }_2\le \sqrt{2}K{N}^2{\mathrm{e}}^{-\sqrt{\pi d\beta N}}.$$

(36)

Now, substituting (36) into (33), we have

$$\parallel C^{n+1}-U^{n+1}{\parallel }_2\le {\displaystyle \frac{4\sqrt{2}dK(1+C_0)}{\alpha \pi B_{n,n}}}{N}^3{\mathrm{e}}^{-\sqrt{\pi d\beta N}}.$$

(37)

Based on (32) and (37), we get

$$\underset{x\in [a,b]}{sup}|V_m^{n+1}\left(x\right)-W_m^{n+1}\left(x\right)|\le C_4{N}^3{\mathrm{e}}^{-\sqrt{\pi d\beta N}},$$

(38)

where $C_4={\displaystyle \frac{4\sqrt{2}dK(1+C_0)}{\alpha \pi B_{n,n}}}.$ □

Theorem 2.

Suppose $V^{n+1}\left(x\right)$ be the analytical solution of (29), $V_m^{n+1}\left(x\right)$ be its Sinc approximation defined by (20). Then, there exists a constant $C_7$ that doesn’t depend on N, such that

$$\underset{x\in [a,b]}{sup}|V^{n+1}\left(x\right)-V_m^{n+1}\left(x\right)|\le C_7{N}^3{\mathrm{e}}^{-\sqrt{\pi d\beta N}}.$$

Proof. 

Using the triangular inequality, we get

$$|V^{n+1}\left(x\right)-V_m^{n+1}\left(x\right)|\le |V^{n+1}\left(x\right)-W_m^{n+1}\left(x\right)|+|W_m^{n+1}\left(x\right)-V_m^{n+1}\left(x\right)|.$$

(39)

Using Theorem 1, we can get

$$|V^{n+1}\left(x\right)-W_m^{n+1}\left(x\right)|\le C_5{N}^3{\mathrm{e}}^{-\sqrt{\pi d\beta N}}.$$

(40)

where $C_5$ is a constant independent of N. Based on Theorem 2, there exists a constant $C_6$ that does not depend on N such that

$$|W_m^{n+1}\left(x\right)-V_m^{n+1}\left(x\right)|\le C_6{N}^3{\mathrm{e}}^{-\sqrt{\pi d\beta N}},$$

(41)

Finally, we conclude

$$\underset{x\in [a,b]}{sup}|V^{n+1}\left(x\right)-V_m^{n+1}\left(x\right)|\le C_7{N}^3{\mathrm{e}}^{-\sqrt{\pi d\beta N}},$$

(42)

where $C_7=max\{C_5,C_6\}$. □

5. Numerical Results

In this section, some numerical calculations are performed to demonstrate the validity and accuracy of our method. In numerical examples, we set parameters $d=\frac{\pi }{2}$ and $\beta =1$ and then $h=\frac{\pi }{\sqrt{2N}}$. All numerical computations are carried out using Matlab 7.14 running on a Lenovo PC (Lenovo, Quarry Bay, Hong Kong) with a 1.6 GHz Intel Core i5-4200 CPU (Intel Corporation, Santa Clara, CA, USA) and 4 GB RAM installed.

To illustrate the accuracy of our method, the error analysis is calculated according to the maximum norm errors, defined as:

$$e_{\infty }(h,\tau )=\underset{0\le n\le N}{max}{\parallel V^n-v^n\parallel }_{\infty }.$$

Furthermore, the temporal convergence order can be expressed by

$$rat{e}_1={log}_2\left({\displaystyle \frac{e_{\infty }(h,2\tau )}{e_{\infty }(h,\tau )}}\right).$$

Example 1.

We consider Equations (1)–(3) with the analytical solution

$$v(x,t)=t^{2}x(x-1),$$

where $0<x<1$, $0<t<1$ and

$$f(x,t)={\displaystyle \frac{2}{\Gamma (3-\alpha )}}x(x-1)t^{2-\alpha }-2t^2-{\displaystyle \frac{4\Gamma (1/2)}{\Gamma (7/2)}}{t}^{\frac{5}{2}}.$$

Table 1. The maximum norm errors and temporal convergence orders with $N=32$.

$\mathit{\alpha }$	$\mathit{\tau }$	${\mathit{e}}_{\mathit{\infty }}(\mathit{h},\mathit{\tau })$	${\mathit{Rate}}_1$
0.1	1/10	2.71422 $\times {10}^{-4}$	∗
	1/20	6.89322 $\times {10}^{-5}$	1.97729
	1/40	1.74371 $\times {10}^{-5}$	1.98507
	1/80	4.40355 $\times {10}^{-6}$	1.98542
0.3	1/10	2.71798 $\times {10}^{-4}$	∗
	1/20	6.90718 $\times {10}^{-5}$	1.97637
	1/40	1.74383 $\times {10}^{-5}$	1.98584
	1/80	4.41692 $\times {10}^{-6}$	1.98115
0.5	1/10	2.78117 $\times {10}^{-4}$	∗
	1/20	7.06728 $\times {10}^{-5}$	1.97647
	1/40	1.78836 $\times {10}^{-5}$	1.98252
	1/80	4.51750 $\times {10}^{-6}$	1.98504
0.7	1/10	2.90936 $\times {10}^{-4}$	∗
	1/20	7.37990 $\times {10}^{-5}$	1.97903
	1/40	1.86602 $\times {10}^{-5}$	1.98365
	1/80	4.71342 $\times {10}^{-6}$	1.98512

The asterisk (∗) symbol indicates that the temporal convergence order cannot be calculated.

represents the maximum norm errors and the temporal convergence order for $N=32$ and $\alpha =0.1,0.3,0.5,0.7$ with different values of time step size.

Table 2. Comparison of the maximum norm errors and temporal convergence orders with $N=100$.

$\mathit{\alpha }$	$\mathit{\tau }$	${\mathit{e}}_{\mathit{\infty }}(\mathit{h},\mathit{\tau })$	${\mathit{Rate}}_1$	${\mathit{e}}_{\mathit{\infty }}(\mathit{h},\mathit{\tau })$ [19]	${\mathit{Rate}}_1$ [19]
0.6	1/10	2.7821 $\times {10}^{-4}$	∗	2.4541 $\times {10}^{-4}$	∗
	1/20	7.0565 $\times {10}^{-5}$	1.9791	9.4738 $\times {10}^{-5}$	1.3732
	1/40	1.7991 $\times {10}^{-5}$	1.9717	3.8919 $\times {10}^{-5}$	1.2835
	1/80	4.5933 $\times {10}^{-6}$	1.9696	1.6656 $\times {10}^{-5}$	1.2244
	1/160	1.1502 $\times {10}^{-6}$	1.9976	7.3463 $\times {10}^{-6}$	1.1810
	1/320	2.8781 $\times {10}^{-7}$	1.9987	3.3191 $\times {10}^{-6}$	1.1462

The asterisk (∗) symbol indicates that the temporal convergence order cannot be calculated.

shows a comparative study for the presented method and the method in [19]. It can be observed from the table that the numerical results are better than the method in [19]. The maximum norm errors for $\alpha =0.8$ and $\tau =\frac{1}{1000}$ with different values of N are plotted in Figure 1. At the same time, it is clear from the figure that the presented scheme converges at an exponential rate as N increases. From these diagrams, it can be seen that the results are in excellent agreement with the theoretical analysis.

Example 2.

We consider Equations (1)–(3) with the analytical solution

$$v(x,t)=t sin\left(\pi x\right),$$

where $0<x<pi$, $0<t<1$ and

$$f(x,t)=sin\left(\pi x\right)\left({\displaystyle \frac{1}{\Gamma (2-\alpha )}}{t}^{1-\alpha }+{\pi }^{2}t+{\displaystyle \frac{4}{3}}{\pi }^2{t}^{\frac{3}{2}}\right).$$

Table 3. The maximum norm errors and temporal convergence orders with $N=128$.

$\mathit{\alpha }$	$\mathit{\tau }$	${\mathit{e}}_{\mathit{\infty }}(\mathit{h},\mathit{\tau })$	${\mathit{Rate}}_1$
0.2	1/10	1.4168 $\times {10}^{-6}$	∗
	1/20	3.5892 $\times {10}^{-7}$	1.9809
	1/40	9.0501 $\times {10}^{-8}$	1.9877
	1/80	2.2729 $\times {10}^{-8}$	1.9934
0.4	1/10	6.0618 $\times {10}^{-6}$	∗
	1/20	1.5345 $\times {10}^{-6}$	1.9819
	1/40	3.8691 $\times {10}^{-7}$	1.9877
	1/80	9.7218 $\times {10}^{-8}$	1.9918
0.6	1/10	2.1210 $\times {10}^{-5}$	∗
	1/20	5.3633 $\times {10}^{-6}$	1.9836
	1/40	1.3508 $\times {10}^{-6}$	1.9894
	1/80	3.3889 $\times {10}^{-7}$	1.9949
0.8	1/10	6.8883 $\times {10}^{-5}$	∗
	1/20	1.7426 $\times {10}^{-5}$	1.9829
	1/40	4.3913 $\times {10}^{-6}$	1.9885
	1/80	1.1017 $\times {10}^{-6}$	1.9949

The asterisk (∗) symbol indicates that the temporal convergence order cannot be calculated.

represents the maximum norm errors and the temporal convergence order for $N=128$ and $\alpha =0.2,0.4,0.6,0.8$ with different values of time step size. The maximum norm errors for $\alpha =0.4$ and $\tau =\frac{1}{512}$ with different values of N are plotted in Figure 2. The figure also shows that our presented scheme converges at an exponential rate as N increases. Figure 3 depicts the graph of the numerical solution and the exact solution with $\alpha =0.5$, $\tau =\frac{1}{512}$ and $N=64$. These figures confirm that the proposed method solution is in good agreement with the exact solution.

6. Conclusions

In the present article, we have presented and analyzed an efficient numerical algorithm for solving FPIEs with a weakly singular kernel. In this technique, the WSGD operator is applied for discretization of the time fractional derivative and the Sinc collocation method is used for discretization of the space derivative. Convergence analysis of our scheme is theoretically proven, and it is shown that the numerical solution converges to the exact solution at the exponential rate in space. Numerical experiments were provided to verify the theoretical results. In the future, we intend to extend the method for solving the higher space dimension equation, which is straightforward, in view of the potential applications.