A higher-order uniform accuracy scheme for nonlinear $\psi$-Volterra integral equations in two dimension with weakly singular kernel

1.   Introduction

Volterra integral equations (VIEs) are important in many fields and they have been studied extensively. The purpose of this paper is to investigate the higher-order uniform accuracy numerical solution of nonlinear $\psi$-VIEs with weakly singular kernel in two dimension,

where $(s, t)\in H, 0 < \sigma_{1}, \sigma_{2} < 1$, $\nu(s, t)$ is an unknown function defined in $H = [a, b]\times [c, d]$. $g, \omega$ are known functions and defined in $H$ and $\Omega = H\times H\times R$, respectively. $\psi$ is a monotonically increasing function and $\psi' > 0$.

Various numerical methods for solving the spectral form of the $\psi$-VIEs (1.1) have already been constructed by the direct quadrature methods ^[1,2], spectral methods ^[3,4], difference methods ^[5], operational matrix methods ^[6,7], block-by-block method^[8], and Taylor polynomials^[9]. More numerical methods are available for readers to read ^{[10,11,12,13,14]}. For the numerical solution for the general $\psi$-VIEs (1.1), many researchers have studied this topic and achieved some research results. For example, in ^[16], the Chebyshev wavelets were employed to establish a computational procedure for $\psi$-VIEs. In ^[17], they established a Maxwell model with the $\psi$-Caputo fractional derivative. In ^[18], they presented a novel pharmacokinetic/pharmacodynamic model with $\psi$-Caputo fractional derivatives for the induction phase of anesthesia. In ^[19], they presented a generalization of the $\psi$-Hilfer fractional derivative. In ^[20], they investigated the existence, uniqueness, and Ulam-Hyers stabilities of mild solutions of $\psi$-fractional differential equations. They investigated the multiple $\psi$-type stability of fractional-order quaternion-valued neural networks in ^[21]. In ^[22], they studied a new class of impulsive boundary value problems with the generalized $\psi$-Caputo fractional derivatives. In ^[23], the stability of the autonomous linear $\alpha \in (0, 1)$ order $\psi$-Caputo fractional differential system was investigated. In ^[24], they presented a numerical method for solving the $\psi$-Caputo fractional differential equations. In ^[25], they solved a nonsingularity kernel VIEs by the fourth order Lagrange basis function.

Therefore, the main aim of this work is to propose an efficient high order time uniform numerical scheme for nonlinear $\psi$-VIEs with weakly singular kernel in two dimension (1.1) by using biquadratic interpolation in each sub-domain with the convergence order $O(h_{s}^{3+\sigma_{1} }+h_{t}^{3+\sigma_{2}})$ with $0 < \sigma_{1}, \sigma_{2} < 1$. The modified block-by-block method is introduced to discrete nonlinear $\psi$-VIEs with weakly singular kernel in two dimension. The proposed scheme can avoid convergence accuracy degeneracy near the two boundary layers by coupled calculating the numerical solutions, which achieves the sharp convergence order without the fine mesh near the boundary layer. The proposed scheme is efficient and does not require coupled solutions in the other sub-domains. The numerical scheme will be established in this article to provide a paradigm for establishing a high-order scheme of the nonlinear $\psi$-VIEs with weakly singular kernel in two dimension and analyzing its convergence and stability of a high-order numerical scheme. It can also provide readers with a feasible method for constructing high order time numerical schemes and their theoretical analysis for similar nonlinear $\psi$-VIEs with weakly singular kernel in two dimension.

The article follows this structure: A higher-order uniform accuracy numerical scheme is introduced in Section 2. Section 3 estimates the truncation error of the higher-order uniform accuracy numerical scheme. Section 4 analyzes the convergence of the higher-order uniform accuracy numerical scheme. In Section 5, we demonstrate the efficiency of the higher-order uniform accuracy numerical scheme and supports our theoretical findings through two numerical experiments. Finally, some conclusions are drawn in Section 6.

2.   Higher-order numerical approach for nonlinear $\psi$-fractional VIEs in two dimension

In this section, we will consider the following nonlinear $\psi$-fractional VIEs in two dimension (1.1) based on the idea of ^[26,27]. For $\omega(s, t, \tau, \eta, \nu(\tau, \eta))$, we assume that it satisfies the following Lipschitz continuity for the fifth variable

Next, we will approximate the nonlinear $\psi$-fractional VIEs in two dimension of (1.1) by using biquadratic Lagrange interpolation. Denote positive integers $L, K$, and $h_{s} = \frac{b-a}{2L}$, $h_{t} = \frac{d-c}{2K}$. Let $s_{i} = a+ih_{s}(0\leq i \leq 2L), t_{j} = c+jh_{t} (0\leq j \leq 2K)$, $\nu_{i, j}$ be the numerical solution at $(s_{i}, t_{j})$, and $\omega^{i, j}(\tau, \eta, \nu(\tau, \eta)) = \omega(s_{i}, t_{j}, \tau, \eta, \nu(\tau, \eta))$, $g_{i, j} = g(s_{i}, t_{j})$. By using simple calculations, one can obtain that $\nu_{i, 0} = g(s_{i}, 0)$, $\nu_{0, j} = g(0, t_{j})$.

Let $\hat{f}_{p, i}(s), p = 0, 1, 2, 0\leq i \leq 2L-2$ and $f_{q, j}(t), q = 0, 1, 2, 0\leq j \leq 2K-2$ be the quadratic Lagrange basis functions at points ${{s}_{i}}, {{s}_{i+1}}, {{s}_{i+2}}$ and ${{t}_{j}}, {{t}_{j+1}}, {{t}_{j+2}}$, respectively,

Let $\nu_{i, j}$, $\nu_{i, 2k+1}$, $\nu_{i, 2k+2}$, $\nu_{2l+1, j}$, and $\nu_{2l+2, j}$ be known, where $0\le i\le 2l, 1\le l\le L-1, 0\le j\le2k, 1\le k\le K-1, i, j, k, l$ are integers. We will construct an approximate scheme for $\nu_{2l+1, 2k+1}$, $\nu_{2l+1, 2k+2}$, $\nu_{2l+2, 2k+1}$, and $\nu_{2l+2, 2k+2}$.

For $\nu(s_{2l+1}, t_{2k+1})$, we have

For $D_{1}$, one can obtain that

where $A_{2l+1}^{p, 0}$, $B_{2k+1}^{q, 0}$ are defined by

Similarly, for ${{D}_{2}}, {{D}_{3}}$, and ${{D}_{4}}$, we have

where $A_{2l+1}^{p, i}$ and $B_{2k+1}^{q, j}$ are defined as follows:

and $A_{2l+1}^{p, 0}$, $B_{2k+1}^{q, 0}$ are defined by (2.4) and (2.5), respectively.

Bringing (2.3) and (2.6)–(2.8) into (2.2), we have

For $\nu(s_{2l+2}, t_{2k+1})$, $\nu(s_{2l+1}, t_{2k+2})$, and $\nu(s_{2l+2}, t_{2k+2})$, we have

where $A_{2l+2}^{p, i}$ and $B_{2k+2}^{q, j}$ are defined as follows:

and $A_{2l+1}^{p, 0}$, $A_{2l+1}^{p, i}$, $B_{2k+1}^{q, 0}$, and $B_{2k+1}^{q, j}$ are defined in (2.4), (2.9), (2.5), and (2.10).

Similar to the same line of (2.11), we can obtain the other point's numerical scheme. A detailed derivation is given in Appendix A.

Combining with (A.1), (A.4)–(A.6), (A.9), (A.10)–(A.12), (A.13)–(A.16), (2.11)–(2.14), we obtain a high-order numerical scheme for (1.1) as $i_1, j_1 = 1, 2,$

3.   Truncation error

Let $E_{i, j}$ be the truncation error of scheme (2.17) at point $(s_{i}, t_{j})$ as follows:

where $\bar{\nu}_{i, j}$ is the numerical solution for $\nu(s_{i}, t_{j})$. For instance, at the point $(s_{2l+1}, t_{2k+1})$, $\bar{\nu}_{i, j}$ is defined as follows:

In order to analyze error estimation for the scheme (2.17), we will introduce the following lemmas.

Lemma 3.1. For all $\psi \in C^{1}(I)$, an increasing function such that $\psi'(s) > 0$, we have

where

Proof. Using the integral mean value theorem and Lagrange mean value theorem, one can get

where $\varrho_{1}\in(\psi(a), \psi(s_{1}))$, $\theta _{2l+1}\in(s_{1}, s_{2l+1})$, $\theta _{1}\in(a, s_{1})$, and $r$ is shown in (3.3).

Similarly, $\int_{c}^{t_{1}}(\psi (t_{2k+1})-\psi (\eta))^{\sigma_{2}-1} d\psi(\eta)\le2^{1-\sigma_{2}}\bar{r}^{\sigma_{2}} (2k+1)^{\sigma_{2}-1}h _{t}^{\sigma_{2}},$ where $\bar{r}$ is shown in (3.3). Therefore, Lemma 3.1 has been proven. □

Lemma 3.2. Let $E_{i, j}$ be defined by (3.1), $\omega (\cdot, \cdot, \cdot, \cdot, \nu(\cdot, \cdot))\in C^{4}([a, b]\times[c, d])$, and $\nu(\cdot, \cdot)\in C^{4}([a, b]\times[c, d])$, then

where $C$ is a constant independent of $h_{s}, h_{t}$.

Proof. For a detailed proof, see Appendix B. □

4.   Convergence analysis

Lemma 4.1. The coefficients are estimated as follows:

where $C$ is independent of $h _{s}$ and $h _{t}$. In equation (A.2), (A.7), (2.4), (2.9), (2.15) and (A.3), (A.8), (2.5), (2.10), (2.16), we have the definition of $A_{1}^{p, 0}$, $A_{2}^{p, 0}$, $A_{2l+1}^{p, 0}$, $A_{2l+1}^{p, i}$, $A_{2l+2}^{p, i}$ and $B_{1}^{q, 0}$, $B_{2}^{q, 0}$, $B_{2k+1}^{q, 0}$, $B_{2k+1}^{q, j}$, $B_{2k+2}^{q, j}$, where $p, q = 0, 1, 2$, $i = 0, 1, \dots, l-1$, $j = 0, 1, \dots, k-1$, $l = 1, 2, \dots, L-1$, and $k = 1, 2, \dots, K-1$.

Proof. By employing the Lagrange mean value theorem, we have

where $r$ is defined by (3.3), $\xi _{1}\in (a, s_{1})$, and

where $\xi' _{0}\in (a, s_{2})$.

In the same way, we can get

For the estimation of $|A_{2l+1}^{0, 0}|$, we can use the integral mean value theorem and Lagrange mean value theorem,

Similarly, we have

We can use the integral mean value theorem and Lagrange mean value theorem to get

where $\xi _{2i-1}\in(s_{2i-1}, s_{2i+1})$.

In the same way, we have

When $i = l$,

Similarly, we have $|A_{2l+2}^{p, l}| \le C 2^{\sigma_1} h_s^{\sigma_1}$. Using the same approach, we can estimate $|B_{1}^{q, 0}|, |B_{2}^{q, 0}|, |B_{2k+1}^{q, 0}|, |B_{2k+1}^{q, j}|, |B_{2k+2}^{q, j}|, |B_{2k+1}^{q, k}|, |B_{2k+2}^{q, k}|$ to be

The proof of Lemma 4.1 is then completed. □

Theorem 4.1. Let $\nu(s_{i}, t_{j})$ be the exact solution of formula (1.1) and $\nu_{i, j}$ be the numerical solution. If $w$ satisfies condition (2.1), $\omega (\cdot, \cdot, \cdot, \cdot, \nu(\cdot, \cdot))\in C^{4}([a, b]\times[c, d])$, and $\nu(\cdot, \cdot)\in C^{4}([a, b]\times[c, d])$, the step size satisfies the following formula:

We have

where $C$ is independent of $h_{s}$ and $h _{t}$.

Proof. $\forall i = 0, 1, \dots, 2L;j = 0, 1, \dots, 2K$; let $e_{i, j} = \nu(s_{i}, t_{j})-\nu_{i, j} = \nu(s_{i}, t_{j})-\bar{\nu}_{i, j}+\bar{\nu}_{i, j}-\nu_{i, j} = \bar{\nu}_{i, j}-\nu_{i, j}+E_{i, j}$. When $i, j = 1, 2$, we have $e_{i, j}$ as follows:

According to (2.1) and Lemma 4.1, there is

Combining the above four inequalities, we have

That is,

When $i\ge 3$, $j = 1, 2$, we have

where $A_{2l+1}^{p, 0}$, $A_{2l+1}^{p, i}$ and $B_{1}^{q, 0}$, $B_{2}^{q, 0}$ are defined in (2.4), (2.9) and (A.3), (A.8), and it satisfies with Lemma 4.1 that we have

where $k_{1} = 1, 2$.

For the sake of convenience, we take

Therefore, we get

For (4.4), we have

For inequality (4.6), use Gronwall inequality ^[28]. We have

Using the same method for (4.5), we can get

When $i = 1, 2, j \ge 3$ computing $e_{i, j}$ is analogous to computing $e_{i, j}$ for $i \ge 3, j = 1, 2$; hence, we omitted it. According to the results above, combined with Lemma 3.2, we have

Next, for $e_{i, j}$, $i, j\ge3$, we obtain

where $A_{2l+1}^{p, 0}, A_{2l+1}^{p, i}, A_{2l+2}^{p, i}$ and $B_{2k+1}^{q, 0}, B_{2k+1}^{q, j}, B_{2k+2}^{q, j}$ are defined in (2.4), (2.9), (2.15) and (2.5), (2.10), (2.16), and they are satisfied with Lemma 4.1. We have the following estimates for $e_{2l+1, 2k+1}$,

Therefore, we can rearrange into the preferred form:

If letting $\left \| e_{i} \right \| = \underset{\begin{smallmatrix} 0\le j \le 2K \\ \end{smallmatrix}}{\mathop{\max\limits}}\ |e_{i, j}|, \left \| E_{i} \right \| = \underset{\begin{smallmatrix} 0\le j \le 2K \\ \end{smallmatrix}}{\mathop{\max\limits}}\ |E_{i, j}|$, then we have

Since the above formula for all $k = 1, 2, \dots, K-1$, we can infer that it is universally applicable that

Hence,

Through the use of the discrete Gronwall inequality ^[28], the inequality (4.8) becomes

Combining with Lemma 3.2, the above estimate becomes

Using the same approach, we can derive the result

Combining (4.3), (4.7), and (4.9)–(4.10), we complete the proof of Theorem 4.1. □

5.   Numerical example

In this section, we will verify the correctness of the theory by using several numerical examples. We employ the proposed algorithm to solve the values of the VIEs with different values of $\psi(s)$ and the error and convergence order. All numerical examples will be implemented using Matlab software.

Example 5.1. The two-dimensional linear $\psi$-type VIEs under consideration are as follows.

where $(s, t)\in[a, b]\times[c, d]$, and

with the exact solution $\nu(s, t) = (\psi(s))^{4}(\psi(t))^{4}$.

In this paper, we give two groups of data where $\sigma_{1}$ and $\sigma_{2}$ take different values, which are $\sigma_{1} = 0.4$, $\sigma_{2} = 0.6$ and $\sigma_{1} = 0.3$, $\sigma_{2} = 0.5$, respectively. The steps in $s$-direction and $t$-direction are divided into $h_{s} = \frac{b-a}{2L}$, $h_{t} = \frac{d-c}{2K}$. Error is defined as follows:

when $\psi(s) = \log(s)$, and take $[a, b]\times[c, d] = [1, 2]\times[1, 2]$. To evaluate the level of convergence achieved by the numerical schemes developed for different values of $\sigma_{1}$ and $\sigma_{2}$, we conducted tests as presented in Tables 1 and 2. The convergence order was determined using the formula $log_{2}(\frac{e_{h}}{e_{h/2}})$. According to the analysis conducted in this study, the theoretical convergence order of the numerical scheme developed is expressed as $O(h_{s}^{3+\sigma_{1}}+h_{t}^{3+\sigma_{2}})$. It can be observed that when $h_{s}$ significantly exceeds $h_{t}$, the convergence order simplifies to $O(h_{s}^{3+\sigma_{1}})$. The tables presented in this paper (Table 1 and Table 2) consider different values for parameters $L$ and $K$. In Table 1, we set $K = 2L$ with a corresponding order of $3+\sigma_{1}$, while in Table 2, we set $L = 2K$ with a corresponding order of $3+\sigma_{2}$. The values chosen for these tables are as follows: $\sigma_{1} = 0.4, \sigma_{2} = 0.6$ and $\sigma_{1} = 0.3, \sigma_{2} = 0.5$. From our observations in Table 1, it can be inferred that when $\sigma_{1} = 0.4, \sigma_{2} = 0.6$, the achieved order closely approximates 3.4; similarly, when $\sigma_{1} = 0.3, \sigma_{2} = 0.5$, the achieved order is close to 3.3, which aligns well with the theoretical expectation of $3+\sigma_{1}$. Similar conclusions can also be drawn from Table 2, with the convergence order close to the theoretical order $3+\sigma_{2}$.

The Figure 1 below shows the error distribution of $\psi(s) = \log(s)$ when $K = 2L$ and $L = 128$. The left diagram shows the discrepancy distribution of $\sigma_{1} = 0.4, \sigma_{2} = 0.6$, and the right diagram shows the discrepancy distribution of $\sigma_{1} = 0.3, \sigma_{2} = 0.5$. From the following two error graphs of Figure 1, we can see that when $\sigma_{1} = 0.4, \sigma_{2} = 0.6$ and $\sigma_{1} = 0.3, \sigma_{2} = 0.5$, although the errors are different, the error graphs are similar, and the errors can reach $10^{-10}$, so it can be known from the error graphs that the numerical method used can better approximate $\nu(s, t)$.

When $\psi(s) = s^{\gamma}$, where $\gamma = 1, 2$, here we're going to take $\gamma = 2$ and take $[a, b]\times[c, d] = [0, 1]\times[0, 1]$. We list in the table the maximum errors and orders obtained for different values of $\sigma_{1}$, $\sigma_{2}$, $h_{s}$, and $h_{t}$. In this case, the values of $\sigma_{1}$ and $\sigma_{2}$ are the same as in the $\psi(s) = \log(s)$, and the values of $K$ and $L$ are going to be the same as $\psi(s) = \log(s)$. We can come to the same conclusion. When $K = 2L$, the orders of Table 3 are close to 3.4 and 3.3, respectively, which is close to the theoretical order $3+\sigma_{1}$; when $L = 2K$, the orders of Table 4 are close to 3.6 and 3.5, respectively. That is, close to the theoretical order $3+\sigma_{2}$.

The Figure 2 below shows the error distribution of $\psi(s) = s^{\gamma}$, where $\gamma = 2$ when $K = 2L$ and $L = 128$. The left diagram shows the discrepancy distribution of $\sigma_{1} = 0.4, \sigma_{2} = 0.6$, and the right diagram shows the discrepancy distribution of $\sigma_{1} = 0.3, \sigma_{2} = 0.5$. From the following two error graphs, we can see that when the values of $\sigma_{1}$ and $\sigma_{2}$ are different, although the errors are not the same, the error graphs are similar, the errors can reach $10^{-7}$, and the effect is better, so it can be seen from the error graph that the numerical method used can better approximate $\nu(s, t)$.

Example 5.2. Consider the following nonlinear equation:

where $(s, t)\in[a, b]\times[c, d]$, and

with the exact solution is $\nu(s, t) = (\psi(s))^{2}(\psi(t))^{3}$. We're still divided into two cases: first of all, $\psi(s) = \log(s)$, and second, $\psi(s) = s^{\gamma}$.

Now, let's consider the case $\psi(s) = \log(s)$, and take $[a, b]\times[c, d] = [1, 2]\times[1, 2]$. For this example, we repeat the calculation procedure of Example 5.1 using the appropriate scheme. The variation of its maximum error and order with $K = 2L$ or $L = 2K$ is shown in Table 5 or Table 6, and the values of $L$ and $K$ are the same as those in Example 5.1. These two tables again verify that the order is close to the theoretical order 3.4, 3.3, i.e., $3+\sigma_{1}$, and 3.6, 3.5, i.e., $3+\sigma_{2}$.

Now, let's think about $\psi(s) = s^{\gamma}$, $\gamma = 2$ and take $[a, b]\times[c, d] = [0, 1]\times[0, 1]$. From the following two tables, we know that the errors gradually become smaller and the order also tends to stabilize, which is close to 3.4, 3.3 and 3.6, 3.5.

6.   Conclusions

In this paper, the modified block-by-block method is used to solve the numerical solution of fractional $\psi$-Volterra integral equations. The convergence of order of the high order numerical scheme is analyzed rigorously by using the discrete Gronwall inequality. Through experimental verification, we find that this method has high accuracy and stability and can effectively deal with high dimensionality fractional $\psi$-Volterra integral equations.

We furthermore discussed the sophisticated numerical scheme of high order for high dimensionality fractional $\psi$-Volterra integral equations with singular solution by using graded mesh. We also noticed that when the grid division is very fine, the computational complexity of the format is very high. Therefore, we will plan to introduce the fast high-order algorithms for solving the high dimensionality fractional $\psi$-Volterra integral equations by using the fast Fourier transform.

Thus, we can apply it to a wider range of practical problems. There are still some limitations and shortcomings in this study, that is, the calculation amount is too big. In the future, we can continue to study that the modified block-by-block method can be further improved to enhance its computational efficiency and stability ^[29]. Additionally, the modified block-by-block method can be applied to practical problems and its application value in engineering and science can be explored.

Acknowledgments

Z. Q. Wang was supported by National Natural Science Foundation of China (NSFC) (Grant No. 11961009) and Foundation of Guizhou Science and Technology Department (Grant No. QHKJC-ZK[2024]YB497). J. Y. Cao was supported by National NSFC (Grant Nos. 12361083, 62341115) and Science research fund support project of the Guizhou Minzu University (Grant No. GZMUZK[2023]CXTD05). Z. Q. Wang and J. Y. Cao were supported by NSF of the Department of Education of Guizhou Province (Grant Nos. QJJ2023012, QJJ2023062, QJJ2023061).

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare that they have no competing interests.

Appendix A.   Construction of numerical scheme at points $(s_i, t_j)$ for $i=1, 2$ or $j=1, 2$.

For $\nu(s_{1}, t_{1})$, one has

where

Similar to (A.1), one can obtain $\nu(s_{2}, t_{1}), \nu(s_{1}, t_{2})$, and $\nu(s_{2}, t_{2})$ as follows:

where

Therefore, $\nu_{1, 1}$, $\nu_{2, 1}$, $\nu_{1, 2}$, and $\nu_{2, 2}$ can be simultaneously solved from (A.1), (A.4), (A.5), and (A.6).

For $\nu(s_{2l+r_{1}}, t_{r_{2}}), l = 1, 2, \dots, L-1$ and $\nu(s_{r_{1}}, t_{2k+r_{2}}), r_{1}, r_{2} = 1, 2, k = 1, 2, \dots, K-1$, under assumption that $\nu_{i, r_{2}}, i = 0, 1, \dots, 2l$ and $\nu_{r_{1}, j}, j = 0, 1, \dots, 2k$ already known, we have

where $B_{1}^{q, 0}$ is given by (A.3) and $A_{2l+1}^{p, 0}, A_{2l+1}^{p, i}$ are defined by (2.4) and (2.9).

For $\nu(s_{2l+2}, t_{1})$, $\nu(s_{2l+1}, t_{2})$, and $\nu(s_{2l+2}, t_{2})$, we have

where $B_{1}^{q, 0}$, $B_{2}^{q, 0}$, $A_{2l+1}^{p, 0}$, $A_{2l+1}^{p, i}$, and $A_{2l+2}^{p, i}$ are given by (A.3), (A.8), (2.4), (2.9), and (2.15), respectively.

For $\nu(s_{1}, t_{2k+1})$, similar to (A.10), one can obtain

where $A_{1}^{p, 0}$, $A_{2}^{p, 0}$ are given by (A.2), (A.7). $B_{2k+1}^{q, 0}$, $B_{2k+1}^{q, j}$, and $B_{2k+2}^{q, j}$ are defined by (2.5), (2.10), and (2.16).

Appendix B.   The proof of Lemma 3.2.

Proof. Let $Q_{1}, Q_{2}$ be defined as follows:

For $E_{2l+1, 2k+1}$, we have

According to the Taylor theorem, there is

where $(\varepsilon _{1}(\tau), \xi_{1}(\eta)) \in [a, s_{1}]\times[c, t_{1}]$.

For $(\tau, \eta) \in [a, s_{1}]\times[t_{2j-1}, t_{2j+1}]$, there is $(\varepsilon _{2}(\tau), \xi_{j1}(\eta)) \in [a, s_{1}]\times[t_{2j-1}, t_{2j+1}]$,

In the same way, we have $(\tau, \eta) \in [s_{2i-1}, s_{2i+1}]\times[c, t_{1}]$, and $(\varepsilon _{i1}(\tau), \xi_{2}(\eta)) \in [s_{2i-1}, s_{2i+1}]\times[c, t_{1}]$,

For $\forall (\tau, \eta) \in [s_{2i-1}, s_{2i+1}]\times[t_{2j-1}, t_{2j+1}]$, there is $(\varepsilon _{i2}(\tau), \xi_{j2}(\eta)) \in [s_{2i-1}, s_{2i+1}]\times[t_{2j-1}, t_{2j+1}]$,

By putting (B.4) into the first term to the right of formula (B.3) gives the following result:

Using direct and simple calculations, we have

Similarily,

By Lemma 3.1 and (B.9), we have

where $Q_{1}$ is defined by (B.1).

By Lemma 3.1 and (B.10), we have

According to (B.11) and (B.12), we get $E_{2l+1, 2k+1}^{1}$ as follows:

Bring equation (B.5) to the second term to the right of equation (B.3), and we have the following:

Using Lemma 3.1 and equation (B.9), we have

For $E_{2}^{2}$, we estimate this by adding one term and subtracting one term,

By definition of $\hat{f}_{p, i}(\tau)$ and $\tau \in (s_{i}, s_{i+2})$, $p = 0, 1, 2;i = 0, 1, \dots, 2l$, we know $|\hat{f}_{p, i}(\tau)| \le 1$; hence, $|\sum\limits_{p = 0}^{2}\hat{f}_{p, i}(\tau)| \le 3$. We get

Next, we will reckon $L_{11}$. To begin, we reckoned $\int_{t_{2j-1}}^{t_{2j}}\prod\limits_{q = 0}^{2}(\psi(\eta)-\psi(t_{2j-1+q}))d\psi(\eta)$. By the Lagrange mean value theorem, the integral mean value theorem, and continuity of the function on a closed interval, we have

where $\zeta_{i}^{j}\in(t_{2j-1}, t_{2j}), i = 0, 1, 2$; $\gamma_{j}\in(t_{2j-1}, t_{2j})$, $\widehat{\zeta_{j}}\in(t_{2j-1}, t_{2j})$.

Using the same method of (B.18), we estimate $\int_{t_{2j}}^{t_{2j+1}}\prod\limits_{q = 0}^{2}(\psi(\eta)-\psi(t_{2j-1+q}))d\psi(\eta)$ and obtain

Combining (B.18) and (B.19), we obtain

First, we will estimate $y1$,

where $\bar{\eta} _{j}\in (t_{2j-1}, t_{2j})$, $\bar{\bar{\eta}}_{j}\in (t_{2j}, t_{2j+1})$, and $\eta _{j}\in(\bar{\eta} _{j}, \bar{\bar{\eta}}_{j})$.

Second, we will estimate $y2$. If $\psi'(\eta)$ satisfies the Lipschitz condition, i.e.,

then, we have

Bring (B.21) and (B.22) into (B.20), and we have the estimation of $L_{11}$ as follows:

For $L_{12}$, using (3.3) and (B.10), we have

Combine (B.23) and (B.24) to (B.17). Therefore, $L1$ is as follows:

Using the mean value theorem and Lemma 3.1, we can estimate $L2$ as follows:

where $Q_{2}$ is defined in (B.2).

Combine (B.25) and (B.26) to obtain $E_{2}^{2}$,

Take (B.15), (B.27) and put it into (B.14), and we have

Next, we will estimate $E_{2l+1, 2k+1}^{3}$ and $E_{2l+1, 2k+1}^{4}$. Similar to the estimate for $E_{2l+1, 2k+1}^{1}$, $E_{2l+1, 2k+1}^{2}$, by putting (B.6) into the third term on the right side of formula (B.3), we have

Bring (B.7) to the fourth item on the right side of (B.3), and we have

Plug (B.13), (B.28), (B.29), and (B.30) into (B.3). We get

where $C$ is independent of $h _{s}$ and $h _{t}$.

The estimation of the other cases for $E_{i, j}$ are similar to $E_{2l+1, 2k+1}$. Here, we will omit them. Therefore, we have already proven Lemma 3.2.