A Higher-Order Numerical Scheme for Two-Dimensional Nonlinear Fractional Volterra Integral Equations with Uniform Accuracy

Abstract

In this paper, based on the modified block-by-block method, we propose a higher-order numerical scheme for two-dimensional nonlinear fractional Volterra integral equations with uniform accuracy. This approach involves discretizing the domain into a large number of subdomains and using biquadratic Lagrangian interpolation on each subdomain. The convergence of the high-order numerical scheme is rigorously established. We prove that the numerical solution converges to the exact solution with the optimal convergence order $O(h_x^{4-\alpha }+h_y^{4-\beta })$ for $0<\alpha ,\beta <1$. Finally, experiments with four numerical examples are shown, to support the theoretical findings and to illustrate the efficiency of our proposed method.

1. Introduction

Volterra integral equations (VIEs) have significant applications in various fields of applied science and engineering, such as superfluidity, elasticity, electromagnetism, electrostatics, potential theory, geophysics, etc. Singular and weakly singular integral equations are of particular interest, since they are used to solve inverse boundary value problems whose domains are fractal curves, where classical calculus cannot be used. Abel equations and other fractional-order integral equations have been studied extensively and are used in modeling various phenomena in biophysics, viscoelasticity and electrical circuits [1,2]. In this paper, we consider using the biquadrature formula to solve the following two-dimensional nonlinear VIEs with a fractional-order weakly singular kernel:

(1)$u(x,y)=f(x,y)+{\int }_a^x{\int }_c^y\frac{\kappa (x,y,s,t,u(s,t))}{{(x-s)}^{\alpha }{(y-t)}^{\beta }}dtds,(x,y)\in D,0<\alpha ,\beta <1,$

where $f(x,y)$ and $\kappa (x,y,s,t,u(s,t))$ are given continuous functions defined on $D=[a,b]\times [c,d]$, $\Omega =D\times D\times R$, and $u(x,y)$ is an unknown function defined on D. In addition, we assume that the VIEs (1) have a smooth solution $u(x,y)$ and that $\kappa (x,y,s,t,u(s,t)$ satisfies the following condition

(2)$\begin{array}{c}\hfill |\kappa (x,y,s,t,u_1(s,t))-\kappa (x,y,s,t,u_2(s,t))|\le L|u_1(s,t)-u_2(s,t)|,L>0.\end{array}$

The classical block-by-block method is an efficient numerical algorithm for VIEs. A general block-by-block method for one-dimensional linear VIEs with a nonsingular kernel function was given in [3]. In [4], systems of one-dimensional nonlinear VIEs with a nonsingular kernel function were solved by the block-by-block method. A new high-order numerical scheme for fractional ordinary differential equations was given in [5], called the modified block-by-block method. In [6], radial basis functions were used to solve for the second kind of nonlinear VIEs in which the kernel function satisfies the Lipschitz condition. The authors of [7] used the Runge–Kutta method and the block-by-block method to solve the second kind of nonlinear VIEs with a continuous kernel. In [8], a new transformation was used to solve the VIEs by using a Laplace transform. The Bernstein approximation method was used for VIEs of the third kind in [9]. In [10], the spectral collocation method with graded meshes was used to solve the second kind of VIEs with a weakly singular kernel. Two-dimensional VIEs were studied by an Euler-type numerical scheme in [11]. A method based on two-dimensional Euler polynomials combined with the Gauss–Jacobi quadrature formula was used to solve two-dimensional VIEs with fractional-order weakly singular kernels in [12]. General linear methods were implemented with a variable step size for VIEs with a sufficiently smooth kernel function, and the corresponding MATLAB code was developed in [13]. VIEs with a nonlinear weakly singular kernel function were solved by an hp-version collocation method in conjunction with Jacobi polynomials in [14]. The optimal homotopy asymptotic method was used for linear and nonlinear two-dimensional VIEs in [15]. In [16], an iterated multi-Galerkin method was used to solve VIEs with a weakly singular kernel, and the proof of convergence rates was given by the projected methods. Multivariate Bernstein polynomials were used to solve multidimensional linear and nonlinear VIEs with fractional weakly singular kernel functions in [17]. In [18], a Jacobi spectral collocation method was proposed for a class of nonlinear VIEs with kernels of the form $x^{\beta }{(z-x)}^{-\alpha }g(y(x))$, where $\alpha \in (0,1),\beta >0$ and $g(y)$ is a nonlinear function. Numerical solutions for weakly singular VIEs were presented based on Chebyshev and Legendre pseudo-spectral Galerkin methods in [19]. Some more recent developments are shown in [20,21,22]. The study of numerical solutions for high-dimensional VIEs with singular nonlinear kernel functions remains an important research issue. To date, there have been few reports on high-order numerical algorithms for two-dimensional VIEs with singular nonlinear kernel functions.

In this paper, we construct a new high-order numerical solution for two-dimensional nonlinear fractional VIEs by using the techniques of [5] with uniform accuracy. To avoid degeneracy near the two boundary layers, the proposed scheme couples the solutions at the two boundary layers. Thus, no smaller meshes are needed to achieve the sharp numerical order. Such coupling is not required in the other subdomains. The convergence analysis is based on a novel technology that couples the ideas of [5] and the Gronwall inequality. Hence, the optimal convergence order $O(h_x^{4-\alpha }+h_y^{4-\beta })$ for $0<\alpha ,\beta <1$ can be achieved for a sufficiently smooth solution and a general nonlinear kernel function.

This paper is arranged as follows. In Section 2, the higher-order numerical scheme is proposed. The estimation of the truncation errors of the constructed higher-order numerical scheme is given in Section 3. A convergence analysis of the higher-order numerical scheme is given in Section 4. In Section 5, four numerical experiments are presented to illustrate the efficiency of the the high-order numerical approach and to support our theoretical findings. Finally, some conclusions arising from this work are drawn in Section 6.

2. Higher-Order Numerical Scheme of Two-Dimensional Nonlinear Fractional VIEs

In this section, we consider the approximate evaluation of two-dimensional nonlinear fractional VIEs. In order to construct the numerical scheme of Equation (1), the region D is divided into $2M\times 2N$ equal subdomains with size $h_x=\frac{b-a}{2M},h_y=\frac{d-c}{2N}$, where M and N are positive integers. We denote $x_i=a+i{h}_x$,$y_j=c+j{h}_y,i=0,1,\cdots ,2M;$$j=0,1,\cdots ,2N$, and the numerical solution of formula (1) at point $(x_i,y_j)$ is denoted by $u_{i,j}$, where $u_{i,0}=f(x_i,0),u_{0,j}=f(0,y_j)$. For convenience of narration, we let ${\kappa }_{i,j}(s,t,u(s,t))=\kappa (x_i,y_j,s,t,u(s,t)),f_{i,j}=f(x_i,y_j)$.

Firstly, we propose a high-order scheme for the nonlinear fractional VIEs. Let ${\phi }_{k,i}(s),$$k=0,1,2;i\in \mathbb{N}$ and ${\varphi }_{k,j}(t),k=0,1,2;j\in \mathbb{N}$ be the basis functions of quadratic interpolation polynomials at points $x_i,x_{i+1},x_{i+2}$ and $y_j,y_{j+1},y_{j+2}$, respectively, where ${\phi }_{k,i}(s)$ and ${\varphi }_{k,j}(t),k=0,1,2$, are defined as follows

$\begin{array}{cc}\hfill & {\phi }_{0,i}(s)=\frac{(s-x_{i+1})(s-x_{i+2})}{2h_x^2},{\phi }_{1,i}(s)=\frac{(s-x_i)(s-x_{i+2})}{-h_x^2},{\phi }_{2,i}(s)=\frac{(s-x_i)(s-x_{i+1})}{2h_x^2};\hfill \\ \hfill & {\varphi }_{0,j}(t)=\frac{(t-y_{j+1})(t-y_{j+2})}{2h_y^2},{\varphi }_{1,j}(t)=\frac{(t-y_j)(t-y_{j+2})}{-h_y^2},{\varphi }_{2,j}(t)=\frac{(t-y_j)(t-y_{j+1})}{2h_y^2}.\hfill \end{array}$

We estimate that $u(x,y)$ at point $(x_1,y_1)$ has the form

(3)$\begin{array}{ccc}\hfill u(x_1,y_1)& =& f_{1,1}+{\int }_a^{x_1}{\int }_c^{y_1}\frac{{\kappa }_{1,1}(s,t,u(s,t))}{{(x_1-s)}^{\alpha }{(y_1-t)}^{\beta }}dtds\hfill \\ \hfill & \approx & f_{1,1}+{\int }_a^{x_1}{\int }_c^{y_1}\frac{1}{{(x_1-s)}^{\alpha }{(y_1-t)}^{\beta }}{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\phi }_{i,0}(s){\varphi }_{j,0}(t){\kappa }_{1,1}(x_i,y_j,u_{i,j})dtds\hfill \\ & =& f_{1,1}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\omega }_1^{i,0}{\widehat{\omega }}_1^{j,0}{\kappa }_{1,1}(x_i,y_j,u_{i,j}),\hfill \end{array}$

with

(4)$\begin{array}{c}\hfill {\omega }_1^{i,0}={\int }_a^{x_1}\frac{{\phi }_{i,0}(s)}{{(x_1-s)}^{\alpha }}ds,i=0,1,2;{\widehat{\omega }}_1^{j,0}={\int }_c^{y_1}\frac{{\varphi }_{j,0}(t)}{{(y_1-t)}^{\beta }}dt,j=0,1,2.\end{array}$

To compute $u(x_2,y_1),u(x_1,y_2)$ and $u(x_2,y_2)$, we use the following approximation:

(5)$\begin{array}{cc}\hfill & u(x_2,y_1)=f_{2,1}+{\int }_a^{x_2}{\int }_c^{y_1}\frac{{\kappa }_{2,1}(s,t,u(s,t))}{{(x_2-s)}^{\alpha }{(y_1-t)}^{\beta }}dtds\approx f_{2,1}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{A}_2^{i,0}{\widehat{\omega }}_1^{j,0}{\kappa }_{2,1}(x_i,y_j,u_{i,j}),\hfill \end{array}$

(6)$\begin{array}{cc}\hfill & u(x_1,y_2)=f_{1,2}+{\int }_a^{x_1}{\int }_c^{y_2}\frac{{\kappa }_{1,2}(s,t,u(s,t))}{{(x_1-s)}^{\alpha }{(y_2-t)}^{\beta }}dtds\approx f_{1,2}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\omega }_1^{i,0}{\widehat{A}}_2^{j,0}{\kappa }_{1,2}(x_i,y_j,u_{i,j}),\hfill \end{array}$

(7)$\begin{array}{cc}\hfill & u(x_2,y_2)=f_{2,2}+{\int }_a^{x_2}{\int }_c^{y_2}\frac{{\kappa }_{2,2}(s,t,u(s,t))}{{(x_2-s)}^{\alpha }{(y_2-t)}^{\beta }}dtds\approx f_{2,2}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{A}_2^{i,0}{\widehat{A}}_2^{j,0}{\kappa }_{2,2}(x_i,y_j,u_{i,j}),\hfill \end{array}$

where

(8)$\begin{array}{c}\hfill A_2^{i,0}={\int }_a^{x_2}\frac{{\phi }_{i,0}(s)}{{(x_2-s)}^{\alpha }}ds,i=0,1,2;{\widehat{A}}_2^{j,0}={\int }_c^{y_2}\frac{{\varphi }_{j,0}(t)}{{(y_2-t)}^{\beta }}dt,j=0,1,2.\end{array}$

Note that computing $u_{1,1}$ through (3) requires the values of $\kappa$ (or indirectly, the values of u) at $x_1,x_2$ and $y_1,y_2$. Particularly, the dependence of $u_{1,1}$ on ${\kappa }_{2,1},{\kappa }_{1,2}$ and ${\kappa }_{2,2}$ means that (3), (5)–(7) must be solved simultaneously with the scheme.

We further estimate $u(x_{2m+l},y_r),m=1,\cdots ,M-1$ and $u(x_l,y_{2n+r}),l,r=1,2,$$n=1,\cdots ,N-1$, assuming $u_{i,r},i=0,1,\cdots ,2m$ and $u_{l,j},j=0,1,\cdots ,2n$ are already known. For $u(x_{2m+1},y_1)$, we have

(9)$\begin{array}{cc}\hfill u(x_{2m+1},y_1)=& f_{2m+1,1}+{\int }_a^{x_1}{\int }_c^{y_1}\frac{{\kappa }_{2m+1,1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_1-t)}^{\beta }}dtds\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{\kappa }_{2m+1,1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_1-t)}^{\beta }}dtds\hfill \\ \hfill \approx & f_{2m+1,1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\int }_a^{x_1}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_1-t)}^{\beta }}\hfill \\ \hfill & \times {\phi }_{k,0}(s){\varphi }_{q,0}(t){\kappa }_{2m+1,1}(x_k,y_q,u_{k,q})dtds\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_1-t)}^{\beta }}\hfill \\ \hfill & \times {\phi }_{k,2i-1}(s){\varphi }_{q,0}(t){\kappa }_{2m+1,1}(x_{2i-1+k},y_q,u_{2i-1+k,q})dtds\hfill \\ \hfill =& f_{2m+1,1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{\omega }}_1^{q,0}{\kappa }_{2m+1,1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{\omega }}_1^{q,0}{\kappa }_{2m+1,1}(x_{2i-1+k},y_q,u_{2i-1+k,q}),\hfill \end{array}$

where ${\widehat{\omega }}_1^{q,0}$ is in (4) and

(10)$\begin{array}{cc}\hfill & {\omega }_{2m+1}^{k,0}={\int }_a^{x_1}{(x_{2m+1}-s)}^{-\alpha }{\phi }_{k,0}(s)ds,k=0,1,2,\hfill \end{array}$

(11)$\begin{array}{cc}\hfill & A_{2m+1}^{k,i}={\int }_{x_{2i-1}}^{x_{2i+1}}{(x_{2m+1}-s)}^{-\alpha }{\phi }_{k,2i-1}(s)ds,k=0,1,2;i=1,\cdots ,m.\hfill \end{array}$

We estimate $u(x_{2m+2},y_1)$ and $u(x_{2m+l},y_2),l=1,2$ using the following approximations

$\begin{array}{cc}\hfill u(x_{2m+2},y_1)=& f_{2m+2,1}+{\displaystyle \sum _{i=0}^m}{\int }_{x_{2i}}^{x_{2i+2}}{\int }_c^{y_1}\frac{{\kappa }_{2m+2,1}(s,t,u(s,t))}{{(x_{2m+2}-s)}^{\alpha }{(y_1-t)}^{\beta }}dtds\hfill \end{array}$

(12)$\begin{array}{cc}\hfill \approx & f_{2m+2,1}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{\omega }}_1^{q,0}{\kappa }_{2m+2,1}(x_{2i+k},y_q,u_{2i+k,q}),\hfill \\ \hfill u(x_{2m+1},y_2)=& f_{2m+1,2}+{\int }_a^{x_1}{\int }_c^{y_2}\frac{{\kappa }_{2m+1,2}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_2-t)}^{\beta }}dtds\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_2}\frac{{\kappa }_{2m+1,2}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_2-t)}^{\beta }}dtds\hfill \\ \hfill \approx & f_{2m+1,2}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_2^{q,0}{\kappa }_{2m+1,2}(x_k,y_q,u_{k,q})\hfill \end{array}$

(13)$\begin{array}{cc}\hfill & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_2^{q,0}{\kappa }_{2m+1,2}(x_{2i-1+k},y_q,u_{2i-1+k,q}),\hfill \\ \hfill u(x_{2m+2},y_2)=& f_{2m+2,2}+{\displaystyle \sum _{i=0}^m}{\int }_{x_{2i}}^{x_{2i+2}}{\int }_c^{y_2}\frac{{\kappa }_{2m+2,2}(s,t,u(s,t))}{{(x_{2m+2}-s)}^{\alpha }{(y_2-t)}^{\beta }}dtds\hfill \end{array}$

(14)$\begin{array}{cc}\hfill \approx & f_{2m+2,2}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{A}}_2^{q,0}{\kappa }_{2m+2,2}(x_{2i+k},y_q,u_{2i+k,q}),\hfill \end{array}$

where ${\widehat{\omega }}_1^{q,0}$ is in (4), ${\widehat{A}}_2^{q,0}$ is in (8), ${\omega }_{2m+1}^{k,0}$ and $A_{2m+1}^{k,i}$ are in (10) and (11), respectively, and

(15)$\begin{array}{c}\hfill A_{2m+2}^{k,i}={\int }_{x_{2i}}^{x_{2i+2}}{(x_{2m+2}-s)}^{-\alpha }{\phi }_{k,2i}(s)ds,k=0,1,2;i=0,1,\cdots ,m.\end{array}$

Next we use the same method to estimate $u(x_l,y_{2n+r}),l,r=1,2$ and directly obtain

(16)$\begin{array}{cc}\hfill u(x_1,y_{2n+1})=& f_{1,2n+1}+{\int }_a^{x_1}{\int }_c^{y_1}\frac{{\kappa }_{1,2n+1}(s,t,u(s,t))}{{(x_1-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^n}{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\kappa }_{1,2n+1}(s,t,u(s,t))}{{(x_1-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill \approx & f_{1,2n+1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_1^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{1,2n+1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_1^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{1,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q}),\hfill \end{array}$

(17)$\begin{array}{cc}\hfill u(x_2,y_{2n+1})=& f_{2,2n+1}+{\int }_a^{x_2}{\int }_c^{y_1}\frac{{\kappa }_{2,2n+1}(s,t,u(s,t))}{{(x_2-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^n}{\int }_a^{x_2}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\kappa }_{2,2n+1}(s,t,u(s,t))}{{(x_2-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill \approx & f_{2,2n+1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_2^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2,2n+1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_2^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q}),\hfill \end{array}$

(18)$\begin{array}{cc}\hfill u(x_1,y_{2n+2})=& f_{1,2n+2}+{\displaystyle \sum _{j=0}^n}{\int }_a^{x_1}{\int }_{y_{2j}}^{y_{2j+2}}\frac{{\kappa }_{1,2n+2}(s,t,u(s,t))}{{(x_1-s)}^{\alpha }{(y_{2n+2}-t)}^{\beta }}dtds\hfill \\ \hfill \approx & f_{1,2n+2}+{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_1^{k,0}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{1,2n+2}(x_k,y_{2j+q},u_{k,2j+q}),\hfill \end{array}$

(19)$\begin{array}{cc}\hfill u(x_2,y_{2n+2})=& f_{2,2n+2}+{\displaystyle \sum _{j=0}^n}{\int }_a^{x_2}{\int }_{y_{2j}}^{y_{2j+2}}\frac{{\kappa }_{2,2n+2}(s,t,u(s,t))}{{(x_2-s)}^{\alpha }{(y_{2n+2}-t)}^{\beta }}dtds\hfill \\ \hfill \approx & f_{2,2n+2}+{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_2^{k,0}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2,2n+2}(x_k,y_{2j+q},u_{k,2j+q}),\hfill \end{array}$

where ${\omega }_1^{k,0}$ and $A_2^{k,0}$ are in (4) and (8), respectively, and

(20)$\begin{array}{cc}\hfill & {\widehat{\omega }}_{2n+1}^{q,0}={\int }_c^{y_1}{(y_{2n+1}-t)}^{-\beta }{\varphi }_{q,0}(t)dt,q=0,1,2,\hfill \end{array}$

(21)$\begin{array}{cc}\hfill & {\widehat{A}}_{2n+1}^{q,j}={\int }_{y_{2j-1}}^{y_{2j+1}}{(y_{2n+1}-t)}^{-\beta }{\varphi }_{q,2j-1}(t)dt,q=0,1,2;j=1,2,\cdots ,n,\hfill \end{array}$

(22)$\begin{array}{cc}\hfill & {\widehat{A}}_{2n+2}^{q,j}={\int }_{y_{2j}}^{y_{2j+2}}{(y_{2n+2}-t)}^{-\beta }{\varphi }_{q,2j}(t)dt,q=0,1,2;j=0,1,\cdots ,n.\hfill \end{array}$

Next we approximate $u(x_i,y_j)$ for $i>2,j>2$. We first assume that $u_{i,j},u_{i,2n+1},u_{i,2n+2},$ $u_{2m+1,j},u_{2m+2,j},i=0,1,\cdots ,2m;j=0,1,\cdots ,2n;m=1,\cdots ,$$M-1;n=1,\cdots ,N-1$ are known, and then we construct the approximation of $u(x_{2m+p},y_{2n+q}),p,q=1,2$.

Firstly, for $u(x_{2m+1},y_{2n+1})$, we have

(23)$\begin{array}{ccc}\hfill u(x_{2m+1},y_{2n+1})& =& f_{2m+1,2n+1}+{\int }_a^{x_1}{\int }_c^{y_1}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & & +{\displaystyle \sum _{j=1}^n}{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ & =& f_{2m+1,2n+1}+I_1+I_2+I_3+I_4,\hfill \end{array}$

where $I_1$ is the integral of the first block $[a,x_1]\times [c,y_1]$ and can be calculated by using biquadratic interpolation as follows

(24)$\begin{array}{cc}\hfill & I_1={\int }_a^{x_1}{\int }_c^{y_1}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & \approx {\int }_a^{x_1}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\phi }_{k,0}(s){\varphi }_{q,0}(t){\kappa }_{2m+1,2n+1}(x_k,y_q,u_{k,q})dtds\hfill \\ \hfill & ={\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_k,y_q,u_{k,q}),\hfill \end{array}$

where ${\omega }_{2m+1}^{k,0}$ is shown in (10) and ${\widehat{\omega }}_{2n+1}^{q,0}$ is defined in (20).

For $I_2$, we can obtain

(25)$\begin{array}{cc}\hfill & I_2={\displaystyle \sum _{j=1}^n}{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & \approx {\displaystyle \sum _{j=1}^n}{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\phi }_{k,0}(s){\varphi }_{q,2j-1}(t){\kappa }_{2m+1,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q})dtds\hfill \\ \hfill & ={\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q}),\hfill \end{array}$

where ${\omega }_{2m+1}^{k,0}$ is shown in (10) and ${\widehat{A}}_{2n+1}^{q,j}$ is defined in (20).

For $I_3$, we have

(26)$\begin{array}{cc}\hfill & I_3={\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & \approx {\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}\hfill \\ \hfill & \times {\phi }_{k,2i-1}(s){\varphi }_{q,0}(t){\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_q,u_{2i-1+k,q})dtds\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_q,u_{2i-1+k,q}),\hfill \end{array}$

where $A_{2m+1}^{k,i}$ and ${\widehat{\omega }}_{2n+1}^{q,0}$ are shown in (11) and (20), respectively.

For $I_4$, along the same lines, we have

(27)$\begin{array}{cc}\hfill & I_4={\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\kappa }_{2m+1,2n+1}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & \approx {\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}\hfill \\ \hfill & \times {\phi }_{k,2i-1}(s){\varphi }_{q,2j-1}(t){\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_{2j-1+q},u_{2i-1+k,2j-1+q})dtds\hfill \\ \hfill & ={\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_{2j-1+q},u_{2i-1+k,2j-1+q}),\hfill \end{array}$

where $A_{2m+1}^{k,i}$ is shown in (11) and ${\widehat{A}}_{2n+1}^{q,j}$ is defined in (21).

Bringing (24)–(27) into (23), we obtain

(28)$\begin{array}{ccc}\hfill & & u_{2m+1,2n+1}=f_{2m+1,2n+1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_q,u_{2i-1+k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_{2j-1+q},u_{2i-1+k,2j-1+q}).\hfill \end{array}$

Next, we construct the numerical scheme for $u(x_{2m+2},y_{2n+1})$. We first split the integration domain, and then use the piecewise biquadratic interpolation to approximate the calculation in the corresponding subdomains. Therefore, the numerical scheme of $u(x_{2m+2},y_{2n+1})$ can be given by the following

(29)$\begin{array}{ccc}\hfill & & u(x_{2m+2},y_{2n+1})=f_{2m+2,2n+1}+{\displaystyle \sum _{i=0}^m}{\int }_{x_{2i}}^{x_{2i+2}}{\int }_c^{y_1}\frac{{\kappa }_{2m+2,2n+1}(s,t,u(s,t))}{{(x_{2m+2}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ \hfill & & +{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{j=1}^n}{\int }_{x_{2i}}^{x_{2i+2}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\kappa }_{2m+2,2n+1}(s,t,u(s,t))}{{(x_{2m+2}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds\hfill \\ & \approx & f_{2m+2,2n+1}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+2,2n+1}(x_{2i+k},y_q,u_{2i+k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+2,2n+1}(x_{2i+k},y_{2j-1+q},u_{2i+k,2j-1+q}),\hfill \end{array}$

where $A_{2m+2}^{k,i}$, ${\widehat{\omega }}_{2n+1}^{q,0}$ and ${\widehat{A}}_{2n+1}^{q,j}$ are in (15), (20) and (21), respectively.

As a consequence, $u(x_{2m+1},y_{2n+2})$ and $u(x_{2m+2},y_{2n+2})$ can be approximated as follows

(30)$\begin{array}{ccc}\hfill & & u(x_{2m+1},y_{2n+2})=f_{2m+1,2n+2}+{\displaystyle \sum _{j=0}^n}{\int }_a^{x_1}{\int }_{y_{2j}}^{y_{2j+2}}\frac{{\kappa }_{2m+1,2n+2}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+2}-t)}^{\beta }}dtds\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=0}^n}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j}}^{y_{2j+2}}\frac{{\kappa }_{2m+1,2n+2}(s,t,u(s,t))}{{(x_{2m+1}-s)}^{\alpha }{(y_{2n+2}-t)}^{\beta }}dtds\hfill \\ & & \approx f_{2m+1,2n+2}+{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2m+1,2n+2}(x_k,y_{2j+q},u_{k,2j+q})\hfill \end{array}$

(31)$\begin{array}{ccc}\hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2m+1,2n+2}(x_{2i-1+k},y_{2j+q},u_{2i-1+k,2j+q}),\hfill \\ \hfill & & u(x_{2m+2},y_{2n+2})=f_{2m+2,2n+2}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{j=0}^n}{\int }_{x_{2i}}^{x_{2i+2}}{\int }_{y_{2j}}^{y_{2j+2}}\frac{{\kappa }_{2m+2,2n+2}(s,t,u(s,t))}{{(x_{2m+2}-s)}^{\alpha }{(y_{2n+2}-t)}^{\beta }}dtds\hfill \\ \hfill & & \approx f_{2m+2,2n+2}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2m+2,2n+2}(x_{2i+k},y_{2j+q},u_{2i+k,2j+q}),\hfill \end{array}$

where ${\omega }_{2m+1}^{k,0}$, $A_{2m+1}^{k,i}$, $A_{2m+2}^{k,i}$ and ${\widehat{A}}_{2n+2}^{q,j}$ are shown in (10), (11), (15) and (22), respectively.

To summarize, by combining (3), (5)–(7), (9), (12)–(14), (16)–(19) and (28)–(31) we compose the high-order numerical scheme of Equation (1) as follows

(32)$\begin{array}{ccc}\hfill u_{1,1}& =& f_{1,1}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\omega }_1^{i,0}{\widehat{\omega }}_1^{j,0}{\kappa }_{1,1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill u_{2,1}& =& f_{2,1}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{A}_2^{i,0}{\widehat{\omega }}_1^{j,0}{\kappa }_{2,1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill u_{1,2}& =& f_{1,2}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\omega }_1^{i,0}{\widehat{A}}_2^{j,0}{\kappa }_{1,2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill u_{2,2}& =& f_{2,2}+{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{A}_2^{i,0}{\widehat{A}}_2^{j,0}{\kappa }_{2,2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill u_{2m+1,1}& =& f_{2m+1,1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{\omega }}_1^{q,0}{\kappa }_{2m+1,1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{\omega }}_1^{q,0}{\kappa }_{2m+1,1}(x_{2i-1+k},y_q,u_{2i-1+k,q}),\hfill \\ \hfill u_{2m+2,1}& =& f_{2m+2,1}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{\omega }}_1^{q,0}{\kappa }_{2m+2,1}(x_{2i+k},y_q,u_{2i+k,q}),\hfill \\ \hfill u_{2m+1,2}& =& f_{2m+1,2}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_2^{q,0}{\kappa }_{2m+1,2}(x_k,y_q,u_{k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_2^{q,0}{\kappa }_{2m+1,2}(x_{2i-1+k},y_q,u_{2i-1+k,q}),\hfill \\ \hfill u_{2m+2,2}& =& f_{2m+2,2}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{A}}_2^{q,0}{\kappa }_{2m+2,2}(x_{2i+k},y_q,u_{2i+k,q}),\hfill \\ \hfill u_{1,2n+1}& =& f_{1,2n+1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_1^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{1,2n+1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_1^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{1,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q}),\hfill \\ \hfill u_{2,2n+1}& =& f_{2,2n+1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_2^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2,2n+1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_2^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q}),\hfill \\ \hfill u_{1,2n+2}& =& f_{1,2n+2}+{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_1^{k,0}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{1,2n+2}(x_k,y_{2j+q},u_{k,2j+q}),\hfill \\ \hfill u_{2,2n+2}& =& f_{2,2n+2}+{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_2^{k,0}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2,2n+2}(x_k,y_{2j+q},u_{k,2j+q}),\hfill \\ \hfill u_{2m+1,2n+1}& =& f_{2m+1,2n+1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_k,y_q,u_{k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_k,y_{2j-1+q},u_{k,2j-1+q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_q,u_{2i-1+k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_{2j-1+q},u_{2i-1+k,2j-1+q}),\hfill \\ \hfill u_{2m+2,2n+1}& =& f_{2m+2,2n+1}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+2,2n+1}(x_{2i+k},y_q,u_{2i+k,q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+2,2n+1}(x_{2i+k},y_{2j-1+q},u_{2i+k,2j-1+q}),\hfill \\ \hfill u_{2m+1,2n+2}& =& f_{2m+1,2n+2}+{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2m+1,2n+2}(x_k,y_{2j+q},u_{k,2j+q})\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2m+1,2n+2}(x_{2i-1+k},y_{2j+q},u_{2i-1+k,2j+q}),\hfill \\ \hfill u_{2m+2,2n+2}& =& f_{2m+2,2n+2}+{\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{j=0}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+2}^{k,i}{\widehat{A}}_{2n+2}^{q,j}{\kappa }_{2m+2,2n+2}(x_{2i+k},y_{2j+q},u_{2i+k,2j+q}).\hfill \end{array}$

3. Estimation of the Truncation Errors

Now, we analyze the estimation of the truncation error of (32). We define the truncation error at the point $(x_i,y_j)$ by

(33)$\begin{array}{c}\hfill r_{i,j}:=u(x_i,y_j)-{\overline{u}}_{i,j},\end{array}$

where ${\overline{u}}_{i,j}$ is an approximate value of $u(x_i,y_j)$, which is substituted into the exact solution of (32). For example, we define ${\overline{u}}_{i,j}$ at the point $(x_{2m+1},y_{2n+1})$ as follows

(34)$\begin{array}{ccc}\hfill & & {\overline{u}}_{2m+1,2n+1}=f_{2m+1,2n+1}+{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_k,y_q,u(x_k,y_q))\hfill \\ & & +{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\omega }_{2m+1}^{k,0}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_k,y_{2j-1+q},u(x_k,y_{2j-1+q}))\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{\omega }}_{2n+1}^{q,0}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_q,u(x_{2i-1+k},y_q))\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{A}_{2m+1}^{k,i}{\widehat{A}}_{2n+1}^{q,j}{\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_{2j-1+q},u(x_{2i-1+k},y_{2j-1+q})).\hfill \end{array}$

For ease of notation, let $\frac{{\partial }^3\kappa }{\partial s^3}\doteq {\partial }_s^3\kappa$. Then, the estimation of $r_{i,j}$ is as follows.

Lemma 1.

Let $r_{i,j}$ be the truncation error defined in (33). If $\kappa (·,·,·,·,u(·,·))\in C^4([a,b]\times [c,d])$, then it holds that

$\begin{array}{c}\hfill |r_{i,j}|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),\end{array}$

where C depends only on $\alpha ,\beta ,G_1,G_2$, and $G_1,G_2$ are defined as follows:

(35) $\begin{array}{c}\hfill G_1=\underset{\begin{array}{c}x,s\in [a,b]\\ y,t\in [c,d]\end{array}}{max}(|{\partial }_s^3\kappa (x,y,s,t,u(s,t))|,|{\partial }_t^3\kappa (x,y,s,t,u(s,t))|),\end{array}$

(36) $\begin{array}{c}\hfill G_2=\underset{\begin{array}{c}x,s\in [a,b]\\ y,t\in [c,d]\end{array}}{max}(|{\partial }_s^4\kappa (x,y,s,t,u(s,t))|,|{\partial }_t^4\kappa (x,y,s,t,u(s,t))|).\end{array}$

Proof. 

Firstly, we estimate the truncation error $r_{2m+1,2n+1}$. According to (33), we define the truncation error at the point $(x_{2m+1},y_{2n+1})$. Combining (28), (33) and (34), we have

(37)$\begin{array}{ccc}\hfill & & r_{2m+1,2n+1}=u\left(x_{2m+1},y_{2n+1}\right)-{\overline{u}}_{2m+1,2n+1}\hfill \\ \hfill & =& {\int }_a^{x_1}{\int }_c^{y_1}{(x_{2\mathrm{m}+1}-s)}^{-\alpha }{(y_{2n+1}-t)}^{-\beta }[{\kappa }_{2m+1,2n+1}(s,t,u(s,t))\hfill \\ \hfill & & -{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\phi }_{k,0}(s){\varphi }_{q,0}(t){\kappa }_{2m+1,2n+1}(x_k,y_q,u(x_k,y_q))]dtds\hfill \\ \hfill & & +{\displaystyle \sum _{j=1}^n}{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}{(x_{2\mathrm{m}+1}-s)}^{-\alpha }{(y_{2n+1}-t)}^{-\beta }[{\kappa }_{2m+1,2n+1}(s,t,u(s,t))\hfill \\ \hfill & & -{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\phi }_{k,0}(s){\varphi }_{q,2j-1}(t){\kappa }_{2m+1,2n+1}(x_k,y_{2j-1+q},u(x_k,y_{2j-1+q}))]dtds\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}{(x_{2\mathrm{m}+1}-s)}^{-\alpha }{(y_{2n+1}-t)}^{-\beta }[{\kappa }_{2m+1,2n+1}(s,t,u(s,t))\hfill \\ \hfill & & -{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\phi }_{k,2i-1}(s){\varphi }_{q,0}(t){\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_q,u(x_{2i-1+k},y_q))]dtds\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}{(x_{2m+1}-s)}^{-\alpha }{(y_{2n+1}-t)}^{-\beta }[{\kappa }_{2m+1,2n+1}(s,t,u(s,t))\hfill \\ \hfill & & -{\displaystyle \sum _{k=0}^2}{\displaystyle \sum _{q=0}^2}{\phi }_{k,2i-1}(s){\varphi }_{q,2j-1}(t){\kappa }_{2m+1,2n+1}(x_{2i-1+k},y_{2j-1+q},u(x_{2i-1+k},y_{2j-1+q}))]dtds\hfill \\ \hfill & =& {\int }_a^{x_1}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{R}_{1}dtds+{\displaystyle \sum _{j=1}^n}{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{R}_{2}dtds\hfill \\ \hfill & & +{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{R}_{3}dtds+{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{R}_{4}dtds\hfill \\ & \doteq & r_{2m+1,2n+1}^{(1)}+r_{2m+1,2n+1}^{(2)}+r_{2m+1,2n+1}^{(3)}+r_{2m+1,2n+1}^{(4)}.\hfill \end{array}$

Using Taylor’s theorem, one can obtain that for all $(s,t)\in [a,x_1]\times [c,y_1]$,

$\begin{array}{ccc}\hfill R_1& =& \frac{1}{3!}{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_1(s),t,u({\xi }_1(s),t)){\displaystyle \prod _{k=0}^2}(s-x_k)\hfill \\ & & +{\displaystyle \sum _{k=0}^2}\frac{{\phi }_{k,0}(s)}{3!}{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_k,{\eta }_1(t),u(x_k,{\eta }_1(t))){\displaystyle \prod _{q=0}^2}(t-y_q),\hfill \end{array}$

where $({\xi }_1(s),{\eta }_1(t))\in [a,x_1]\times [c,y_1]$. For all $(s,t)\in [a,x_1]\times [y_{2j-1},y_{2j+1}]$, there exists $({\xi }_2(s),{\eta }_j(t))\in [a,x_1]\times [y_{2j-1},y_{2j+1}]$, such that

$\begin{array}{ccc}\hfill R_2& =& \frac{1}{3!}{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_2(s),t,u({\xi }_2(s),t)){\displaystyle \prod _{k=0}^2}(s-x_k)\hfill \\ & & +{\displaystyle \sum _{k=0}^2}\frac{{\phi }_{k,0}(s)}{3!}{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_k,{\eta }_j(t),u(x_k,{\eta }_j(t))){\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q}).\hfill \end{array}$

Similarly, for all $(s,t)\in [x_{2i-1},x_{2i+1}]\times [c,y_1]$, there exists $({\xi }_i(s),{\eta }_2(t))\in [x_{2i-1},x_{2i+1}]\times [c,y_1]$, such that

$\begin{array}{ccc}\hfill R_3& =& \frac{1}{3!}{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_i(s),t,u({\xi }_i(s),t)){\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})\hfill \\ & & +{\displaystyle \sum _{k=0}^2}\frac{{\phi }_{k,2i-1}(s)}{3!}{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_{2i-1+k},{\eta }_2(t),u(x_{2i-1+k},{\eta }_2(t))){\displaystyle \prod _{q=0}^2}(t-y_q),\hfill \end{array}$

and for all $(s,t)\in [x_{2i-1},x_{2i+1}]\times [y_{2j-1},y_{2j+1}]$, there exists $({\xi }_{i1}(s),{\eta }_{j2}(t))\in [x_{2i-1},x_{2i+1}]\times [y_{2j-1},y_{2j+1}]$, such that

$\begin{array}{ccc}\hfill R_4& =& \frac{1}{3!}{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_{i1}(s),t,u({\xi }_{i1}(s),t)){\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})\hfill \\ & & +{\displaystyle \sum _{k=0}^2}\frac{{\phi }_{k,2i-1}(s)}{3!}{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_{2i-1+k},{\eta }_{j2}(t),u(x_{2i-1+k},{\eta }_{j2}(t))){\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q}).\hfill \end{array}$

We can directly obtain $r_{2m+1,2n+1}^{(1)}$ as follows

(38)$\begin{array}{ccc}\hfill & & |r_{2m+1,2n+1}^{(1)}|\le |{\int }_a^{x_1}{\int }_c^{y_1}\frac{{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_1(s),t,u({\xi }_1(s),t)){\displaystyle \prod _{k=0}^2}(s-x_k)}{{3!(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds|\hfill \\ \hfill & & +|{\int }_a^{x_1}{\int }_c^{y_1}{\displaystyle \sum _{k=0}^2}{\phi }_{k,0}(s)\frac{{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_k,{\eta }_1(t),u(x_k,{\eta }_1(t)))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{q=0}^2}(t-y_q)dtds|\hfill \\ & & \doteq R_1^{(1)}+R_1^{(2)}.\hfill \end{array}$

Next, we estimate each item on the right side of (38), and we have

(39)$\begin{array}{ccc}\hfill R_1^{(1)}& \le & G_1{h}_x^3|{\int }_a^{x_1}{\int }_c^{y_1}{(x_{2m+1}-s)}^{-\alpha }{(y_{2n+1}-t)}^{-\beta }dtds|\hfill \\ \hfill & =& G_1{h}_x^{4-\alpha }{h}_y^{1-\beta }\frac{\left({(2m+1)}^{1-\alpha }-{(2m)}^{1-\alpha }\right)\left({(2n+1)}^{1-\beta }-{(2n)}^{1-\beta }\right)}{(1-\alpha )(1-\beta )}\hfill \\ & =& G_1{h}_x^{4-\alpha }{h}_y^{1-\beta }\frac{1}{{\varsigma }^{\alpha }}\frac{1}{{\xi }^{\beta }}\le G_1{h}_x^{4-\alpha }{h}_y^{1-\beta },\varsigma \in (2m,2m+1),\xi \in (2n,2n+1),\hfill \end{array}$

(40)$\begin{array}{ccc}\hfill R_1^{(2)}& \le & G_1{h}_y^3|{\int }_a^{x_1}{\int }_c^{y_1}{(x_{2m+1}-s)}^{-\alpha }{(y_{2n+1}-t)}^{-\beta }dtds|\le G_1{h}_x^{1-\alpha }{h}_y^{4-\beta },\hfill \end{array}$

where $G_1$ is defined by (35).

Combining (39) and (40) yields

(41)$\begin{array}{c}\hfill |r_{2m+1,2n+1}^{(1)}|\le G_1(h_x^{4-\alpha }{h}_y^{1-\beta }+h_x^{1-\alpha }{h}_y^{4-\beta }).\end{array}$

We use the same technique for $r_{2m+1,2n+1}^{(2)}$, and it is easy to see that

(42)$\begin{array}{cc}\hfill & |r_{2m+1,2n+1}^{(2)}|\le {\displaystyle \sum _{j=1}^n}|{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_2(s),t,u({\xi }_2(s),t))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{k=0}^2}(s-x_k)dtds|\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^n}|{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\displaystyle \sum _{k=0}^2}{\phi }_{k,0}(s){\partial }_t^3{\kappa }_{2m+1,2n+1}(x_k,{\eta }_j(t),u(x_k,{\eta }_j(t))){\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds|\hfill \\ \hfill & \dot{=}{R}_2^{(1)}+R_2^{(2)}.\hfill \end{array}$

For the first item on the right side of (42), we can obtain

(43)$\begin{array}{cc}\hfill R_2^{(1)}& \le G_1{h}_x^3{\displaystyle \sum _{j=1}^n}|{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}{(x_{2m+1}-s)}^{-\alpha }{(y_{2n+1}-t)}^{-\beta }dtds|\hfill \\ \hfill & \le G_1{h}_x^{4-\alpha }{\displaystyle \sum _{j=1}^n}{\int }_{y_{2j-1}}^{y_{2j+1}}{(y_{2n+1}-t)}^{-\beta }dt\le G_1{h}_x^{4-\alpha }{\int }_{y_1}^{y_{2n+1}}{(y_{2n+1}-t)}^{-\beta }dt\hfill \\ \hfill & =G_1{h}_x^{4-\alpha }\frac{{((2n)h_y)}^{1-\beta }}{1-\beta }\le \frac{d^{1-\beta }}{1-\beta }{G}_1{h}_x^{4-\alpha }.\hfill \end{array}$

Next, we make a detailed estimate of $R_2^{(2)}$ with the following form

(44)$\begin{array}{cc}\hfill & R_2^{(2)}\le \hfill \\ \hfill & {\displaystyle \sum _{j=1}^n}|{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\displaystyle \sum _{k=0}^2}{\phi }_{k,0}(s){\partial }_t^3{\kappa }_{2m+1,2n+1}(x_k,{\eta }_j({\tilde{t}}_j),u(x_k,{\eta }_j({\tilde{t}}_j))){\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds|\hfill \\ \hfill & +{\displaystyle \sum _{j=1}^n}|{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{\displaystyle \sum _{k=0}^2}{\phi }_{k,0}(s)[\frac{1}{3!}{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_k,{\eta }_j(t),u(x_k,{\eta }_j(t)))\hfill \\ \hfill & -\frac{1}{3!}{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_k,{\eta }_j({\tilde{t}}_j),u(x_k,{\eta }_j({\tilde{t}}_j)))]{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})dtds|\hfill \\ \hfill & \dot{=}{A}_1+A_2,\hfill \end{array}$

where ${\tilde{t}}_j=y_{2j}$. Based on the definition of ${\phi }_{k,i}(s),s\in (x_i,x_{i+2}),k=0,1,2;i\in \mathbb{N}$, we know that $|{\phi }_{k,i}(s)|\le 1$, and therefore $|{\displaystyle \sum _{k=0}^2}{\phi }_{k,i}(s)|\le 3$, and we have

(45)$\begin{array}{ccc}\hfill A_1& \le & G_1{\displaystyle \sum _{j=1}^n}|{\int }_a^{x_1}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{\displaystyle \sum _{k=0}^2}{\phi }_{k,0}(s){\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})dtds|\hfill \\ \hfill & =& G_1|{\int }_a^{x_1}\frac{{\displaystyle \sum _{k=0}^2}{\phi }_{k,0}(s)}{{(x_{2m+1}-s)}^{\alpha }}ds|\times {\displaystyle \sum _{j=1}^n}|{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt|\hfill \\ \hfill & \le & 3G_1{h}_x^{1-\alpha }{\displaystyle \sum _{j=1}^n}|{\int }_{y_{2j-1}}^{y_{2j}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt+{\int }_{y_{2j}}^{y_{2j+1}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt|\hfill \\ \hfill & \le & 3G_1{h}_x^{1-\alpha }{\displaystyle \sum _{j=1}^{n-2}}|{\int }_{y_{2j-1}}^{y_{2j}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt+{\int }_{y_{2j}}^{y_{2j+1}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt|\hfill \\ \hfill & & +3G_1{h}_x^{1-\alpha }\left(|{\int }_{y_{2n-3}}^{y_{2n-1}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2n-3+q})}{{(y_{2n+1}-t)}^{\beta }}dt|+|{\int }_{y_{2n-1}}^{y_{2n+1}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2n-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt|\right)\hfill \\ & \doteq & 3G_1{h}_x^{1-\alpha }{A}_{11}+3G_1{h}_x^{1-\alpha }{A}_{12}.\hfill \end{array}$

We make the following detailed estimates of $A_{11}$ and $A_{12}$, which can be directly obtained.

(46)$\begin{array}{ccc}\hfill A_{11}& =& {\displaystyle \sum _{j=1}^{n-2}}|{(y_{2n+1}-{\widehat{t}}_j)}^{-\beta }{\int }_{y_{2j-1}}^{y_{2j}}{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})dt\hfill \\ \hfill & & +{(y_{2n+1}-{\overline{t}}_j)}^{-\beta }{\int }_{y_{2j}}^{y_{2j+1}}{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})dt|\hfill \\ \hfill & =& {\displaystyle \sum _{j=1}^{n-2}}|\frac{1}{4}{(y_{2n+1}-{\widehat{t}}_j)}^{-\beta }{h}_y^4-\frac{1}{4}{(y_{2n+1}-{\overline{t}}_j)}^{-\beta }{h}_y^4|\hfill \\ \hfill & =& \frac{1}{4}{h}_y^4{\displaystyle \sum _{j=1}^{n-2}}|{(y_{2n+1}-{\widehat{t}}_j)}^{-\beta }-{(y_{2n+1}-{\overline{t}}_j)}^{-\beta }|\hfill \\ \hfill & =& \frac{1}{4}{h}_y^4{\displaystyle \sum _{j=1}^{n-2}}|-\beta {(y_{2n+1}-t_j)}^{-\beta -1}(\overline{t_j}-\widehat{t_j})|\le \frac{\beta }{4}{h}_y^4{\displaystyle \sum _{j=1}^{n-2}}|2{(y_{2n+1}-y_{2j+1})}^{-\beta -1}{h}_y|\hfill \\ \hfill & \le & \frac{\beta }{4}{h}_y^4{\int }_{y_3}^{y_{2n-1}}{(y_{2n+1}-t)}^{-\beta -1}dt=\frac{1}{4}{h}_y^4\left(2^{-\beta }{h}_y^{-\beta }-{(2n-2)}^{-\beta }{h}_y^{-\beta }\right)\hfill \\ & \le & \frac{2}{4}{h}_y^4{h}_y^{-\beta }=\frac{1}{2}{h}_y^{4-\beta },\hfill \end{array}$

(47)$\begin{array}{ccc}\hfill A_{12}& \le & h_y^3\left({\int }_{y_{2n-3}}^{y_{2n-1}}{(y_{2n+1}-t)}^{-\beta }dt+{\int }_{y_{2n-1}}^{y_{2n+1}}{(y_{2n+1}-t)}^{-\beta }dt\right)\hfill \\ & =& h_y^3{\int }_{y_{2n-3}}^{y_{2n+1}}{(y_{2n+1}-t)}^{-\beta }dt=\frac{4^{1-\beta }}{1-\beta }{h}_y^{4-\beta },\hfill \end{array}$

where ${\widehat{t}}_j\le t_j\le {\overline{t}}_j,y_{2j-1}\le {\widehat{t}}_j\le y_{2j}\le {\overline{t}}_j\le y_{2j+1}$. Then, we insert (46) and (47) into (45) and obtain $A_1$ as follows

(48)$\begin{array}{c}\hfill A_1\le 3\left(\frac{1}{2}+\frac{4^{1-\beta }}{1-\beta }\right)G_1{h}_x^{1-\alpha }{h}_y^{4-\beta }.\end{array}$

For $A_2$, we discover

(49)$\begin{array}{cc}\hfill & A_2\le G_2{h}_y{\int }_a^{x_1}\frac{|{\displaystyle \sum _{k=0}^2}{\phi }_{k,0}(s)|}{3!{(x_{2m+1}-s)}^{\alpha }}ds\times {\displaystyle \sum _{j=1}^n}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{|{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})|}{{(y_{2n+1}-t)}^{\beta }}dt,\hfill \\ \hfill & \le G_2{h}_x^{1-\alpha }{h}_y^4{\displaystyle \sum _{j=1}^n}{\int }_{y_{2j-1}}^{y_{2j+1}}{(y_{2n+1}-t)}^{-\beta }dt\le \frac{G_2{d}^{1-\beta }}{(1-\beta )}{h}_x^{1-\alpha }{h}_y^4,\hfill \end{array}$

where $G_2$ is defined by ().

According to (48) and (49), (44) becomes

(50)$\begin{array}{c}\hfill R_2^{(2)}\le 3\left(\frac{1}{2}+\frac{4^{1-\beta }}{1-\beta }\right)G_1{h}_x^{1-\alpha }{h}_y^{4-\beta }+\frac{G_2{d}^{1-\beta }}{(1-\beta )}{h}_x^{1-\alpha }{h}_y^4.\end{array}$

Therefore, combining (43) and (50) with (42), we can obtain the following result:

(51)$\begin{array}{cc}\hfill & |r_{2m+1,2n+1}^{(2)}|\le \frac{d^{1-\beta }}{1-\beta }{G}_1{h}_x^{4-\alpha }+\left[3\left(\frac{1}{2}+\frac{4^{1-\beta }}{1-\beta }\right)G_1{h}_y^{4-\beta }+\frac{G_2{d}^{1-\beta }}{1-\beta }{h}_y^4\right]h_x^{1-\alpha }.\hfill \end{array}$

Now, we estimate $r_{2m+1,2n+1}^{(3)}$, and we can obtain

$\begin{array}{cc}\hfill & |r_{2m+1,2n+1}^{(3)}|\le |{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_i(s),t,u({\xi }_i(s),t))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})dtds|\hfill \\ \hfill & +|{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}{\displaystyle \sum _{k=0}^2}{\phi }_{k,2i-1}(s)\frac{{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_{2i-1+k},{\eta }_2(t),u(x_{2i-1+k},{\eta }_2(t)))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{q=0}^2}(t-y_q)dtds|\hfill \\ \hfill & \dot{=}{R}_3^{(1)}+R_3^{(2)}.\hfill \end{array}$

For $R_3^{(1)}$, we have

$\begin{array}{cc}\hfill & R_3^{(1)}\le {\displaystyle \sum _{i=1}^m}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_i({\tilde{s}}_i),t,u({\xi }_i({\tilde{s}}_i),t))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})dtds|\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^m}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{3!{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})\times \hfill \\ \hfill & ({\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_i(s),t,u({\xi }_i(s),t))-{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_i({\tilde{s}}_i),t,u({\xi }_i({\tilde{s}}_i),t)))dtds|,\hfill \end{array}$

where ${\tilde{s}}_i=x_{2i}$. Taking the first term and the second term on the right side of the above formula as $B_1$ and $B_2$, we have

$\begin{array}{ccc}\hfill B_1& \le & G_1{h}_y^{1-\beta }{\displaystyle \sum _{i=1}^m}|{\int }_{x_{2i-1}}^{x_{2i}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})}{{(x_{2m+1}-s)}^{\alpha }}ds+{\int }_{x_{2i}}^{x_{2i+1}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})}{{(x_{2m+1}-s)}^{\alpha }}ds|\hfill \\ & \le & G_1{h}_y^{1-\beta }{\displaystyle \sum _{i=1}^{m-2}}|{\int }_{x_{2i-1}}^{x_{2i}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})}{{(x_{2m+1}-s)}^{\alpha }}ds+{\int }_{x_{2i}}^{x_{2i+1}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})}{{(x_{2m+1}-s)}^{\alpha }}ds|\hfill \\ & & +G_1{h}_y^{1-\beta }\left(|{\int }_{x_{2m-3}}^{x_{2m-1}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2m-3+k})}{{(x_{2m+1}-s)}^{\alpha }}ds|+|{\int }_{x_{2m-1}}^{x_{2m+1}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2m-1+k})}{{(x_{2m+1}-s)}^{\alpha }}ds|\right)\hfill \\ & \doteq & G_1{h}_y^{1-\beta }{B}_{11}+G_1{h}_y^{1-\beta }{B}_{12}.\hfill \end{array}$

For $B_{11}$ and $B_{12}$, we use the same method as (46) and (47), obtaining

$\begin{array}{c}\hfill B_{11}\le \frac{1}{2}{h}_x^{4-\alpha },B_{12}\le \frac{4^{1-\alpha }}{1-\alpha }{h}_x^{4-\alpha }.\end{array}$

Hence, for $B_1$, we can directly obtain

$\begin{array}{c}\hfill B_1\le \left(\frac{1}{2}+\frac{4^{1-\alpha }}{1-\alpha }\right)G_1{h}_x^{4-\alpha }{h}_y^{1-\beta }.\end{array}$

For $B_2$, we can directly obtain the result

$\begin{array}{c}\hfill B_2\le G_2{h}_x^4{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}dtds\le \frac{b^{1-\alpha }}{1-\alpha }{G}_2{h}_x^4{h}_y^{1-\beta }.\end{array}$

Next, we estimate $R_3^{(2)}$ as follows

$\begin{array}{ccc}\hfill R_3^{(2)}& \le & G_1{\displaystyle \sum _{i=1}^m}|{\int }_{x_{2i-1}}^{x_{2i+1}}\frac{{\displaystyle \sum _{k=0}^2}{\phi }_{k,2i-1}(s)}{3!{(x_{2m+1}-s)}^{\alpha }}ds|·|{\int }_0^{y_1}\frac{{\displaystyle \prod _{q=0}^2}(t-y_q)}{{(y_{2n+1}-t)}^{\beta }}dt|\hfill \\ & \le & G_1{h}_y^3{\displaystyle \sum _{i=1}^m}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_c^{y_1}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}dtds\le \frac{b^{1-\alpha }{G}_1}{(1-\alpha )}{h}_y^{4-\beta }.\hfill \end{array}$

Therefore, we obtain that

(52)$\begin{array}{c}\hfill |r_{2m+1,2n+1}^{(3)}|\le \left[\left(\frac{1}{2}+\frac{4^{1-\alpha }}{1-\alpha }\right)G_1{h}_x^{4-\alpha }+\frac{b^{1-\alpha }}{1-\alpha }{G}_2{h}_x^4\right]h_y^{1-\beta }+\frac{b^{1-\alpha }{G}_1}{1-\alpha }{h}_y^{4-\beta }.\end{array}$

Finally, we make a detailed estimation of $r_{2m+1,2n+1}^{(4)}$, which has the form

$\begin{array}{cc}\hfill & |r_{2m+1,2n+1}^{(4)}|\le \hfill \\ & {\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_{i1}(s),t,u({\xi }_{i1}(s),t))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})dtds|\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{3!{(y_{2n+1}-t)}^{\beta }}{\displaystyle \sum _{k=0}^2}{\phi }_{k,2i-1}(s)\hfill \\ \hfill & \times {\partial }_t^3{\kappa }_{2m+1,2n+1}(x_{2i-1+k},{\eta }_{j2}(t),u(x_{2i-1+k},{\eta }_{j2}(t))){\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})dtds|\hfill \\ \hfill & \dot{=}{R}_4^{(1)}+R_4^{(2)}.\hfill \end{array}$

Next, we estimate $R_4^{(1)}$ and $R_4^{(2)}$ one by one. For $R_4^{(1)}$, we have

$\begin{array}{ccc}& & R_4^{(1)}\le \hfill \\ & & {\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_{i1}({\tilde{s}}_i),t,u({\xi }_{i1}({\tilde{s}}_i),t)){\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}dtds|\hfill \\ & & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{3!{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})\hfill \\ & & \times ({\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_{i1}(s),t,u({\xi }_{i1}(s),t))-{\partial }_s^3{\kappa }_{2m+1,2n+1}({\xi }_{i1}({\tilde{s}}_i),t,u({\xi }_{i1}({\tilde{s}}_i),t)))dtds|\hfill \\ & \doteq & P_1+P_2,\hfill \end{array}$

where ${\tilde{s}}_i=x_{2i}$.

For $P_1$, the same processing method can be used as for $B_1$, to obtain

$\begin{array}{ccc}\hfill P_1& \le & G_1{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})dtds|\hfill \\ & \le & \frac{d^{1-\beta }{G}_1}{1-\beta }{\displaystyle \sum _{i=1}^m}|{\int }_{x_{2i-1}}^{x_{2i}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})}{{(x_{2m+1}-s)}^{\alpha }}ds+{\int }_{x_{2i}}^{x_{2i+1}}\frac{{\displaystyle \prod _{k=0}^2}(s-x_{2i-1+k})}{{(x_{2m+1}-s)}^{\alpha }}ds|\hfill \\ & \le & \frac{d^{1-\beta }}{1-\beta }\left(\frac{1}{2}+\frac{4^{1-\alpha }}{1-\alpha }\right)G_1{h}_x^{4-\alpha }.\hfill \end{array}$

For $P_2$, we can obtain

$\begin{array}{c}\hfill P_2\le G_2{h}_x^4{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{{(y_{2n+1}-t)}^{\beta }}dtds\le \frac{b^{1-\alpha }{d}^{1-\beta }}{(1-\alpha )(1-\beta )}{G}_2{h}_x^4.\end{array}$

Finally, we estimate $R_4^{(2)}$ as follows

$\begin{array}{cc}\hfill & R_4^{(2)}\le {\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{(x_{2m+1}-s)}^{-\alpha }}{3!{(y_{2n+1}-t)}^{\beta }}{\displaystyle \sum _{k=0}^2}{\phi }_{k,2i-1}(s)\times \hfill \\ \hfill & {\partial }_t^3{\kappa }_{2m+1,2n+1}(x_{2i-1+k},{\eta }_{j2}({\tilde{t}}_j)),u(x_{2i-1+k},{\eta }_{j2}({\tilde{t}}_j)){\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})dtds|\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^m}{\displaystyle \sum _{j=1}^n}|{\int }_{x_{2i-1}}^{x_{2i+1}}{\int }_{y_{2j-1}}^{y_{2j+1}}{\displaystyle \sum _{k=0}^2}{\phi }_{k,2i-1}(s)[\frac{{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_{2i-1+k},{\eta }_{j2}(t)),u(x_{2i-1+k},{\eta }_{j2}(t))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}\hfill \\ \hfill & -\frac{{\partial }_t^3{\kappa }_{2m+1,2n+1}(x_{2i-1+k},{\eta }_{j2}({\tilde{t}}_j)),u(x_{2i-1+k},{\eta }_{j2}({\tilde{t}}_j))}{3!{(x_{2m+1}-s)}^{\alpha }{(y_{2n+1}-t)}^{\beta }}]{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})dtds|,\hfill \end{array}$

where ${\tilde{t}}_j=y_{2j}$. The two terms of the above formula are denoted $J_1$ and $J_2$, respectively. Using the same method as for $A_1$ to deal with $J_1$, we can obtain

$\begin{array}{ccc}\hfill J_1& \le & G_1{\displaystyle \sum _{i=1}^m}|{\int }_{x_{2i-1}}^{x_{2i+1}}\frac{{\displaystyle \sum _{k=0}^2}{\phi }_{k,2i-1}(s)}{3!{(x_{2m+1}-s)}^{\alpha }}ds|\times {\displaystyle \sum _{j=1}^n}|{\int }_{y_{2j-1}}^{y_{2j+1}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt|\hfill \\ & \le & \frac{G_1{b}^{1-\alpha }}{1-\alpha }{\displaystyle \sum _{j=1}^n}|{\int }_{y_{2j-1}}^{y_{2j}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt+{\int }_{y_{2j}}^{y_{2j+1}}\frac{{\displaystyle \prod _{q=0}^2}(t-y_{2j-1+q})}{{(y_{2n+1}-t)}^{\beta }}dt|\hfill \\ & \le & \frac{G_1{b}^{1-\alpha }}{1-\alpha }\left(\frac{1}{2}+\frac{4^{1-\beta }}{1-\beta }\right)h_y^{4-\beta },\hfill \\ \hfill J_2& \le & G_2{h}_y^4{\displaystyle \sum _{i=1}^m}|{\int }_{x_{2i-1}}^{x_{2i+1}}\frac{{\displaystyle \sum _{k=0}^2}{\phi }_{k,2i-1}(s)}{3!{(x_{2m+1}-s)}^{\alpha }}ds|\times {\displaystyle \sum _{j=1}^n}{\int }_{y_{2j-1}}^{y_{2j+1}}{(y_{2n+1}-t)}^{-\beta }dt\hfill \\ & \le & \frac{G_2{b}^{1-\alpha }{d}^{1-\beta }}{(1-\alpha )(1-\beta )}{h}_y^4.\hfill \end{array}$

Therefore, we obtain

(53)$\begin{array}{ccc}\hfill |r_{2m+1,2n+1}^{(4)}|& \le & G_1\left[\left(\frac{1}{2}+\frac{4^{1-\alpha }}{1-\alpha }\right)\frac{d^{1-\beta }}{1-\beta }{h}_x^{4-\alpha }+(\frac{1}{2}+\frac{4^{1-\beta }}{1-\beta })\frac{b^{1-\alpha }}{1-\alpha }{h}_y^{4-\beta }\right]\hfill \\ & & +\frac{b^{1-\alpha }{d}^{1-\beta }{G}_2}{(1-\alpha )(1-\beta )}(h_x^4+h_y^4).\hfill \end{array}$

Bringing (41) and (51)–(53) into (37) gives

$\begin{array}{cc}\hfill & |r_{2m+1,2n+1}|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),\hfill \end{array}$

where C depends only on $G_1,G_2,\alpha$ and $\beta$.

We can prove the truncation error at the other steps in a similar way to the truncation error at the point $(x_{2m+1},y_{2n+1})$. We can draw the conclusion that the truncation errors $r_{i,j}$ satisfy

(54)$\begin{array}{c}\hfill |r_{i,j}|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),i=1,2,\cdots ,2M;j=1,2,\cdots ,2N.\end{array}$

The proof is then completed. □

4. Convergence Analysis

To simplify the notation, we rewrite the numerical scheme by introducing the following coefficients

(55)$\begin{array}{ccc}\hfill & & {\widehat{B}}_i=\frac{{\omega }_1^{i,0}}{h_x^{1-\alpha }},{\tilde{B}}_i=\frac{A_2^{i,0}}{h_x^{1-\alpha }},i=0,1,2;{\overline{B}}_0^m=\frac{{\omega }_{2m+1}^{0,0}}{h_x^{1-\alpha }},{\overline{B}}_1^m=\frac{{\omega }_{2m+1}^{1,0}+A_{2m+1}^{0,1}}{h_x^{1-\alpha }},\hfill \\ \hfill & & {\overline{B}}_2^m=\frac{{\omega }_{2m+1}^{2,0}+A_{2m+1}^{1,1}}{h_x^{1-\alpha }},B_0^m=\frac{A_{2m+2}^{0,0}}{h_x^{1-\alpha }},B_{2m+2}^m=\frac{A_{2m+2}^{2,m}}{h_x^{1-\alpha }},\hfill \\ \hfill & & B_{2i+1}^m=\frac{A_{2m+2}^{1,i}}{h_x^{1-\alpha }},i=0,1,\cdots ,m;B_{2i}^m=\frac{A_{2m+2}^{2,i-1}+A_{2m+2}^{0,i}}{h_x^{1-\alpha }},i=1,2,\cdots ,m;\hfill \\ \hfill & & {\widehat{D}}_j=\frac{{\widehat{\omega }}_1^{j,0}}{h_y^{1-\beta }},{\tilde{D}}_j=\frac{{\widehat{A}}_2^{j,0}}{h_y^{1-\beta }},j=0,1,2;{\overline{D}}_0^n=\frac{{\widehat{\omega }}_{2n+1}^{0,0}}{h_y^{1-\beta }},{\overline{D}}_1^n=\frac{{\widehat{\omega }}_{2n+1}^{1,0}+{\widehat{A}}_{2n+1}^{0,1}}{h_y^{1-\beta }},\hfill \\ \hfill & & {\overline{D}}_2^n=\frac{{\widehat{\omega }}_{2n+1}^{2,0}+{\widehat{A}}_{2n+1}^{1,1}}{h_y^{1-\beta }},D_0^n=\frac{{\widehat{A}}_{2n+2}^{0,0}}{h_y^{1-\beta }},D_{2n+2}^n=\frac{{\widehat{A}}_{2n+2}^{2,n}}{h_y^{1-\beta }},\hfill \\ \hfill & & D_{2j+1}^n=\frac{{\widehat{A}}_{2n+2}^{1,j}}{h_y^{1-\beta }},j=0,1,\cdots ,n;D_{2j}^n=\frac{{\widehat{A}}_{2n+2}^{2,j-1}+{\widehat{A}}_{2n+2}^{0,j}}{h_y^{1-\beta }},j=1,2,\cdots ,n.\hfill \end{array}$

Carefully observing (11), (15), (21) and (22), it is obvious that $A_{2m+1}^{k,i}=A_{2m+2}^{k,i},k=0,1,2,i=1,2,\cdots ,m;$${\widehat{A}}_{2n+1}^{q,j}={\widehat{A}}_{2n+2}^{q,j},q=0,1,2,j=1,2,\cdots ,n$. Hence, the numerical scheme (32) can be reformulated with an equivalent form as follows

(56)$\begin{array}{cc}\hfill & u_{1,1}=f_{1,1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\widehat{B}}_i{\widehat{D}}_j{\kappa }_{1,1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2,1}=f_{2,1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\tilde{B}}_i{\widehat{D}}_j{\kappa }_{2,1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{1,2}=f_{1,2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\widehat{B}}_i{\tilde{D}}_j{\kappa }_{1,2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2,2}=f_{2,2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\tilde{B}}_i{\tilde{D}}_j{\kappa }_{2,2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+1,1}=f_{2m+1,1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\overline{B}}_i^m{\widehat{D}}_j{\kappa }_{2m+1,1}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^2}{B}_{i+1}^m{\widehat{D}}_j{\kappa }_{2m+1,1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+2,1}=f_{2m+2,1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^2}{B}_i^m{\widehat{D}}_j{\kappa }_{2m+2,1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+1,2}=f_{2m+1,2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\overline{B}}_i^m{\tilde{D}}_j{\kappa }_{2m+1,2}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^2}{B}_{i+1}^m{\tilde{D}}_j{\kappa }_{2m+1,2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+2,2}=f_{2m+2,2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^2}{B}_i^m{\tilde{D}}_j{\kappa }_{2m+2,2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{1,2n+1}=f_{1,2n+1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\widehat{B}}_i{\overline{D}}_j^n{\kappa }_{1,2n+1}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=3}^{2n+1}}{\widehat{B}}_i{D}_{j+1}^n{\kappa }_{1,2n+1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2,2n+1}=f_{2,2n+1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\tilde{B}}_i{\overline{D}}_j^n{\kappa }_{2,2n+1}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=3}^{2n+1}}{\tilde{B}}_i{D}_{j+1}^n{\kappa }_{2,2n+1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{1,2n+2}=f_{1,2n+2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^{2n+2}}{\widehat{B}}_i{D}_j^n{\kappa }_{1,2n+2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2,2n+2}=f_{2,2n+2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^{2n+2}}{\tilde{B}}_i{D}_j^n{\kappa }_{2,2n+2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+1,2n+1}=f_{2m+1,2n+1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\overline{B}}_i^m{\overline{D}}_j^n{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^2}{B}_{i+1}^m{\overline{D}}_j^n{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=3}^{2n+1}}{\overline{B}}_i^m{D}_{j+1}^n{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=3}^{2n+1}}{B}_{i+1}^m{D}_{j+1}^n{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+2,2n+1}=f_{2m+2,2n+1}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^2}{B}_i^m{\overline{D}}_j^n{\kappa }_{2m+2,2n+1}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=3}^{2n+1}}{B}_i^m{D}_{j+1}^n{\kappa }_{2m+2,2n+1}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+1,2n+2}=f_{2m+1,2n+2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^{2n+2}}{\overline{B}}_i^m{D}_j^n{\kappa }_{2m+1,2n+2}(x_i,y_j,u_{i,j})\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^{2n+2}}{B}_{i+1}^m{D}_j^n{\kappa }_{2m+1,2n+2}(x_i,y_j,u_{i,j}),\hfill \\ \hfill & u_{2m+2,2n+2}=f_{2m+2,2n+2}+h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^{2n+2}}{B}_i^m{D}_j^n{\kappa }_{2m+2,2n+2}(x_i,y_j,u_{i,j}),\hfill \end{array}$

where $m=1,2,\cdots ,M-1$; $n=1,2,\cdots ,N-1$.

Based on the idea of [5], we can prove the following lemma.

Lemma 2.

The coefficients ${\overline{B}}_i^m,B_i^m$ and ${\overline{D}}_j^n,D_j^n$, are defined in (55), and they satisfy

$\begin{array}{cc}\hfill & |{\overline{B}}_i^m|\le C{(2m+2-i)}^{-\alpha },i=0,1,2,|B_i^m|\le C{(2m+3-i)}^{-\alpha },i=0,1,\cdots ,2m+2,\hfill \\ \hfill & |{\overline{D}}_j^n|\le C{(2n+2-j)}^{-\beta },j=0,1,2,|D_j^n|\le C{(2n+3-j)}^{-\beta },j=0,1,\cdots ,2n+2.\hfill \end{array}$

In the following, we will use Lemmas 1 and 2 to analyze the convergence of the scheme (56), as in the following theorem.

Theorem 1.

Let u be the exact solution of Equation (1), $u_{i,j},i=0,1,\cdots ,2M$, $j=0,1,\cdots ,2N$, be the numerical solution of (56). If $\kappa (·,·,·,·,u(·,·))\in C^4([a,b]\times [c,d]\times R)$ satisfies the condition (2), and the step size $h_x$, $h_y$ satisfies

(57) $\begin{array}{c}\hfill L{h}_x^{1-\alpha }{h}_y^{1-\beta }\left|B_{2m+2}^m\right|\left|D_{2n+2}^n\right|<1,2L{h}_x^{1-\alpha }{h}_y^{1-\beta }\left|B_{2m+2}^m\right|<1,2L{h}_x^{1-\alpha }{h}_y^{1-\beta }\left|D_{2n+2}^n\right|<1,CL{h}_x^{1-\alpha }\frac{d^{1-\beta }}{1-\beta }<1,\end{array}$

then the following error estimates hold

(58) $\begin{array}{c}\hfill |u(x_i,y_j)-u_{i,j}|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),i=1,2,\cdots ,2M;j=1,2,\cdots ,2N,\end{array}$

where C depends only on $L,\kappa ,b,d,\alpha$ and β.

Proof. 

Let $e_{i,j}=u(x_i,y_j)-u_{i,j},i=0,1,\cdots ,2M;j=0,1,\cdots ,2N$. It is easily seen that $e_{i,0}=e_{0,j}=0$. Firstly, $e_{i,j},i,j=1,2$, satisfy

$\begin{array}{cc}\hfill & e_{1,1}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\widehat{B}}_i{\widehat{D}}_j[{\kappa }_{1,1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{1,1}(x_i,y_j,u_{i,j})]+r_{1,1},\hfill \\ \hfill & e_{1,2}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\widehat{B}}_i{\tilde{D}}_j[{\kappa }_{1,2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{1,2}(x_i,y_j,u_{i,j})]+r_{1,2},\hfill \\ \hfill & e_{2,1}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\tilde{B}}_i{\widehat{D}}_j[{\kappa }_{2,1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2,1}(x_i,y_j,u_{i,j})]+r_{2,1},\hfill \\ \hfill & e_{2,2}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\tilde{B}}_i{\tilde{D}}_j[{\kappa }_{2,2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2,2}(x_i,y_j,u_{i,j})]+r_{2,2},\hfill \end{array}$

where ${\widehat{B}}_i,{\tilde{B}}_i,{\widehat{D}}_j$ and ${\tilde{D}}_j$ are defined in (55). By direct calculation, it is easy to show that ${\widehat{B}}_i,{\tilde{B}}_i,i=0,1,2$, and ${\widehat{D}}_j,{\tilde{D}}_j,j=0,1,2,$ are bounded and $\kappa$ satisfies (2). Hence, it follows that

$\begin{array}{cc}\hfill & |e_{1,1}|\le CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}|e_{i,j}|+|r_{1,1}|,|e_{1,2}|\le CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}|e_{i,j}|+|r_{1,2}|,\hfill \\ \hfill & |e_{2,1}|\le CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}|e_{i,j}|+|r_{2,1}|,|e_{2,2}|\le CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}|e_{i,j}|+|r_{2,2}|.\hfill \end{array}$

By combing these four inequalities, and taking into account (57), we discover that

$\begin{array}{cc}\hfill & |e_{i,j}|\le CL(|r_{1,1}|+|r_{1,2}|+|r_{2,1}|+|r_{2,2}|),i,j=1,2.\hfill \end{array}$

Secondly, $e_{i,j},i\ge 3;j=1,2,$ satisfies

$\begin{array}{cc}\hfill & e_{2m+1,1}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\overline{B}}_i^m{\widehat{D}}_j[{\kappa }_{2m+1,1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,1}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^2}{B}_{i+1}^m{\widehat{D}}_j[{\kappa }_{2m+1,1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,1}(x_i,y_j,u_{i,j})]+r_{2m+1,1},\hfill \\ \hfill & e_{2m+2,1}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^2}{B}_i^m{\widehat{D}}_j[{\kappa }_{2m+2,1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+2,1}(x_i,y_j,u_{i,j})]+r_{2m+2,1},\hfill \\ \hfill & e_{2m+1,2}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\overline{B}}_i^m{\tilde{D}}_j[{\kappa }_{2m+1,2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^2}{B}_{i+1}^m{\tilde{D}}_j[{\kappa }_{2m+1,2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2}(x_i,y_j,u_{i,j})]+r_{2m+1,2},\hfill \\ \hfill & e_{2m+2,2}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^2}{B}_i^m{\tilde{D}}_j[{\kappa }_{2m+2,2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+2,2}(x_i,y_j,u_{i,j})]+r_{2m+2,2},\hfill \end{array}$

where ${\overline{B}}_i^m$ and $B_i^m$ are in (55) and satisfy Lemma 2. Since the above four equations are coupled, they must be solved simultaneously. For simplicity, we choose to set $\left||{\widehat{e}}_i\right.∥=max\{|e_{i,1}|,|e_{i,2}|,i=0,1,\cdots ,2M\}$, $\left||{\widehat{r}}_i\right.∥=max\{|r_{i,1}|,|r_{i,2}|,i=0,1,\cdots ,2M\}$, and then we can convert them into two groups of similar formulas, one of which is

$\begin{array}{cc}\hfill ∥{\widehat{e}}_{2m+1}∥& \le 2CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥{\widehat{e}}_i∥+2L{h}_x^{1-\alpha }{h}_y^{1-\beta }|B_{2m+2}^m|∥{\widehat{e}}_{2m+1}∥+∥{\widehat{r}}_{2m+1}∥,\hfill \\ \hfill ∥{\widehat{e}}_{2m+2}∥& \le 2CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+1}}{(2m+3-i)}^{-\alpha }∥{\widehat{e}}_i∥+2L{h}_x^{1-\alpha }{h}_y^{1-\beta }|B_{2m+2}^m|∥{\widehat{e}}_{2m+2}∥+∥{\widehat{r}}_{2m+2}∥.\hfill \end{array}$

For $∥{\widehat{e}}_{2m+1}∥$, we have

$\begin{array}{cc}\hfill ∥{\widehat{e}}_{2m+1}∥& \le 2CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥{\widehat{e}}_i∥+C∥{\widehat{r}}_{2m+1}∥\hfill \\ \hfill & \le 2CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+1-i)}^{-\alpha }∥{\widehat{e}}_i∥+C∥{\widehat{r}}_{2m+1}∥.\hfill \end{array}$

After applying the Gronwall’s inequality [23] to the above equation, we obtain

$\begin{array}{cc}\hfill ∥{\widehat{e}}_{2m+1}∥& \le C∥{\widehat{r}}_{2m+1}∥E_{1-\alpha }(2CL{h}_y^{1-\beta }\Gamma (1-\alpha ){((2m+1)h_x)}^{1-\alpha })\hfill \\ \hfill & \le C∥{\widehat{r}}_{2m+1}∥E_{1-\alpha }(2CL{h}_y^{1-\beta }\Gamma (1-\alpha )b^{1-\alpha }).\hfill \end{array}$

Using the same method for $∥{\widehat{e}}_{2m+2}∥$, we directly obtain

$\begin{array}{c}\hfill ∥{\widehat{e}}_{2m+2}∥\le C∥{\widehat{r}}_{2m+2}∥E_{1-\alpha }(2CL{h}_y^{1-\beta }\Gamma (1-\alpha )b^{1-\alpha }).\end{array}$

The process of calculating $e_{i,j},i=1,2;j\ge 3,$ is similar to that for $e_{i,j},i\ge 3;j=1,2$, and therefore it is omitted.

Taking the above result, together with Lemma 1, we have

(59)$\begin{array}{cc}\hfill & |e_{i,j}|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),i\ge 1;j=1,2,\hfill \end{array}$

(60)$\begin{array}{cc}\hfill & |e_{i,j}|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),i=1,2;j\ge 3.\hfill \end{array}$

Next, $e_{i,j},i,j\ge 3$ satisfy

$\begin{array}{cc}\hfill & e_{2m+1,2n+1}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^2}{\overline{B}}_i^m{\overline{D}}_j^n[{\kappa }_{2m+1,2n+1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^2}{B}_{i+1}^m{\overline{D}}_j^n[{\kappa }_{2m+1,2n+1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=3}^{2n+1}}{\overline{B}}_i^m{D}_{j+1}^n[{\kappa }_{2m+1,2n+1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=3}^{2n+1}}{B}_{i+1}^m{D}_{j+1}^n[{\kappa }_{2m+1,2n+1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2n+1}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +r_{2m+1,2n+1},\hfill \\ \hfill & e_{2m+2,2n+1}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^2}{B}_i^m{\overline{D}}_j^n[{\kappa }_{2m+2,2n+1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+2,2n+1}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=3}^{2n+1}}{B}_i^m{D}_{j+1}^n[{\kappa }_{2m+2,2n+1}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+2,2n+1}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +r_{2m+2,2n+1},\hfill \\ \hfill & e_{2m+1,2n+2}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^2}{\displaystyle \sum _{j=0}^{2n+2}}{\overline{B}}_i^m{D}_j^n[{\kappa }_{2m+1,2n+2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2n+2}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=3}^{2m+1}}{\displaystyle \sum _{j=0}^{2n+2}}{B}_{i+1}^m{D}_j^n[{\kappa }_{2m+1,2n+2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+1,2n+2}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +r_{2m+1,2n+2},\hfill \\ \hfill & e_{2m+2,2n+2}=h_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m+2}}{\displaystyle \sum _{j=0}^{2n+2}}{B}_i^m{D}_j^n[{\kappa }_{2m+2,2n+2}(x_i,y_j,u(x_i,y_j))-{\kappa }_{2m+2,2n+2}(x_i,y_j,u_{i,j})]\hfill \\ \hfill & +r_{2m+2,2n+2},\hfill \end{array}$

where ${\overline{B}}_i^m,B_i^m,{\overline{D}}_j^n$ and $D_j^n$ are defined in (55) and satisfy Lemma 2. We estimate $e_{2m+1,2n+1}$ as shown below.

(61)$\begin{array}{cc}\hfill & \left|e_{2m+1,2n+1}\right|\le CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{\displaystyle \sum _{j=0}^{2n}}{(2m+2-i)}^{-\alpha }{(2n+2-j)}^{-\beta }|e_{i,j}|\hfill \\ \hfill & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{j=0}^{2n}}{(2n+2-j)}^{-\beta }|e_{2m+1,j}|+CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }|e_{i,2n+1}|\hfill \\ \hfill & +L{h}_x^{1-\alpha }{h}_y^{1-\beta }\left|B_{2m+2}^m\right|\left|D_{2n+2}^n\right|\left|e_{2m+1,2n+1}\right|+\left|r_{2m+1,2n+1}\right|,\hfill \end{array}$

which can be rearranged into the preferred form

$\begin{array}{ccc}\hfill \left|e_{2m+1,2n+1}\right|& \le & CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{\displaystyle \sum _{j=0}^{2n}}{\left(2m+2-i\right)}^{-\alpha }{\left(2n+2-j\right)}^{-\beta }\left|e_{i,j}\right|\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{j=0}^{2n}}{\left(2n+2-j\right)}^{-\beta }\left|e_{2m+1,j}\right|\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{\left(2m+2-i\right)}^{-\alpha }\left|e_{i,2n+1}\right|+C\left|r_{2m+1,2n+1}\right|.\hfill \end{array}$

If we let $∥e_i∥=\underset{0\le j\le 2N}{max}|e_{i,j}|$, $∥r_i∥=\underset{0\le j\le 2N}{max}|r_{i,j}|$, then we obtain

$\begin{array}{ccc}\hfill \left|e_{2m+1,2n+1}\right|& \le & CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{\displaystyle \sum _{j=0}^{2n}}{\left(2m+2-i\right)}^{-\alpha }{\left(2n+2-j\right)}^{-\beta }∥e_i∥\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{j=0}^{2n}}{\left(2n+2-j\right)}^{-\beta }∥e_{2m+1}∥\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{\left(2m+2-i\right)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥\hfill \\ & \le & CL{h}_x^{1-\alpha }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥{\displaystyle \sum _{j=0}^{2n}}({((2n+2-j)h_y)}^{-\beta }{h}_y)\hfill \\ & & +CL{h}_x^{1-\alpha }∥e_{2m+1}∥{\displaystyle \sum _{j=0}^{2n}}({((2n+2-j)h_y)}^{-\beta }{h}_y)\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥\hfill \\ & \le & CL{h}_x^{1-\alpha }{\displaystyle \sum _{i=0}^{2m}}{\left(2m+2-i\right)}^{-\alpha }∥e_i∥{\int }_{y_0}^{y_{2n+2}}{\left(y_{2n+2}-t\right)}^{-\beta }dt\hfill \\ & & +CL{h}_x^{1-\alpha }∥e_{2m+1}∥{\int }_{y_0}^{y_{2n+2}}{\left(y_{2n+2}-t\right)}^{-\beta }dt\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥\hfill \\ & =& CL{h}_x^{1-\alpha }\frac{y_{2n+2}^{1-\beta }}{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥+CL{h}_x^{1-\alpha }\frac{y_{2n+2}^{1-\beta }}{1-\beta }∥e_{2m+1}∥\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥\hfill \\ & \le & CL{h}_x^{1-\alpha }\frac{d^{1-\beta }}{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥+CL{h}_x^{1-\alpha }\frac{d^{1-\beta }}{1-\beta }∥e_{2m+1}∥\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥.\hfill \end{array}$

As the above formula is applicable to all $n=1,2,\cdots ,N-1$, we can obtain

$\begin{array}{ccc}\hfill ∥e_{2m+1}∥& \le & CL{h}_x^{1-\alpha }\frac{d^{1-\beta }}{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{\left(2m+2-i\right)}^{-\alpha }∥e_i∥+CL{h}_x^{1-\alpha }\frac{d^{1-\beta }}{1-\beta }∥e_{2m+1}∥\hfill \\ & & +CL{h}_x^{1-\alpha }{h}_y^{1-\beta }{\displaystyle \sum _{i=0}^{2m}}{\left(2m+2-i\right)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥.\hfill \end{array}$

Therefore,

(62)$\begin{array}{ccc}\hfill ∥e_{2m+1}∥& \le & \left(CL{h}_x^{1-\alpha }\frac{d^{1-\beta }}{1-\beta }+CL{h}_x^{1-\alpha }{h}_y^{1-\beta }\right){\displaystyle \sum _{i=0}^{2m}}{(2m+2-i)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥\hfill \\ & \le & \left(CL{h}_x^{1-\alpha }\frac{d^{1-\beta }}{1-\beta }+CL{h}_x^{1-\alpha }{h}_y^{1-\beta }\right){\displaystyle \sum _{i=0}^{2m}}{(2m+1-i)}^{-\alpha }∥e_i∥+C∥r_{2m+1}∥.\hfill \end{array}$

Using the discrete Gronwall inequality [23], we obtain

$\begin{array}{ccc}\hfill ∥e_{2m+1}∥& \le & C∥r_{2m+1}∥E_{1-\alpha }\left((CL\frac{d^{1-\beta }}{1-\beta }+CL{h}_y^{1-\beta )})\Gamma (1-\alpha ){((2m+1)h_x)}^{1-\alpha }\right)\hfill \\ \hfill & \le & C∥r_{2m+1}∥E_{1-\alpha }\left(CL\Gamma (1-\alpha )b^{1-\alpha }(\frac{d^{1-\beta }}{1-\beta }+h_y^{1-\beta )})\right).\hfill \end{array}$

Combining the above estimations with Lemma 1 yields

(63)$\begin{array}{c}\hfill \left|e_{2m+1,2n+1}\right|\le C(h_x^{4-\alpha }+h_y^{4-\beta }).\end{array}$

In the same way, we can obtain the results

(64)$\begin{array}{c}\hfill \left|e_{2m+2,2n+1}\right|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),\left|e_{2m+p,2n+2}\right|\le C(h_x^{4-\alpha }+h_y^{4-\beta }),p=1,2.\end{array}$

We combine (59), (60), (63) and (64), then the proof is completed. □

5. Numerical Examples

In this section, the numerical scheme is used to solve the two-dimensional fractional Volterra integral equations, and we propose four calculation examples to prove its effectiveness. All the numerical examples were implemented using MATLAB on a ThinkCentre computer with a 3.40 GHz Intel(R) Core(TM) i5-7500 CPU and 8.00 GB of RAM. The time complexity of the proposed numerical scheme on discrete grids was $O({(MN)}^2)$.

Example 1.

Consider the following two-dimensional linear fractional Volterra integral equations

$\begin{array}{cc}\hfill & u(x,y)=f(x,y)+{\int }_0^x{\int }_0^y\frac{(xy+s+t)u(s,t)}{{(x-s)}^{\alpha }{(y-t)}^{\beta }}dtds,(x,y)\in [0,1]\times [0,1],\hfill \end{array}$

where

$\begin{array}{cc}\hfill & f(x,y)=x^4{y}^4-576x^{6-\alpha }{y}^{6-\beta }{\displaystyle \prod _{i=1}^5}\frac{1}{(i-\alpha )}{\displaystyle \prod _{j=1}^5}\frac{1}{(j-\beta )}\hfill \\ \hfill & -2880x^{6-\alpha }{y}^{5-\beta }{\displaystyle \prod _{i=1}^6}\frac{1}{(i-\alpha )}{\displaystyle \prod _{j=1}^5}\frac{1}{(j-\beta )}-2880x^{5-\alpha }{y}^{6-\beta }{\displaystyle \prod _{i=1}^5}\frac{1}{(i-\alpha )}{\displaystyle \prod _{j=1}^6}\frac{1}{(j-\beta )}.\hfill \end{array}$

The corresponding exact solution is $u(x,y)=x^4{y}^4.$

In our test, we chose $\alpha =0.3,\beta =0.6$ and $\alpha =0.1,\beta =0.7$, respectively. The step size in the x-direction was $h_x=\frac{1}{2M}$ and the step size in the y-direction was $h_y=\frac{1}{2N}$. The definition of error was as follows

$\begin{array}{c}\hfill e_h^u=\underset{\begin{array}{c}i=1,\cdots ,2M\\ j=1,\cdots ,2N\end{array}}{max}\left|u(x_i,y_j)-u_{i,j}\right|.\end{array}$

In this example, we wished to check how well the constructed numerical format converged for different choices of α and β. The test results are given in Table 1 and Table 2, and the convergence order was calculated as $lo{g}_2(\frac{e_{2h}^u}{e_h^u})$. From the results of the theoretical analysis in Theorem 1, we know that the theoretical convergence order of the constructed numerical scheme is $O(h_x^{4-\alpha }+h_y^{4-\beta })$. When $h_x$ is sufficiently small with respect to $h_y$, we have $O(h_x^{4-\alpha }+h_y^{4-\beta })=O(h_y^{4-\beta })$. In Table 1 and Table 2, we take $\alpha =0.3,\beta =0.6$ or $\alpha =0.1,\beta =0.7$ and $M=2N$. In this case, our scheme’s theoretical order is $4-\beta$. In Table 1 and Table 2, we show the corresponding convergence order and CPU time when N takes a series of values and $\alpha ,\beta$ are given specific values. In Table 1, when $\alpha =0.3,\beta =0.6$, the tested order is close to the theoretical order of 3.4. In Table 2, when $\alpha =0.1,\beta =0.7$, the order of the test is close to the theoretical order of 3.3. This is consistent with the theoretical prediction. However, from Table 1 and Table 2, we find that the CPU time quickly increases from 0.14 s to about 2 h as N increases.

Similarly, we test the convergence order using the proposed numerical scheme in another way. When we take $N=\left[M^{\frac{4-\alpha }{4-\beta }}\right]$, where $[·]$ indicates rounding up, the numerical scheme’s theoretical order $O(h_x^{4-\alpha }+h_y^{4-\beta })=O(h_x^{4-\alpha })$. The test results are given in Table 3. When $\alpha =0.3,\beta =0.6$, the numerical result of the test is close to 3.7, and when $\alpha =0.1,\beta =0.7$, the numerical result of the test is close to 3.9. It may be noted from Table 1, Table 2 and Table 3 that the high-order numerical scheme is convergent and has good approximation.

Example 2.

Consider the following two-dimensional nonlinear fractional Volterra integral equations:

$\begin{array}{cc}\hfill & u(x,y)=f(x,y)+{\int }_0^x{\int }_0^y\frac{(xy+s+t)u^2(s,t)}{{(x-s)}^{\alpha }{(y-t)}^{\beta }}dtds,(x,y)\in [0,1]\times [0,1],\hfill \end{array}$

where

$\begin{array}{cc}\hfill & f(x,y)=x^3{y}^2-17280x^{8-\alpha }{y}^{6-\beta }{\displaystyle \prod _{i=1}^7}\frac{1}{(i-\alpha )}{\displaystyle \prod _{j=1}^5}\frac{1}{(j-\beta )}\hfill \\ \hfill & -120960x^{8-\alpha }{y}^{5-\beta }{\displaystyle \prod _{i=1}^8}\frac{1}{(i-\alpha )}{\displaystyle \prod _{j=1}^5}\frac{1}{(j-\beta )}-86400x^{7-\alpha }{y}^{6-\beta }{\displaystyle \prod _{i=1}^7}\frac{1}{(i-\alpha )}{\displaystyle \prod _{j=1}^6}\frac{1}{(j-\beta )},\hfill \end{array}$

and the exact solution is $u(x,y)=x^3{y}^2.$

In this example, the meanings described in Table 4 and Table 5 are similar to those in Table 1, Table 2 and Table 3 in Example 1. In Table 4, we can easily see that when $\alpha =0.3,\beta =0.6$ and $\alpha =0.1,\beta =0.7$, the numerical results of the test are very close to the theoretical values of 3.4 and 3.3, respectively. In Table 5, when $\alpha =0.3,\beta =0.6$, the numerical result of the test is close to 3.7, and when $\alpha =0.1,\beta =0.7$, the numerical result of the test is close to 3.9. This shows very good consistency with the theoretical prediction.

Next, we plot the error distribution. The error distribution of $M=2N,N=128$ is shown in Figure 1, where the error distribution of $\alpha =0.3,\beta =0.6$ is on the left and the error distribution of $\alpha =0.1,\beta =0.7$ is on the right. From Figure 1, we find that the errors can be as small as ${10}^{-7}$ and ${10}^{-8}$.

Example 3.

Assume the below two-dimensional fractional Volterra integral equation [24]:

$\begin{array}{cc}\hfill u(x,y)=& x^2(y^2-y)-0.0828x^{2}y(3y-4)\hfill \\ \hfill & +\frac{1}{\Gamma (\frac{2}{3})\Gamma (\frac{2}{3})}{\int }_0^x{\int }_0^y{(x-s)}^{\frac{-1}{3}}{(y-t)}^{\frac{-1}{3}}{(xy)}^{-\frac{2}{3}}u(s,t)dtds,\hfill \end{array}$

where the exact solution is $u(x,y)=x^2(y^2-y)$.

In this example, we take $M=N=32$ and compare our scheme with the numerical method in [24]. The numerical solutions of different points $(x,y)\in [0,1]\times [0,1]$ are given in Table 6, where “Haar solution” denotes the numerical solutions in Table 2 of [24], and “Numerical solution” denotes the numerical solutions obtained by our method. Since we are using a quadratic polynomial, while a zero polynomial is used in [24], combined with the data in Table 6, we can clearly see that the error produced by using our proposed method is much smaller.

Example 4.

We consider the following two-dimensional nonlinear fractional Volterra integral equation [25]:

$\begin{array}{cc}\hfill u(x,y)& =\sqrt{y}(\frac{-1}{180}{x}^3{y}^{\frac{7}{2}}+\sqrt{\frac{x}{3}})\hfill \\ \hfill & +\frac{1}{\Gamma (\frac{3}{2})\Gamma (\frac{5}{2})}{\int }_0^x{\int }_0^y{(x-s)}^{\frac{1}{2}}{(y-t)}^{\frac{3}{2}}\sqrt{xyt}{u}^2(s,t)dtds,\hfill \end{array}$

with the exact solution $u(x,y)=\frac{\sqrt{3xy}}{3}$.

In this example, we take $M=N=8,32$ and give the exact and approximate solutions at the point $(x,y)\in [0,1]\times [0,1]$ in Table 7, where "Method solution" represents the approximate solution in Table 1 of [25] and "Numerical solution" represents the numerical solution using our proposed method. The maximum error is also given in Table 7. According to the calculation results shown in Table 7, we can see that the approximate solution obtained using our proposed method is closer to the exact solution at the corresponding point, and the maximum error is smaller.

6. Concluding Remarks

Following the modified block-by-block method [5], we proposed an efficient and high-order numerical scheme (32) to approximate the solution of two-dimensional nonlinear fractional Volterra integral equations with singular kernels. The main idea of the high-order numerical scheme was to discretize its domain into a number of subdomains and then use biquadratic interpolation on each subdomain. In this high-order numerical scheme, only the two boundary layers were coupled, and the rest of the blocks were explicit in the scheme, which made our calculations more convenient. The scheme has uniform accuracy, and the optimal convergence order was $O(h_x^{4-\alpha }+h_y^{4-\beta })$. For the high-order numerical scheme, we conducted a detailed error analysis and verified the correctness of the theoretical analysis through numerical examples. The advantages of the numerical scheme (32) are high accuracy, easy calculation and the fact that there is no need for coupling solutions. The disadvantage of the numerical scheme is that the computation time is long, increasing as the number of partitions increases. In the future, according to reference [26,27], we intend to use a fast algorithm to realize the results. Based on the idea of [28], we will use the high-order numerical scheme to solve Fredholm–Hammerstein integral equations of the second kind. Based on the ideas of [29,30,31,32], we expect that the constructed efficient higher-order scheme can be applied to engineering and practical problems.

Author Contributions

Funding acquisition, Z.-Q.W. and J.-Y.C.; investigation, Z.-Q.W. and Q.L.; methodology, Z.-Q.W. and J.-Y.C.; project administration, Z.-Q.W. and J.-Y.C.; software, Q.L.; supervision, Z.-Q.W. and J.-Y.C.; visualization, Q.L.; writing—original draft, Z.-Q.W. and Q.L.; writing—review and editing, Z.-Q.W., Q.L. and J.-Y.C.; All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All the data were computed using our high-order numerical scheme.

Acknowledgments

The authors thank the anonymous referees for their valuable suggestions, which improved the quality of this work significantly.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses or interpretation of data, in the writing of the manuscript or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figure and Tables

Figure 1. Error distribution of [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.] (right) for [Forumla omitted. See PDF.].

Table 1

Maximum errors, decay rate and CPU time with $M=2N$.

N	$\mathit{\alpha }=0.3,\mathit{\beta }=0.6$	Order	CPU Time
8	1.9973659945$\times {10}^{-3}$	-	0.1438033000 s
16	2.0439893391$\times {10}^{-4}$	3.2886391376	1.7739478000 s
32	2.0217387802$\times {10}^{-5}$	3.3377191614	27.3948994000 s
64	1.9623458838$\times {10}^{-6}$	3.3649453461	7.1728836283 min
128	1.8839724669$\times {10}^{-7}$	3.3807295678	1.8958612181 h

Table 2

Maximum errors, decay rate and CPU time with $M=2N$.

N	$\mathit{\alpha }=0.1,\mathit{\beta }=0.7$	Order	CPU Time
8	1.0362949128$\times {10}^{-3}$	-	0.2546744000 s
16	1.1195268795$\times {10}^{-4}$	3.2104735569	1.8277023000 s
32	1.1762915084$\times {10}^{-5}$	3.2505716291	27.406621300 s
64	1.2176430484$\times {10}^{-6}$	3.2720824624	7.1583191750 min
128	1.2502313451$\times {10}^{-7}$	3.2838242818	1.8919108256 h

Table 3

Maximum errors and decay rate with $N=\left[M^{(4-\alpha )/(4-\beta )}\right]$.

M	$\mathit{\alpha }=0.3,\mathit{\beta }=0.6$	Order	$\mathit{\alpha }=0.1,\mathit{\beta }=0.7$	Order
$\phantom{00}8$	2.7093547741$\times {10}^{-3}$	−	8.5837427821$\times {10}^{-4}$	−
$\phantom{0}16$	2.4017036964$\times {10}^{-4}$	3.4958192397	6.7887046124$\times {10}^{-5}$	3.6603986275
$\phantom{0}32$	2.0350262447$\times {10}^{-5}$	3.5609388681	5.0238828559$\times {10}^{-6}$	3.7562615829
$\phantom{0}64$	1.6641178999$\times {10}^{-6}$	3.6122178456	3.6303832318$\times {10}^{-7}$	3.7906090690
128	1.3366509233$\times {10}^{-7}$	3.6380630008	2.5643266843$\times {10}^{-8}$	3.8234698786

Table 4

Maximum errors and decay rate with $M=2N$.

N	$\mathit{\alpha }=0.3,\mathit{\beta }=0.6$	Order	$\mathit{\alpha }=0.1,\mathit{\beta }=0.7$	Order
$\phantom{00}8$	1.4497410528$\times {10}^{-3}$	−	6.7268209314$\times {10}^{-4}$	−
$\phantom{0}16$	1.4126072922$\times {10}^{-4}$	3.3593628803	6.9945531354$\times {10}^{-5}$	3.2656105874
$\phantom{0}32$	1.3572217464$\times {10}^{-5}$	3.3796320927	7.1896243538$\times {10}^{-6}$	3.2822435908
$\phantom{0}64$	1.2887160530$\times {10}^{-6}$	3.3966501211	7.3370971254$\times {10}^{-7}$	3.2926351059
128	1.2153830342$\times {10}^{-7}$	3.4064514616	7.4615217116$\times {10}^{-8}$	3.2976675925

Table 5

Maximum errors and decay rate with $N=\left[M^{(4-\alpha )/(4-\beta )}\right]$.

M	$\mathit{\alpha }=0.3,\mathit{\beta }=0.6$	Order	$\mathit{\alpha }=0.1,\mathit{\beta }=0.7$	Order
$\phantom{00}8$	4.1299227631$\times {10}^{-3}$	−	1.2444932111$\times {10}^{-3}$	−
$\phantom{0}16$	3.5190750852$\times {10}^{-4}$	3.5528465991	9.1544255179$\times {10}^{-5}$	3.7649451968
$\phantom{0}32$	2.9677317197$\times {10}^{-5}$	3.5677637121	6.6471665752$\times {10}^{-6}$	3.7836579403
$\phantom{0}64$	2.4390547260$\times {10}^{-6}$	3.6049666456	4.7740052334$\times {10}^{-7}$	3.7994674546
128	1.9706716024$\times {10}^{-7}$	3.6295628422	3.3836403013$\times {10}^{-8}$	3.8185520322

Table 6

The numerical results in Example 3.

$({\mathit{x}}_{\mathit{i}},{\mathit{y}}_{\mathit{i}})$	Exact Solution	Haar Solution in [24]	Numerical Solution
$(0,0)$	0	0	0
$(0.1,0.1)$	−0.000900	−0.000931	−0.000899999999999
$(0.2,0.2)$	−0.006400	−0.006458	−0.006400000000000
$(0.3,0.3)$	−0.018900	−0.018932	−0.018900000000000
$(0.4,0.4)$	−0.038400	−0.038412	−0.038400000000000
$(0.5,0.5)$	−0.062500	−0.062504	−0.062499999999999
$(0.6,0.6)$	−0.086400	−0.086475	−0.086399999999997
$(0.7,0.7)$	−0.102900	−0.102941	−0.102899999999998
$(0.8,0.8)$	−0.102400	−0.102482	−0.102399999999991
$(0.9,0.9)$	−0.072900	−0.072970	−0.072899999999974
$(1,1)$	0	0	0.000000000000137

Table 7

The numerical results for Example 4.

$\mathit{x}=\mathit{y}$	Exact Solution	Method Solution in [25]		Numerical Solution
		$\mathbf{h}=\mathbf{1}/\mathbf{16}$	$\mathbf{h}=\mathbf{1}/\mathbf{64}$	$\mathbf{h}=\mathbf{1}/\mathbf{16}$	$\mathbf{h}=\mathbf{1}/\mathbf{64}$
0	0	0.018452	0.002162	0	0
0.1	0.057735	0.031135	0.049679	0.062761	0.057660
0.2	0.11547	0.13261	0.117302	0.114850	0.115427
0.3	0.173205	0.147605	0.168105	0.171241	0.173187
0.4	0.23094	0.246768	0.232442	0.230698	0.230858
0.5	0.288675	0.262075	0.285535	0.288662	0.288662
0.6	0.34641	0.360925	0.347581	0.346202	0.346324
0.7	0.404145	0.378535	0.398155	0.404047	0.404057
0.8	0.46188	0.475083	0.462721	0.461606	0.461757
0.9	0.519615	0.501015	0.513525	0.519441	0.519463
0.99	0.571577	0.583533	0.572104	0.570976	0.571549
Max error	0	2.96$\times {10}^{-2}$	8.09$\times {10}^{-3}$	5.02$\times {10}^{-3}$	1.51$\times {10}^{-4}$

References

1. Kimeu, J. Fractional Calculus: Definitions and Applications. Master’s Thesis; Western Kentucky University: Bowling Green, Kentucky, USA, 2009.

2. Micula, S. A numerical method for weakly singular nonlinear Volterra integral equations of the second kind. Symmetry; 2020; 12, 1862.

3. Saberi-Nadjafi, J.; Heidari, M. A generalized block-by-block method for solving linear Volterra integral equations. Appl. Math. Comput.; 2007; 188, pp. 1969-1974. [DOI: https://dx.doi.org/10.1016/j.amc.2006.11.101]

4. Katani, R.; Shahmorad, S. Block by block method for the systems of nonlinear Volterra integral equations. Appl. Math. Model.; 2010; 34, pp. 400-406. [DOI: https://dx.doi.org/10.1016/j.apm.2009.04.013]

5. Cao, J.; Xu, C. A high order schema for the numerical solution of the fractional ordinary differential equations. J. Comput. Phys.; 2013; 238, pp. 154-168. [DOI: https://dx.doi.org/10.1016/j.jcp.2012.12.013]

6. Assari, P.; Dehghan, M. The approximate solution of nonlinear Volterra integral equations of the second kind using radial basis functions. Appl. Numer. Math.; 2018; 131, pp. 140-157. [DOI: https://dx.doi.org/10.1016/j.apnum.2018.05.001]

7. Al-Bugami, A.M.; Al-Juaid, J.G. Runge-Kutta method and block by block method to solve nonlinear Fredholm-Volterra integral equation with continuous Kernel. J. Appl. Math. Phys.; 2020; 8, pp. 2043-2054. [DOI: https://dx.doi.org/10.4236/jamp.2020.89152]

8. Jaabar, S.M.; Hussain, A.H. Solving Volterra integral equation by using a new transformation. J. Interdiscip. Math.; 2021; 24, pp. 735-741. [DOI: https://dx.doi.org/10.1080/09720502.2020.1861787]

9. Usta, F. Bernstein approximation technique for numerical solution of Volterra integral equations of the third kind. Comput. Appl. Math.; 2021; 40, pp. 161-172.

10. Ma, Z.; Alikhanov, A.A.; Huang, C.; Zhang, G. A multi-domain spectral collocation method for Volterra integral equations with a weakly singular kernel. Appl. Numer. Math.; 2021; 167, pp. 218-236. [DOI: https://dx.doi.org/10.1016/j.apnum.2021.05.006]

11. Mckee, S.; Tang, T.; Diogo, T. An Euler-type method for two-dimensional Volterra integral equations of the first kind. IMA J. Numer. Anal.; 2000; 20, pp. 423-440. [DOI: https://dx.doi.org/10.1093/imanum/20.3.423]

12. Wang, Y.; Huang, J.; Wen, X. Two-dimensional Euler polynomials solutions of two-dimensional Volterra integral equations of fractional order. Appl. Numer. Math.; 2021; 163, pp. 77-95. [DOI: https://dx.doi.org/10.1016/j.apnum.2021.01.007]

13. Abdi, A.; Conte, D. Implementation of general linear methods for Volterra integral equations. J. Comput. Appl. Math.; 2021; 386, pp. 1-19. [DOI: https://dx.doi.org/10.1016/j.cam.2020.113261]

14. Dehbozorgi, R.; Nedaiasl, K. Numerical solution of nonlinear weakly singular Volterra integral equations of the first kind: An hp-version collocation approach. Appl. Numer. Math.; 2021; 161, pp. 111-136. [DOI: https://dx.doi.org/10.1016/j.apnum.2020.10.030]

15. Ahsan, S.; Nawaz, R.; Akbar, M. Nisar. K.S.; Baleanu, D. Approximate solutions of nonlinear two-dimensional Volterra integral equations. Math. Methods Appl. Sci.; 2021; 44, pp. 5548-5559. [DOI: https://dx.doi.org/10.1002/mma.7128]

16. Kant, K.; Nelakanti, G. Approximation methods for second kind weakly singular Volterra integral equations. J. Comput. Appl. Math.; 2020; 368, 112531.

17. Pan, Y.; Huang, J.; Ma, Y. Bernstein series solutions of multidimensional linear and nonlinear Volterra integral equations with fractional order weakly singular kernels. Appl. Math. Comput.; 2019; 347, pp. 149-161. [DOI: https://dx.doi.org/10.1016/j.amc.2018.10.022]

18. Allaei, S.S.; 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.

19. Li, X.; Tang, T.; Xu, C. Numerical solutions for weakly singular Volterra integral equations using Chebyshev and Legendre Pseudo-Spectral Galerkin methods. J. Sci. Comput.; 2016; 67, pp. 43-64. [DOI: https://dx.doi.org/10.1007/s10915-015-0069-5]

20. Abdelrahman, M.; Hassan, S.; Alomair, R.; Alsaleh, D. Fundamental solutions for the conformable time fractional Phi-4 and space-time fractional simplified MCH equations. AIMS Math.; 2021; 6, pp. 6555-6568. [DOI: https://dx.doi.org/10.3934/math.2021386]

21. Abdelrahman, M.; Sohaly, M.; Alharbi, Y. Fundamental stochastic solutions for the conformable fractional NLSE with spatiotemporal dispersion via exponential distribution. Phys. Scr.; 2021; 96, 125223. [DOI: https://dx.doi.org/10.1088/1402-4896/ac119c]

22. Wang, K.; Wang, G.; Zhu, H. A new perspective on the study of the fractal coupled Boussinesq-Burger equation in shallow water. Fractals; 2021; 29, 2150122. [DOI: https://dx.doi.org/10.1142/S0218348X2150122X]

23. Dixon, J.; Mckee, S. Weakly singular discrete Gronwall inequalities. Z. Angew. Math. Mech.; 1978; 66, pp. 535-544. [DOI: https://dx.doi.org/10.1002/zamm.19860661107]

24. Abdollahi, Z.; Moghadam, M.M.; Saeedi, H.; Ebadi, M. A computational approach for solving fractional Volterra integral equations based on two-dimensional Haar wavelet method. Int. J. Comput. Math.; 2022; 99, pp. 1488-1504.

25. Najafalizadeh, S.; Ezzati, R. Numerical methods for solving two-dimensional nonlinear integral equations of fractional order by using two-dimensional block pulse operational matrix. Appl. Math. Comput.; 2016; 280, pp. 46-56. [DOI: https://dx.doi.org/10.1016/j.amc.2015.12.042]

26. Micula, S. A fast converging iterative method for Volterra integral equations of the second kind with delayed arguments. Fixed Point Theory.; 2015; 16, pp. 371-380.

27. Gu, X.; Wu, S. A parallel-in-time iterative algorithm for Volterra partial integro-differential problems with weakly singular kernel. J. Comput. Phys.; 2020; 417, 109576. [DOI: https://dx.doi.org/10.1016/j.jcp.2020.109576]

28. Micula, S.; Cattani, C. On a numerical method based on wavelets for Fredholm-Hammerstein integral equations of the second kind. Math. Methods Appl. Sci.; 2018; 41, pp. 9103-9115.

29. Cao, J.; Cai, Z. Numerical analysis of a high-order scheme for nonlinear fractional differential equations with uniform accuracy. Numer. Math. Theor. Meth. Appl.; 2021; 14, pp. 71-112. [DOI: https://dx.doi.org/10.4208/nmtma.OA-2020-0039]

30. Wang, Z.; Cui, J. Second-order two-scale method for bending behavior analysis of composite plate with 3-D periodic configuration and its approximation. Sci. China Math.; 2014; 57, pp. 1713-1732. [DOI: https://dx.doi.org/10.1007/s11425-014-4831-1]

31. Huang, Y.; Gu, X.; Gong, Y.; Li, H.; Zhao, Y.; Carpentieri, B. A fast preconditioned semi-implicit difference scheme for strongly nonlinear space-fractional diffusion equations. Fractal Fract.; 2021; 5, 230. [DOI: https://dx.doi.org/10.3390/fractalfract5040230]

32. Gu, X.; Sun, H.; Zhao, Y.; Zheng, X. An implicit difference scheme for time-fractional diffusion equations with a time-invariant type variable order. Appl. Math. Lett.; 2021; 120, 107270. [DOI: https://dx.doi.org/10.1016/j.aml.2021.107270]