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

: In this paper, based on the modiﬁed 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 4 − α x + h 4 − β y ) for 0 < α , β < 1. Finally, experiments with four numerical examples are shown, to support the theoretical ﬁndings and to illustrate the efﬁciency of our proposed method.


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: u(x, y) = f (x, y) + x a y c κ(x, y, s, t, u(s, t)) (x − s) α (y − t) β dtds, (x, y) ∈ D, 0 < α, β < 1, where f (x, y) and κ(x, y, s, t, u(s, t)) are given continuous functions defined on D = [a, b] × [c, d], Ω = D × D × 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 κ(x, y, s, t, u(s, t) satisfies the following condition |κ(x, y, s, t, u 1 (s, t)) − κ(x, y, s, t, u 2 (s, t))| ≤ L|u 1 (s, t) − u 2 (s, t)|, L > 0.
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 β (z − x) −α g(y(x)), where α ∈ (0, 1), β > 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 4−α x + h 4−β y ) for 0 < α, β < 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.

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 × 2N equal subdomains with size h x = b−a 2M , h y = d−c 2N , where M and N are positive integers. We denote x i = a + ih x ,y j = c + jh y , i = 0, 1, · · · , 2M; j = 0, 1, · · · , 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 κ i,j (s, t, u(s, t)) = κ(x i , y j , s, t, u(s, t)), f i,j = f (x i , y j ).

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 whereū 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ū i,j at the point (x 2m+1 , y 2n+1 ) as follows ).
For ease of notation, let ∂ 3 κ Then, the estimation of r i,j is as follows.
Combining (39) and (40) yields We use the same technique for r (2) 2m+1,2n+1 , and it is easy to see that For the first item on the right side of (42), we can obtain Next, we make a detailed estimate of R 2 with the following form wheret j =y 2j . Based on the definition of ϕ k,i (s), s ∈ (x i , x i+2 ), k = 0, 1, 2; i ∈ N, we know that |ϕ k,i (s)| ≤ 1, and therefore | 2 ∑ k=0 ϕ k,i (s)| ≤ 3, and we have We make the following detailed estimates of A 11 and A 12 , which can be directly obtained.
For R

, we have
Taking the first term and the second term on the right side of the above formula as B 1 and B 2 , we have For B 11 and B 12 , we use the same method as (46) and (47), obtaining Hence, for B 1 , we can directly obtain For B 2 , we can directly obtain the result Next, we estimate R 3 as follows Therefore, we obtain that Finally, we make a detailed estimation of r (4) 2m+1,2n+1 , which has the form |r (4) x 2i+1 x 2i−1 y 2j+1

and R
(2) 4 one by one. For R (1)

, we have
For P 1 , the same processing method can be used as for B 1 , to obtain For P 2 , we can obtain Finally, we estimate R 4 as follows wheret 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 Therefore, we obtain Bringing (41) and (51)-(53) into (37) gives where C depends only on G 1 , G 2 , α and β. 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 The proof is then completed.
Based on the idea of [5], we can prove the following lemma.

Theorem 1. Let u be the exact solution of Equation
then the following error estimates hold where C depends only on L, κ, b, d, α and β.
Secondly, e i,j , i ≥ 3; j = 1, 2, satisfies whereB m i and B m i 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 ||ê i = max{|e i,1 |, |e i,2 |, i = 0, 1, · · · , 2M}, ||r i = max{|r i,1 |, |r i,2 |, i = 0, 1, · · · , 2M}, and then we can convert them into two groups of similar formulas, one of which is For ê 2m+1 , we have After applying the Gronwall's inequality [23] to the above equation, we obtain Using the same method for ê 2m+2 , we directly obtain The process of calculating e i,j , i = 1, 2; j ≥ 3, is similar to that for e i,j , i ≥ 3; j = 1, 2, and therefore it is omitted.
Taking the above result, together with Lemma 1, we have Next, e i,j , i, j ≥ 3 satisfy whereB m i , B m i ,D n j and D n j are defined in (55) and satisfy Lemma 2. We estimate e 2m+1,2n+1 as shown below.
which can be rearranged into the preferred form As the above formula is applicable to all n = 1, 2, . . . , N − 1, we can obtain Therefore, Using the discrete Gronwall inequality [23], we obtain Combining the above estimations with Lemma 1 yields In the same way, we can obtain the results We combine (59), (60), (63) and (64), then the proof is completed.

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 .
The corresponding exact solution is u(x, y) = x 4 y 4 .
In our test, we chose α = 0.3, β = 0.6 and α = 0.1, β = 0.7, respectively. The step size in the x-direction was h x = 1 2M and the step size in the y-direction was h y = 1 2N . The definition of error was as follows 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 Tables 1 and 2 Tables 1 and 2, we take α = 0.3, β = 0.6 or α = 0.1, β = 0.7 and M = 2N. In this case, our scheme's theoretical order is 4 − β. In Tables 1 and 2, we show the corresponding convergence order and CPU time when N takes a series of values and α, β are given specific values. In Table 1, when α = 0.3, β = 0.6, the tested order is close to the theoretical order of 3.4. In Table 2, when α = 0.1, β = 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 Tables 1 and 2, we find that the CPU time quickly increases from 0.14 s to about 2 h as N increases.
The test results are given in Table 3. When α = 0.3, β = 0.6, the numerical result of the test is close to 3.7, and when α = 0.1, β = 0.7, the numerical result of the test is close to 3.9. It may be noted from Tables 1-3 that the high-order numerical scheme is convergent and has good approximation.
, and the exact solution is u(x, y) = x 3 y 2 .
In this example, the meanings described in Tables 4 and 5 are similar to those in Tables 1-3 in Example 1. In Table 4, we can easily see that when α = 0.3, β = 0.6 and α = 0.1, β = 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 α = 0.3, β = 0.6, the numerical result of the test is close to 3.7, and when α = 0.1, β = 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 α = 0.3, β = 0.6 is on the left and the error distribution of α = 0.1, β = 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]: u(x, y) =x 2 (y 2 − y) − 0.0828x 2 y(3y − 4) 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) ∈ [0, 1] × [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.    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.

Concluding Remarks
Following the modified block-by-block method [5], we proposed an efficient and highorder 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 4−α x + h 4−β y ). 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.