High Accuracy Combination Method for Solving the Systems of Nonlinear Volterra Integral and Integro-Differential Equations with Weakly Singular Kernels of the Second Kind

Abstract

This paper presents a high accuracy combination algorithm for solving the systems of nonlinear Volterra integral and integro-differential equations with weakly singular kernels of the second kind. Two quadrature algorithms for solving the systems are discussed, which possess high accuracy order and the asymptotic expansion of the errors. By means of combination algorithm, we may obtain a numerical solution with higher accuracy order than the original two quadrature algorithms. Moreover an a posteriori error estimation for the algorithm is derived. Both of the theory and the numerical examples show that the algorithm is effective and saves storage capacity and computational cost.

1. Introduction

In this paper, we first consider the following system of differential equations: 𝑥^′_𝑖(𝑡)=𝑓_𝑖𝑥(𝑡,𝑥(𝑡),𝑧(𝑡)),_𝑖(0)=𝑥_𝑖0,𝑖=1,…,𝑚, (1.1) where 𝑥(𝑡)=(𝑥₁(𝑡),…,𝑥_𝑚(𝑡)), 𝑧(𝑡)=(𝑧₁₁(𝑡),…,𝑧_1𝑚(𝑡),…,𝑧_𝑚1(𝑡),…,𝑧_𝑚𝑚(𝑡)), and 𝑧_𝑖𝑗(𝑡)=∫^𝑡₀(𝑡−𝑠)^𝛼𝑖𝑗(ln(𝑡−𝑠))^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡,𝑠,𝑥(𝑠))𝑑𝑠, 0>𝛼_𝑖𝑗>−1, 𝛽_𝑖𝑗=0 or 1, 𝑖,𝑗=1,…,𝑚 as well as 𝑘_𝑖𝑗(𝑡,𝑠,𝑥), 𝑖,𝑗=1,…,𝑚, are continuous functions for 0≤𝑠≤𝑡≤𝑇 and 𝑥∈𝑅^𝑚.

The problem (1.1) can be transformed to the following system of nonlinear integral equations: find 𝑧(𝑡) and 𝑥(𝑡) satisfying 𝑧_𝑖𝑗(𝑡)=^𝑡₀(𝑡−𝑠)^𝛼𝑖𝑗(ln(𝑡−𝑠))^𝛽𝑖𝑗𝑘_𝑖𝑗𝑥(𝑡,𝑠,𝑥(𝑠))𝑑𝑠,𝑖,𝑗=1,…,𝑚,_𝑖(𝑡)=𝑥_𝑖0+^𝑡₀𝑓_𝑖(𝑠,𝑥(𝑠),𝑧(𝑠))𝑑𝑠,𝑖=1,…,𝑚, (1.2) which is a special case of nonlinear Volterra integral system with weakly singular kernels of the second kind 𝑢_𝑖(𝑡)=𝑦_𝑖(𝑡)+^𝑚_𝑗=1^𝑡₀(𝑡−𝑠)^𝛼𝑖𝑗(ln(𝑡−𝑠))^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡,𝑠,𝑢(𝑠))𝑑𝑠,𝑖=1,…,𝑚, (1.3) where −1<𝛼_𝑖𝑗<0, and 𝑘_𝑖𝑗(𝑡,𝑠,𝑢) is continuous function on 0≤𝑠≤𝑡≤𝑇 and 𝑢=(𝑢₁,…,𝑢_𝑚)∈𝑅^𝑚.

Nonlinear Volterra integral and integro-differential equations with weakly singular kernels of the second kind play important roles in the mathematical modeling of many physical and biological phenomena, particularly in such fields as heat transfer, nuclear reactor dynamics, and thermoelasticity, in which it is necessary to take into account the effect of past history. Plenty of work has been done to develop and analyze numerical methods for solving the Volterra integral and integro-differential equations with weakly singular kernels of the second kind; see [1–5] and reference therein. The recent progress in this research area has been achieved for the extrapolation method [6], the spline collocation method [7], and the Galerkin method [8, 9]. However, there are few works for solving the systems of these types of equations, which are more important than single equation for many applications. For example, the elastodynamic problems for piezoelectric and pyroelectric hollow cylinders under radial deformation can be successfully transformed into a system of two second kind Volterra integral equations.

The combination method as an accelerating convergence technique for solving integral equations was firstly presented in 1984 [10]. Similar to the extrapolation method, the combination method combines several approximations to obtain an approximation of higher accuracy. One important advantage of the extrapolation method and the combination method is parallel computation since those original approximations can be computed independently. However, the extrapolation method uses a coarse grid and some finer grids. We must do much more work on the finer grids than the coarse one, which lowers the degree of parallelism. On the other hand, the loads of computing the approximations in the combination method are close to each other. Hence the combination method is an efficient parallel method to obtain an approximation of high accuracy with a high degree of parallelism. Recently this method has been used to solve the first kind Abel integral equations [11]. In this paper, we will apply the high accuracy combination method for solving the systems of nonlinear Volterra integral and integro-differential equations with weakly singular kernels of the second kind, for which there are few results, even for the general convergence to be proved in Section 4 of this paper.

In this paper, based on Lyness's [12] modified mid-point rectangular quadrature rule and modified trapezoidal quadrature rule, we will construct two quadrature methods to solve the systems. By means of the asymptotic expansions of the errors for both algorithms, we present a high accuracy combination algorithm. Then by using the generalization of discrete Gronwall inequality [13], the convergence rate, stability, and the asymptotic expansion of the error of the combination approximate solution are proved. Comparing to other algorithms, for example, extrapolation methods [6, 14], the combination algorithm has the advantage of less computation complexity because both of the modified mid-point rectangular quadrature rule and the modified trapezoidal quadrature rule have the same step size. Moreover an a posteriori error estimation is obtained, by which we can rectify the accuracy of our algorithm in processing.

2. Existence and Uniqueness of the Solution

Since the uniqueness, existence, and numerical methods of (1.2) can be decided by (1.1), we only discuss the problem (1.3) in the following.

Let 𝑘_𝑖𝑗(𝑡,𝑠,𝑢) satisfy Lipschitz condition (A): |𝑘_𝑖𝑗(𝑡,𝑠,𝑢)−𝑘_𝑖𝑗(𝑡,𝑠,𝑣)|≤𝐿‖𝑢−𝑣‖_∞∀𝑢,𝑣∈𝑅^𝑚,0≤𝑠≤𝑡≤𝑇,1≤𝑖,𝑗≤𝑚, (2.1) then there is a unique solution in (1.3). In fact 𝛼=min_𝑖,𝑗𝛼_𝑖𝑗>−1, 𝑢(𝑡)=(𝑢₁(𝑡),…,𝑢_𝑚(𝑡))^𝑇, and ‖𝑢(𝑡)‖_∞=max_{1≤𝑖≤𝑚}|𝑢_𝑖(𝑡)|. We can choose 𝜂>0 such that 𝛼−𝜂>−1. Thus we have 𝑁=max_{1≤𝑖,𝑗≤𝑚,0<𝑡<𝑇}𝑡^{𝛼𝑖𝑗−𝛼+𝜂}|(ln𝑡)^𝛽𝑖𝑗|<∞. (2.2) Let (𝐶_[0,𝑇])^𝑚=𝐶_[0,𝑇]×⋯×𝐶_[0,𝑇] be a continuous function space which maps [0,𝑇] into 𝑅^𝑚 and let 𝐹 be a mapping from (𝐶_[0,𝑇])^𝑚 to itself. Then we have 𝐹_𝑖(𝑢)(𝑡)=𝑦_𝑖(𝑡)+^𝑡₀^𝑚_𝑗=1(𝑡−𝑠)^𝛼𝑖𝑗(ln(𝑡−𝑠))^𝛽_𝑖𝑗𝑘_𝑖𝑗(𝑡,𝑠,𝑢(𝑠))𝑑𝑠,1≤𝑖≤𝑚, (2.3) where 𝑦=(𝑦₁,…,𝑦_𝑚)∈(𝐶_[0,𝑇])^𝑚, and |𝐹_𝑖(𝑢)(𝑡)−𝐹_𝑖=(𝑣)(𝑡)|_{𝑚𝑗=1𝑡0}(𝑡−𝑠)^𝛼−𝜂(𝑡−𝑠)^{𝛼𝑖𝑗−𝛼+𝜂}|(ln(𝑡−𝑠))^𝛽_𝑖𝑗||𝑘_𝑖𝑗(𝑡,𝑠,𝑢(𝑠))−𝑘_𝑖𝑗≤(𝑡,𝑠,𝑣(𝑠))|𝑑𝑠_{𝑚𝑗=1𝑡0}(𝑡−𝑠)^𝛼−𝜂𝑁𝐿|𝑢(𝑠)−𝑣(𝑠)|_∞𝑑𝑠=𝑚𝑁𝐿^𝑡₀(𝑡−𝑠)^𝛼−𝜂|𝑢(𝑠)−𝑣(𝑠)|_∞𝑑𝑠,1≤𝑖≤𝑚, (2.4) where 𝑓(𝑡)=𝑡^𝛼−𝜂∈𝐿^′_[0,𝑇], which is due to 𝛼−𝜂>−1.

Lemma 2.1. Suppose that condition (A) is satisfied, then there is a unique solution to (1.3).

Proof. Assume that 𝑢(𝑡)=(𝑢₁(𝑡),…,𝑢_𝑚(𝑡))^𝑇 and 𝑣(𝑡)=(𝑣₁(𝑡),…,𝑣_𝑚(𝑡))^𝑇 are two different solutions to (1.3). Defining 𝑤(𝑡)=𝑢(𝑡)−𝑣(𝑡), we get 𝑤_𝑖(𝑡)=^𝑚_𝑗=1^𝑡₀(𝑡−𝑠)^𝛼𝑖𝑗(ln|𝑡−𝑠|)^𝛽𝑖𝑗(𝑘_𝑖𝑗(𝑡,𝑠,𝑢(𝑠))−𝑘_𝑖𝑗(𝑡,𝑠,𝑣(𝑠)))𝑑𝑠,𝑖=1,…,𝑚,‖𝑤(𝑡)‖_∞≤𝑚𝑁𝐿^𝑡₀𝑡^𝛼−𝜂‖𝑤(𝑠)‖_∞𝑑𝑠. (2.5) Therefore using Gronwall inequality, we get 𝑤(𝑡)=0, which leads to the uniqueness.
In order to prove existence, we use a simple iterative process: let 𝑢⁽⁰⁾(𝑡)=0, 𝑢⁽¹⁾(𝑡)=𝑦(𝑡) and 𝑢_𝑖^(𝑛)(𝑡)=𝑦_𝑖(𝑡)+^𝑡₀^𝑚_𝑗=1(𝑡−𝑠)^𝛼𝑖𝑗(ln(𝑡−𝑠))^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡,𝑠)𝑢^(𝑛−1)(𝑠)𝑑𝑠,𝑛=1,2,…. (2.6) Now we will prove that 𝑢^(𝑛)(𝑡) is a convergent sequence. Let 𝑉^(𝑛)(𝑡)=‖𝑢^(𝑛)(𝑡)−𝑢^(𝑛−1)(𝑡)‖_∞, (2.7) then from (2.6) and (2.4) we get 𝑉^(𝑛+1)(𝑡)≤𝑀^𝑡₀(𝑡−𝑠)^𝛼−𝜂𝑉^(𝑛)(𝑠)𝑑𝑠=𝑀𝑡^𝛽∗𝑉^(𝑛)=𝑀^𝑛𝑡^𝛽∗⋯∗𝑡^𝛽∗𝑉⁽¹⁾, (2.8) where 𝑡^𝛽∗𝑉^(𝑛) denotes convolution, 𝑀=𝑚𝑁𝐿, and 𝛽=𝛼−𝜂.
Taking Laplace transformation, we can easily deduce that 𝑉^(𝑛+1)(𝑡)≤(𝑀Γ(𝛽+1)𝑠^𝛽+1)^𝑛𝑉⁽¹⁾(𝑠). (2.9) Taking inverse Laplace transformation, we get 𝑉^(𝑛+1)(𝑡)≤(𝑀Γ(𝛽+1))^𝑛Γ(𝑛𝛽+𝑛)^𝑡₀𝜏^𝑛(𝛽+1)𝑉⁽⁰⁾(𝑡−𝜏)𝑑𝜏. (2.10) Note that for 𝑛>𝑚, we have 𝑢^(𝑛)(𝑡)−𝑢^(𝑚)(𝑡)=(𝑢^(𝑛)(𝑡)−𝑢^(𝑛−1)(𝑡))+(𝑢^(𝑛−1)(𝑡)−𝑢^(𝑛−2)(𝑡))+⋯+(𝑢^(𝑚+1)(𝑡)−𝑢^(𝑚)(𝑡)) (2.11) or ‖𝑢^(𝑛)(𝑡)−𝑢^(𝑚)(𝑡)‖_∞≤𝑉^(𝑛)(𝑡)+⋯+𝑉^(𝑚+1)(𝑡). (2.12) But using (2.10), we easily prove that ^∞_𝑛=0𝑉^(𝑛+1)(𝑡)≤^∞_𝑛=0(𝑀Γ(𝛽+1))^𝑛Γ(𝑛𝛽+𝑛)^𝑡₀𝜏^𝑛(𝛽+1)𝑉⁽⁰⁾(𝑡−𝜏)𝑑𝜏 (2.13) is a convergent series, which means that lim_{𝑛,𝑚→∞}‖𝑢^(𝑛)(𝑡)−𝑢^(𝑚)(𝑡)‖_∞=0or𝑢^(𝑛)(𝑡)→𝑢(𝑡) (2.14) and 𝑢(𝑡) is the solution to (1.2).

3. The Numerical Methods

In this section two quadrature algorithms which are based on [12] will be given to solve the systems. Consider the weakly singular integral 𝐼(𝐺)=^𝑏_𝑎(𝑏−𝑥)^𝛼𝑔(𝑥)𝑑𝑥,𝐺(𝑥)=(𝑏−𝑥)^𝛼𝑔(𝑥) (3.1) and Navot's modified trapezoidal rule [15] 𝑄^ℎ_𝑇ℎ(𝐺)=2𝐺(𝑎)+ℎ^𝑁−1_𝑗=1𝐺(𝑥_𝑗)−𝜁(−𝛼)ℎ^1+𝛼𝑔(𝑏). (3.2) There is an error estimate as follows: 𝐸^ℎ_𝑇(𝐺)=𝑄^ℎ_𝑇(𝐺)−𝐼(𝐺)=−𝜁(−𝛼−1)ℎ^2+𝛼𝑔^′(𝑏)+𝑂(ℎ²). (3.3) We differentiate with respect to 𝛼 in (3.3) and get 𝐸^ℎ_𝑇(𝐺^′)=𝑄^ℎ_𝑇(𝐺^′)−𝐼(𝐺^′)=[𝜁^′(−𝛼−1)−𝜁(−𝛼−1)lnℎ]ℎ^2+𝛼𝑔^′(𝑏)+𝑂(ℎ²). (3.4) Then we have 𝐼(𝐺^′)=^𝑏_𝑎(𝑏−𝑥)^𝛼(ln(𝑏−𝑥))^𝛽𝑔(𝑥)𝑑𝑥,𝐺^′(𝑥)=(ln(𝑏−𝑥))^𝛽(𝑏−𝑥)^𝛼𝑔(𝑥),𝛽=0,1 (3.5) such that 𝑄^ℎ_𝑇(𝐺^′ℎ)=2𝐺^′(𝑎)+ℎ^𝑁−1_𝑗=1𝐺^′(𝑥_𝑗)−[−𝛽𝜁^′(−𝛼)+𝜁(−𝛼)(lnℎ)^𝛽]ℎ^1+𝛼𝑔(𝑏). (3.6) Thus we have 𝐸^ℎ_𝑇(𝐺^′)=𝑄^ℎ_𝑇(𝐺^′)−𝐼(𝐺^′)=[𝛽𝜁^′(−𝛼−1)−𝜁(−𝛼−1)(lnℎ)^𝛽]ℎ^2+𝛼𝑔^′(𝑏)+𝑂(ℎ²). (3.7) Equations (3.5) and (3.6) are the integral functions with logarithm singularity and their modified trapezoidal rule, and (3.7) is an asymptotic expansion of the error.

From [16], we have the modified mid-point rectangular rule 𝑄^ℎ_𝑀(𝐺)=ℎ^𝑁−1_𝑗=0𝐺(𝑥_𝑗+1/2)−(2^−𝛼−1)𝜁(−𝛼)ℎ^1+𝛼𝑔(𝑏). (3.8) Hence we get 𝐸^ℎ_𝑀(𝐺)=𝑄^ℎ_𝑀(𝐺)−𝐼(𝐺)=−(2^−𝛼−1)𝜁(−𝛼−1)ℎ^2+𝛼𝑔^′(𝑏)+𝑂(ℎ²). (3.9) From (3.4)–(3.7), we get 𝐸^ℎ_𝑀(𝐺^′)=𝑄^ℎ_𝑀(𝐺^′)−𝐼(𝐺^′)=−[−𝛽2^−𝛼−1ln2𝜁(−𝛼−1)−𝛽(2^−𝛼−1−1)𝜁^′(−𝛼−1)+(2^−𝛼−1−1)𝜁(−𝛼−1)(lnℎ)^𝛽]ℎ^2+𝛼𝑔^′(𝑏)+𝑂(ℎ²). (3.10) More generally we have 𝑄^ℎ_𝑀(𝐺^′)=ℎ^𝑁−1_𝑗=0𝐺^′(𝑥_𝑗+1/2)−[−𝛽2^−𝛼ln2𝜁(−𝛼)−𝛽(2^−𝛼−1)𝜁^′(−𝛼)+(2^−𝛼−1)𝜁(−𝛼)(lnℎ)^𝛽]ℎ^1+𝛼𝑔(𝑏). (3.11) Therefore from (3.7) and (3.10), if 𝛽=0, we get (2^−𝛼−1−1)𝐸^ℎ_𝑇(𝐺^′)−𝐸^ℎ_𝑀(𝐺^′)=𝑂(ℎ²) (3.12) or (1−2^−𝛼−12−2^−𝛼−1𝑄^ℎ_𝑇(𝐺^′1)+2−2^−𝛼−1𝑄^ℎ_𝑀(𝐺^′))−𝐼(𝐺^′)=𝑂(ℎ²). (3.13) Note that the combination has a high accuracy order 𝑂(ℎ²), which is higher than 𝑂(ℎ^2+𝛼) in (3.7) and (3.10).

If 𝛼=0, 𝛽=1 in (3.5), then 𝑄^ℎ_𝑇(𝐺^′ℎ)=2𝐺^′(𝑎)+ℎ^𝑁−1_𝑗=1𝐺^′(𝑥_𝑗1)+2ℎln(𝑄2𝜋)ℎ𝑔(𝑏),_ℎ𝑀(𝐺^′)=ℎ^𝑁−1_𝑗=0𝐺^′(𝑥_𝑗+1/21)−2ln2ℎ𝑔(𝑏), (3.14) where 𝜁(0)=−0.5, 𝜁^′(0)=−ln(2𝜋)/2, and (3.13) becomes (13𝑄^ℎ_𝑇(𝐺^′2)+3𝑄^ℎ_𝑀(𝐺^′))−𝐼(𝐺^′)=𝑂(ℎ²). (3.15) Using (3.14) and (3.15), we can construct two algorithms for solving Volterra integral system of equations. Since the kernels of the systems are weakly singular, the following Navot's quadrature rule is used. Setting 𝑡=𝑡_𝑙, we get 𝑢_𝑖(𝑡_𝑙)=𝑦_𝑖(𝑡_𝑙)+^𝑚_𝑗=1^𝑡𝑙₀(𝑡_𝑙−𝑠)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑠|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑠,𝑢(𝑠))𝑑𝑠,𝑖=1,2,…,𝑚. (3.16) Now we recall the following lemma from [8].

Lemma 3.1. Let 𝑔(𝑡)∈𝐶^2𝑟[𝑎,𝑏], 𝐺(𝑥)=(𝑏−𝑡)^𝛼(ln|𝑏−𝑡|)^𝛽𝑔(𝑡), ℎ=(𝑏−𝑎)/𝑁, 𝑡_𝑘=𝑎+𝑘ℎ(𝑘=0,…,𝑁), then modified trapezoidal rule 𝑇_𝑁ℎ(𝐺)=2𝐺(𝑡₀)+ℎ^𝑁−1_𝑗=1−[−𝛽𝜁^′(−𝛼)+𝜁(−𝛼)(lnℎ)^𝛽]𝑔(𝑏)ℎ^1+𝛼 (3.17) has an asymptotic error expansion 𝐸_𝑁(𝐺)=^𝑟−1_𝑗=1𝐵_2𝑗𝐺(2𝑗)!^(2𝑗−1)(𝑎)ℎ^2𝑗+^2𝑟−1_𝑗=1(−1)^𝑗[−𝛽𝜁^′(−𝛼−𝑗)+𝜁(−𝛼−𝑗)(lnℎ)^𝛽]𝑔^(𝑗)ℎ^𝑗+𝛼+1𝑗!+𝑂(ℎ^2𝑟), (3.18) where 𝛼>−1, 𝛽=0,1, 𝜁(𝑡) and 𝜁^′(𝑡) are Riemann-Zeta function and its derivative function, and 𝐵_2𝑗 are Bernoulli numbers.

From the above, we obtain the following discrete system of equations for modified trapezoidal quadrature method: 𝑢_𝑖0=𝑦_𝑖(𝑡₀𝑢),_𝑖𝑙=𝑦_𝑖(𝑡_𝑙)+^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗×(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘)+𝑤_{𝑙𝑙,𝑖𝑗}ℎ^{1+𝛼𝑖𝑗}𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙)],𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁, (3.19) where ℎ=1/𝑁, 𝑡_𝑘=𝑘ℎ, 𝑈_𝑗=(𝑢_1𝑗,𝑢_2𝑗,…,𝑢_𝑚𝑗), 𝑗=1,2,…,𝑚, and 𝑤_𝑙𝑘=12𝛽if𝑘=0,1if1≤𝑘≤𝑙−1,_𝑖𝑗𝜁^′(−𝛼_𝑖𝑗)+𝜁(−𝛼_𝑖𝑗)(lnℎ)^𝛽𝑖𝑗−𝛽_𝑖𝑗𝜁^′(−𝛼_𝑖𝑗)if𝑘=𝑙,𝑙=1,2,…,𝑁. (3.20)

On the other hand, we also obtain the discrete system for modified mid-point rectangular quadrature method 𝑢_𝑖0=𝑦_𝑖(𝑡₀𝑢),_𝑖𝑙=𝑦_𝑖(𝑡_𝑙)+^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘+1/2)^𝛼𝑖𝑗×((ln|𝑡_𝑙−𝑡_𝑘+1/2|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘+𝑈_𝑘+1/2)2)+𝑤_{𝑙𝑙,𝑖𝑗}ℎ^{1+𝛼𝑖𝑗}𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙)],𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁, (3.21) where 𝑤_𝑙𝑘=1if0≤𝑘≤𝑙−1,−[−𝛽_𝑖𝑗2^−𝛼ln2𝜁(−𝛼)−𝛽_𝑖𝑗(2^−𝛼−1)𝜁^′(−𝛼)if𝑘=𝑙,𝑙=1,2,…,𝑁,+(2^−𝛼−1)𝜁(−𝛼)(lnℎ)^𝛽𝑖𝑗]𝛽_𝑖𝑗=0,or1. (3.22)

Because the discrete system is nonlinear diagonal system of equations, we introduce the following iterative algorithms.

Algorithm 3.2 (modified trapezoidal quadrature method). We have the following steps.Step 1. Take 𝜀>0 sufficiently small and let 𝑢⁰(𝑡_𝑗)=(𝑢⁰_1𝑗,…,𝑢⁰_𝑚𝑗)=(𝑦₁(𝑡_𝑗),…,𝑦_𝑚(𝑡_𝑗)),𝑗=1,2,…,𝑁 (3.23) and 𝑛∶=0.Step 2. Compute 𝑢^𝑛+1_𝑖𝑙(𝑙≤𝑁,𝑖=1,2,…,𝑚) in parallel by the following simple iteration: 𝑢^𝑛+1_𝑖𝑙=𝑦_𝑖(𝑡_𝑙)+^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘)+𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈^𝑛_𝑙)], (3.24) where 𝑈_𝑘=(𝑢_1𝑘,…,𝑢_𝑚𝑘), 𝑈^𝑛_𝑙=(𝑢^𝑛_1𝑙,…,𝑢^𝑛_𝑚𝑙), and 𝑤_𝑙𝑘 is defined by (3.20).Step 3. If max_{1≤𝑖≤𝑚}|𝑢^𝑛+1_𝑖𝑙−𝑢^𝑛_𝑖𝑙|≤𝜀, then let 𝑢^𝑛+1_𝑖𝑙=𝑢_𝑖𝑙 and stop the iteration, otherwise set 𝑛∶=𝑛+1, go to Step 2.

Algorithm 3.3 (modified mid-point rectangular quadrature method). We have the following steps.Step 1. Take 𝜀>0 sufficiently small and let 𝑢⁰(𝑡_𝑗)=(𝑢⁰_1𝑗,…,𝑢⁰_𝑚𝑗)=(𝑦₁(𝑡_𝑗),…,𝑦_𝑚(𝑡_𝑗)),𝑗=1,2,…,𝑁 (3.25) and 𝑛∶=0.Step 2. Compute 𝑢^𝑛+1_𝑖𝑙(𝑙≤𝑁,𝑖=1,2,…,𝑚) in parallel by the following simple iteration: 𝑢^𝑛+1_𝑖𝑙=𝑦_𝑖(𝑡_𝑙)+^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗×𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,(𝑈^𝑛_𝑘+𝑈^𝑛_𝑘+1/2)2)+𝑤_{𝑙𝑙,𝑖𝑗}ℎ^{1+𝛼𝑖𝑗}𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈^𝑛_𝑙)], (3.26) where 𝑈_𝑘=(𝑢_1𝑘,…,𝑢_𝑚𝑘), 𝑈^𝑛_𝑙=(𝑢^𝑛_1𝑙,…,𝑢^𝑛_𝑚𝑙), and 𝑤_𝑙𝑘 is defined by (3.22).Step 3. If max_{1≤𝑖≤𝑚}|𝑢^𝑛+1_𝑖𝑙−𝑢^𝑛_𝑖𝑙|≤𝜀, then let 𝑢^𝑛+1_𝑖𝑙=𝑢_𝑖𝑙 and stop the iteration, otherwise set 𝑛∶=𝑛+1, go to Step 2.

4. Convergence and Error Estimation

In this section we will analyze the convergence of Algorithm 3.2 proposed in Section 3. The proof of convergence rate of Algorithm 3.3 is similar. When 𝑡=𝑡_𝑙, the system can be expressed as 𝑢_𝑖(𝑡_𝑙)=𝑦_𝑖(𝑡_𝑙)+^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘)+𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙+)]_𝑚𝑗=1[𝐸_{𝑖𝑗,𝑙,𝑡}((𝑡_𝑙−𝑡)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡,𝑢(𝑡)))]. (4.1) By Lemma 3.1, the remainder has the following estimate: ^𝑚_𝑗=1[𝐸_{𝑖𝑗,𝑙,𝑡}((𝑡_𝑙−𝑡)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡,𝑢(𝑡)))]=𝑂(ℎ^2+𝛼(lnℎ)^𝛽), (4.2) where 𝛼=min_{1≤𝑖,𝑗≤𝑚}𝛼_𝑖𝑗>−1 and 𝛽=max_{1≤𝑖,𝑗≤𝑚}𝛽_𝑖𝑗.

Letting 𝑒_𝑖𝑙=𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙 and {𝑢_𝑖𝑙} be the solution of the discrete equations, we derive 𝑒_𝑖0𝑒=0,_𝑖𝑙=^𝑚_𝑗=1^𝑙−1_𝑘=0ℎ𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗×(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑢₁(𝑡_𝑘),…,𝑢_𝑚(𝑡_𝑘))−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑢_1𝑘,…,𝑢_𝑚𝑘+))_𝑚𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑢₁(𝑡_𝑙),…,𝑢_𝑚(𝑡_𝑙))−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑢_1𝑙,…,𝑢_𝑚𝑙+))_𝑚𝑗=1[𝐸_{𝑖𝑗,𝑙,𝑡}((𝑡_𝑙−𝑡)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡,𝑢₁(𝑡),…,𝑢_𝑚(𝑡)))],𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁. (4.3)

Lemma 4.1. Suppose that the sequence {𝑒_𝑖}^𝑁_𝑖=1 satisfies 𝑒₀=0,|𝑒_𝑖|≤^𝑖−1_𝑗=1𝐵_𝑖𝑗|𝑒_𝑗|+𝐴,1≤𝑖≤𝑁, (4.4) where 𝐵_𝑖𝑗=2𝐿ℎ(𝑡_𝑖−𝑡_𝑗)^𝛼(ln|𝑡_𝑖−𝑡_𝑗|)^𝛽, −1<𝛼≤0, 𝛽=0,1, and ℎ is sufficiently small such that 𝐿ℎ𝑤_𝑖𝑖≤1/2. Then |𝑒_𝑖|≤𝐻𝐴, (4.5) where 𝐻=∑^∞_𝑘=0𝑅^𝑘/(𝑘!)^𝑠, 𝑅=2𝐿(𝑏−𝑎)^𝑠Γ(𝑠)𝑒^1/12𝑠(𝑒/𝑠)ℎ^𝑠, and 𝑠=1+𝛼.

Theorem 4.2. Assume that ℎ is sufficiently small, then the system of the nonlinear discrete equation (3.19) has a unique solution and the simple iteration (3.24) is geometrically convergent.

Proof. Firstly, if 𝑈_𝑙=(𝑢_1𝑙,𝑢_2𝑙,…,𝑢_𝑚𝑙) and 𝑉_𝑙=(𝑣_1𝑙,𝑣_2𝑙,…,𝑣_𝑚𝑙)(𝑙=1,2,…,𝑁) are solutions to (3.19), then {𝑍_𝑙=𝑈_𝑙−𝑉_𝑙} satisfies the inequality |𝑢_𝑖𝑙−𝑣_𝑖𝑙|≤^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗|𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑉_𝑘+)|]_𝑚𝑗=1[𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}|𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑉_𝑙≤)|]_𝑚𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗𝐿‖𝑈_𝑘−𝑉_𝑘‖_∞+𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}𝐿‖𝑈_𝑙−𝑉_𝑙‖_∞],𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁, (4.6) where we use 𝑢_𝑖0=𝑣_𝑖0=𝑦(𝑡_𝑖). Note that ‖𝑍_𝑘‖_∞=max_{1≤𝑗≤𝑚}|𝑧_𝑗𝑘|, then we easily deduce that ‖𝑍₀‖_∞=0,‖𝑍_𝑙‖_∞≤𝑀₁ℎ^𝑙−1_𝑘=0(𝑡_𝑙−𝑡_𝑘)^𝛼(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽‖𝑍_𝑘‖_∞,𝑙=1,2,…,𝑁, (4.7) where 𝑀₁=𝑚max_{0≤𝑘<𝑙≤𝑁}|𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^{𝛼𝑖𝑗−𝛼}(1+𝑐)(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝐿_𝑖𝑗|.
If ℎ is sufficiently small such that 𝑀₁ℎ^1+𝛼≤1/2, we have ‖𝑍₀‖_∞=0;‖𝑍_𝑙‖_∞≤2𝑀₁ℎ^𝑙−1_𝑘=0(𝑡_𝑙−𝑡_𝑘)^𝛼(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽‖𝑍_𝑘‖_∞,𝑙=1,2,…,𝑁. (4.8) Then by Lemma 4.1, we get ‖𝑍_𝑙‖_∞=0, 𝑙=1,2,…,𝑁. Hence the uniqueness is shown.
Secondly, from the iteration (3.24) we have |𝑢^𝑛+1_𝑖𝑙−𝑢^𝑛_𝑖𝑙|=^𝑚_𝑗=1[𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}|𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈^𝑛_𝑙)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙^𝑛−1≤)|]_𝑚𝑗=1[𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}𝐿‖𝑈^𝑛_𝑙−𝑈_𝑙^𝑛−1‖_∞],𝑖=1,2,…,𝑚;𝑙=1,2,…,𝑁. (4.9) Then ‖𝑈_𝑙^𝑛+1−𝑈^𝑛_𝑙‖_∞≤𝑀ℎ^1+𝛼‖𝑈^𝑛_𝑙−𝑈_𝑙^𝑛−1‖_∞≤12‖𝑈^𝑛_𝑙−𝑈_𝑙^𝑛−1‖_∞,𝑙=1,2,…,𝑁, (4.10) where 𝑀=𝑚𝑤_𝑙𝑙𝐿. We assume that ℎ is small enough such that 𝑀ℎ^1+𝛼≤1/2. Thus we prove that the simple iteration (3.24) is geometrically convergent, and its limit is the unique solution of (3.19).

Theorem 4.3. There is a positive constant 𝐶 independent of ℎ such that ‖𝑈_𝑙−𝑈_𝑙‖_∞≤𝐶𝜀ℎ^1+𝛼, (4.11) where 𝑈_𝑙=(𝑢_1𝑙,𝑢_2𝑙,…,𝑢_𝑚𝑙)=(𝑢^𝑛+1_1𝑙,𝑢^𝑛+1_2𝑙,…,𝑢^𝑛+1_𝑚𝑙) is defined in Algorithm 3.2.

Proof. Letting 𝑉_𝑙=𝑈_𝑙−𝑈_𝑙, we get 𝑢_𝑖0−𝑢_𝑖0𝑢=0,_𝑖𝑙−𝑢_𝑖𝑙=^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘+))]_𝑚𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈^𝑛_𝑙=))_𝑚𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘+))]_𝑚𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙^𝑛+1+))_𝑚𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙^𝑛+1)−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈^𝑛_𝑙)),𝑖=1,2,…,𝑚;𝑙=1,2,…,𝑁. (4.12) By Algorithm 3.2 and Lipschitz condition, we have ‖𝑉_𝑙‖_∞≤𝑀₁ℎ^𝑙−1_𝑘=0(𝑡_𝑙−𝑡_𝑘)^𝛼(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽‖𝑉_𝑘‖_∞+2𝑀₁𝜀ℎ^1+𝛼. (4.13) Using Lemma 3.1, we get ‖𝑈_𝑙−𝑈_𝑙‖_∞≤𝐶𝜀ℎ^1+𝛼.

Theorem 4.4. There is a positive constant 𝐶 independent of ℎ such that the errors 𝑒_𝑖𝑙=𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙(𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁) have the following error bound: max_{0≤𝑙≤𝑁;1≤𝑖≤𝑚}|𝑒_𝑖𝑙|≤𝐶ℎ^2+𝛼|lnℎ|^𝛽. (4.14)

Proof. Taking 𝐸_𝑙=max_{1≤𝑖≤𝑚}|𝑒_𝑖𝑙|, we get 𝐸₀𝐸=0,_𝑙≤𝑀₁ℎ^𝑙−1_𝑘=0(𝑡_𝑙−𝑡_𝑘)^𝛼(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝐸_𝑘+𝑀₁ℎ^1+𝛼𝐸_𝑙+𝑂(ℎ^2+𝛼),𝑙=1,2,…,𝑁. (4.15) If ℎ is sufficiently small such that 𝑀₁ℎ^1+𝛼≤1/2, we have 𝐸_𝑙≤2𝑀₁ℎ^𝑙−1_𝑘=0(𝑡_𝑙−𝑡_𝑘)^𝛼(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝐸_𝑘+𝑂(ℎ^2+𝛼),𝑙=1,2,…,𝑁. (4.16) Using inequality (4.5), we can get max_{0≤𝑙≤𝑁;1≤𝑖≤𝑚}|𝑒_𝑖𝑙|≤𝐶ℎ^2+𝛼|lnℎ|^𝛽. (4.17)

Corollary 4.5. Assume that 𝜀=𝑂(ℎ) in Algorithm 3.2, one can obtain the estimate max_{1≤𝑖≤𝑚;1≤𝑙≤𝑁}|𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙|=𝑂(ℎ^2+𝛼(lnℎ)^𝛽). (4.18)

5. Asymptotic Expansion, Combination, and a Posteriori Error Estimate

In the following we only derive the asymptotic expansions of the errors and the a posteriori error estimation of Algorithm 3.2. For Algorithm 3.3 we just simply present the corresponding result.

From Lemma 3.1 and (4.1) by using Taylor's expansion, we have 𝑢_𝑖(𝑡_𝑙)=𝑦_𝑖(𝑡_𝑙)+^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘)+𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙)]+ℎ^2+𝛼(lnℎ)^𝛽𝑚_𝑗=1𝑇⁰_𝑖𝑗(𝑡_𝑙)+ℎ^𝑚2+𝛼_𝑗=1𝑇¹_𝑖𝑗(𝑡_𝑙)+𝑂(ℎ^2+𝛼1|lnℎ|^𝛽1),𝑖=1,…,𝑚, (5.1) where 𝛼₁=min_{1≤𝑖,𝑗≤𝑚}{𝛼_𝑖𝑗∶𝛼_𝑖𝑗>𝛼}, 𝛽₁=max_{1≤𝑖,𝑗≤𝑚}{𝛽_𝑖𝑗∶𝛼_𝑖𝑗=𝛼₁}, and 𝑇⁰_𝑖𝑗𝑑(𝑡)=−𝜁(−𝛼−1)𝑘𝑑𝑡^𝑖𝑗(𝑡_𝑙,𝑡,𝑢(𝑡))if𝛼_𝑖𝑗=𝛼,0if𝛼_𝑖𝑗𝑇>𝛼,_1𝑖𝑗(𝑡)=𝛽𝜁^′𝑑(−𝛼−1)𝑘𝑑𝑡^𝑖𝑗(𝑡_𝑙,𝑡,𝑢(𝑡))if𝛼_𝑖𝑗=𝛼,0if𝛼_𝑖𝑗>𝛼. (5.2)

For 𝑒_𝑖𝑙=𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙, we have 𝑒_𝑖0=𝑢_𝑖(𝑡₀)−𝑢_𝑖0𝑒=0,_𝑖𝑙=𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙=^𝑚_𝑗=1^𝑙−1_𝑘=0ℎ𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^𝛽𝑖𝑗(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑢(𝑡_𝑘))−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑈_𝑘+))_𝑚𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}(𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑢(𝑡_𝑙))−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑙,𝑈_𝑙))+ℎ^2+𝛼(lnℎ)^𝛽𝑇⁰_𝑖(𝑡_𝑙)+ℎ^2+𝛼𝑇¹_𝑖(𝑡_𝑙)+𝑂(ℎ^2+𝛼1|lnℎ|^𝛽1),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁. (5.3) Here ∑^𝑚_𝑗=1𝑇^𝑘_𝑖𝑗=𝑇^𝑘_𝑖(𝑡_𝑙), 𝑘=0,1.

Using Theorem 4.4 and Taylor's expansion, we get 𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑢₁(𝑡_𝑘),…,𝑢_𝑚(𝑡_𝑘))−𝑘_𝑖𝑗(𝑡_𝑙,𝑡_𝑘,𝑢_1𝑘,…,𝑢_𝑚𝑘)=𝑘_{𝑖𝑗,𝑢1}(𝑡_𝑙,𝑡_𝑘,𝑢₁(𝑡_𝑘),…,𝑢_𝑚(𝑡_𝑘))(𝑢₁(𝑡_𝑘)−𝑢_1𝑘)+⋯+𝑘_{𝑖𝑗,𝑢𝑚}(𝑡_𝑙,𝑡_𝑘,𝑢₁(𝑡_𝑘),…,𝑢_𝑚(𝑡_𝑘))(𝑢_𝑚(𝑡_𝑘)−𝑢_𝑚𝑘)+𝑜(ℎ²)=𝑘_{𝑖𝑗,𝑢1}(𝑡_𝑙,𝑡_𝑘,𝑢₁(𝑡_𝑘),…,𝑢_𝑚(𝑡_𝑘))𝑒_1𝑘+⋯+𝑘_{𝑖𝑗,𝑢𝑚}(𝑡_𝑙,𝑡_𝑘,𝑢₁(𝑡_𝑘),…,𝑢_𝑚(𝑡_𝑘))𝑒_𝑚𝑘+𝑜(ℎ²), (5.4) where we assume that 𝑘_𝑖𝑗(𝑠,𝑡,𝑢) is derivable for 𝑢_𝑙, 𝑙=1,…,𝑚 and let 𝑘_{𝑖𝑗,𝑢𝑙}𝜕(𝑠,𝑡,𝑢)=𝜕𝑢_𝑙𝑘_𝑖𝑗(𝑠,𝑡,𝑢),𝑙=1,…,𝑚. (5.5) Then 𝑒_𝑖0=𝑢_𝑖(𝑡₀)−𝑢_𝑖0𝑒=0,_𝑖𝑙=𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙=^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^{𝛽𝑖𝑗𝑚}_𝑟=1𝑘_{𝑖𝑗,𝑢𝑟}(𝑡_𝑙,𝑡_𝑘,𝑢(𝑡_𝑘))𝑒_𝑟𝑘]+^𝑚_𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗𝑚}_𝑟=1𝑘_{𝑖𝑗,𝑢𝑟}(𝑡_𝑙,𝑡_𝑙,𝑢(𝑡_𝑙))𝑒_𝑟𝑙+ℎ^2+𝛼(lnℎ)^𝛽𝑇⁰_𝑖(𝑡_𝑙)+ℎ^2+𝛼𝑇¹_𝑖(𝑡_𝑙)+𝑂(ℎ^2+𝛼1|lnℎ|^𝛽1),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁. (5.6) Obviously if 𝛽=0, then 𝑇¹_𝑖(𝑡_𝑙)=0.

Now we construct the following auxiliary system of linear Volterra integral equations: find {Φ^𝑟_𝑗(𝑠),𝑗=1,2,…,𝑚,𝑟=0,1} satisfying Φ^𝑟_𝑖(𝑡)=𝑇^𝑟_𝑖(𝑡)+^𝑚_𝑗=1^𝑡₀(𝑡−𝑠)^𝛼𝑖𝑗(ln|𝑡−𝑠|)^𝛽𝑖𝑗[^𝑚_𝑝=1𝑘_{𝑖𝑗,𝑢𝑝}(𝑡,𝑠,𝑢(𝑠))Φ^𝑟₁(𝑠)],𝑖=1,2,…,𝑚,𝑟=0,1, (5.7) and their discrete system of equations: find {Φ^𝑟_𝑖𝑙,𝑖=1,2,…,𝑚,𝑙=1,…,𝑁} satisfying Φ^𝑟_𝑖0Φ=0,_𝑟𝑖𝑙=𝑇^𝑟_𝑖(𝑡_𝑙)+^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗(ln|𝑡_𝑙−𝑡_𝑘|)^{𝛽𝑖𝑗𝑚}_𝑝=1𝑘_{𝑖𝑗,𝑢𝑝}(𝑡_𝑙,𝑡_𝑘,𝑢(𝑡_𝑘))Φ^𝑟_𝑝𝑘]+^𝑚_𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗𝑚}_𝑝=1𝑘_{𝑖𝑗,𝑢𝑝}(𝑡_𝑙,𝑡_𝑙,𝑢(𝑡_𝑘))Φ^𝑟_𝑝𝑘,𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁;𝑟=0,1. (5.8) From Theorem 4.4, we have max_{1≤𝑖≤𝑚,1≤𝑙≤𝑁}|Φ^𝑟_𝑖𝑙−Φ^𝑟_𝑖(𝑡_𝑙)|=𝑂(ℎ^2+𝛼1|lnℎ|^𝛽1),𝑟=0,1. (5.9) Substituting (5.8) and (5.7) into (5.6), we get 𝐸_𝑖0𝐸=0,_𝑖𝑙=^𝑚_𝑗=1[ℎ^𝑙−1_𝑘=0𝑤_𝑙𝑘(𝑡_𝑙−𝑡_𝑘)^𝛼𝑖𝑗𝑘_{𝑖𝑗,𝑢𝑙}(𝑡_𝑙,𝑡_𝑘,𝑢_𝑗(𝑡_𝑘))𝐸_𝑗𝑘]+^𝑚_𝑗=1𝑤_𝑙𝑙ℎ^{1+𝛼𝑖𝑗}𝑘_{𝑖𝑗,𝑢𝑙}(𝑡_𝑙,𝑡_𝑙,𝑢_𝑗(𝑡_𝑙))𝐸_𝑗𝑙+𝑂(ℎ^2+𝛼1(lnℎ)^𝛽1),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁. (5.10) Note that 𝐸_𝑖𝑙=𝑒_𝑖𝑙−ℎ^2+𝛼(lnℎ)^𝛽𝑇⁰_𝑖𝑙−ℎ^2+𝛼𝑇¹_𝑖𝑙. Using Lemma 4.1, we get 𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙−ℎ^2+𝛼(lnℎ)^𝛽Φ⁰_𝑖𝑙(𝑡_𝑙)−ℎ^2+𝛼Φ¹_𝑖𝑙(𝑡_𝑙)=𝑂(ℎ^2+𝛼1(lnℎ)^𝛽1),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁. (5.11) From (5.9) we obtain 𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙−ℎ^2+𝛼(lnℎ)^𝛽Φ⁰_𝑖(𝑡_𝑙)−ℎ^2+𝛼Φ¹_𝑖(𝑡_𝑙)=𝑂(ℎ^2+𝛼1|lnℎ|^𝛽1),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁, (5.12) where Φ¹_𝑖(𝑡)=0 if 𝛽=0. Similarly for Algorithm 3.3 we have 𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙−ℎ^2+𝛼(lnℎ)^𝛽Φ⁰_𝑖𝑙−ℎ^2+𝛼Φ¹_𝑖𝑙=𝑂(ℎ^2+𝛼1|lnℎ|^𝛽1) (5.13) or 𝑢_𝑖(𝑡_𝑙)−𝑢_𝑖𝑙−ℎ^2+𝛼(lnℎ)^𝛽Φ⁰_𝑖(𝑡_𝑙)−ℎ^2+𝛼Φ¹_𝑖(𝑡_𝑙)=𝑂(ℎ^2+𝛼1|lnℎ|^𝛽1),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁, (5.14) where Φ¹_𝑖(𝑡)=0 if 𝛽=0 and Φ¹_𝑖(𝑡_𝑙)=−[2^−𝛼−1ln2𝜁(−𝛼−1)𝜁^′(−𝛼−1)+(2^−𝛼−1−1)]Φ¹_𝑖(𝑡_𝑙),Φ⁰_𝑖(𝑡_𝑙)=(2^−𝛼−1−1)Φ⁰_𝑖(𝑡_𝑙). (5.15) From the above discussion, we have proved the following theorem for the combination method.

Theorem 5.1. (1) If 𝛽=0, then 𝑢^𝑐_𝑖𝑙=(1−2^−1−𝛼)𝑢_𝑖𝑙+𝑢_𝑖𝑙2−2^−1−𝛼=𝑢_𝑖(𝑡_𝑙)+𝑂(ℎ^2+𝛼1),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁. (5.16) Furthermore, one has the following a posteriori error estimate: min{𝑢_𝑖𝑙,𝑢_𝑖𝑙}≤𝑢^𝑐_𝑖𝑙≤max{𝑢_𝑖𝑙,𝑢_𝑖𝑙}or|𝑢_𝑖(𝑡_𝑙𝑢)−^𝑖𝑙+𝑢_𝑖𝑙2𝑢|≤|^𝑖𝑙−𝑢_𝑖𝑙2|. (5.17) That is, one can estimate the average errors |𝑢_𝑖(𝑡_𝑙)−(𝑢_𝑖𝑙+𝑢_𝑖𝑙)/2| by |(𝑢_𝑖𝑙−𝑢_𝑖𝑙)/2|.
(2) If 𝛽=1, then 𝑢^𝑐_𝑖𝑙=(1−2^−1−𝛼)𝑢_𝑖𝑙+𝑢_𝑖𝑙2−2^−1−𝛼=𝑢_𝑖(𝑡_𝑙)+𝑂(ℎ^2+𝛼),𝑖=1,2,…,𝑚,𝑙=1,2,…,𝑁. (5.18) Hence one can obtain a high order of accuracy, which is better than 𝑂(ℎ^2+𝛼(lnℎ)). One also easily deduces that the modified mid-point rectangular quadrature method is better than modified trapezoidal quadrature method when 1−2^−1−𝛼→0.

6. Numerical Examples

In this section, we present two numerical examples to illustrate the features of the combination method discussed in this paper. Let 𝐸_𝑇 denote the error of modified trapezoidal quadrature method, 𝐸_𝑀 the error of modified mid-point rectangular quadrature method, and 𝐸_𝐶 the error of combination method.

Example 6.1. Consider the following system of integral equations with algebraic singularity: 9𝑥(𝑠)=5𝑠^5/3+𝑠−^𝑠₀₃1√𝜋𝑠−𝑡(𝑥(𝑡)+𝑦(𝑡))𝑑𝑡,0≤𝑠≤1,𝑦(𝑠)=2𝑠−16𝑠15^5/2+𝑠+^𝑠₀1√(1𝑠−𝑡2𝑥²1(𝑡)+2𝑦²(𝑡)−𝑦(𝑡))𝑑𝑡,0≤𝑠≤1 (6.1) with the exact solution 𝑥(𝑠)=𝑠, 𝑦(𝑠)=𝑠.
The errors of the approximation solutions obtained by Algorithms 3.2 and 3.3 and their combination are presented in Tables 1, 2, 3, and 4. Numerical results show that the combination method has higher order convergence rate than the two original algorithms. Tables 1–4 also show that the error ratios of Algorithms 3.2 and 3.3 are close to the theoretic value, which is 2^2+(−1/2)=2.828 and the error ratio of the combination method is better than the 2^2+(−1/2).

Table 1 

The errors and the a posterior error estimate of 𝑥(𝑠) in Example 6.1 at mesh point (ℎ=1/10, 𝛼=−1/2 in combination coefficient).

Table 2 

The errors and the a posterior error estimate of 𝑥(𝑠) in Example 6.1 at mesh point (ℎ=1/20, 𝛼=−1/2 in combination coefficient).

Table 3 

The errors and the a posterior error estimate of 𝑦(𝑠) in Example 6.1 at mesh point (ℎ=1/10, 𝛼=−1/2 in combination coefficient).

Table 4 

The errors and the a posterior error estimate of 𝑦(𝑠) in Example 6.1 at mesh point (ℎ=1/20, 𝛼=−1/2 in combination coefficient).

Example 6.2. Consider the following system of integro-differential equations with algebraic singularity: 𝑥^′5(𝑠)=28𝑥(𝑠)^7/5−125𝑠84^12/5+1+^𝑠₀₅1(𝑠−𝑡)³(𝑥²1(𝑡)−10𝑥(𝑡))𝑑𝑡,0≤𝑠≤1,𝑥(0)=0 (6.2) with the exact solution 𝑥(𝑠)=𝑠. Let 5𝑢(𝑠)=𝑥(𝑠),𝑣(𝑠)=𝑥28^7/5(𝑠)−125𝑠84^12/5+1+^𝑠₀₅1(𝑠−𝑡)³(𝑥²1(𝑡)−10𝑥(𝑡))𝑑𝑡. (6.3) Then the original system of integro-differential equations can be transformed to the following system of Volterra integral equations: 𝑢(𝑠)=^𝑠₀5𝑣(𝑡)𝑑𝑡,𝑢(0)=0,𝑣(𝑠)=𝑢28^7/5(𝑠)−125𝑠84^12/5+1+^𝑠₀₅1(𝑠−𝑡)³(𝑢²1(𝑡)−10𝑢(𝑡))𝑑𝑡,𝑣(0)=1. (6.4) The two equations contain algebraic singularity with the coefficient 𝛼=−3/5. The exact solution is 𝑢(𝑠)=𝑠,𝑣(𝑠)=1. (6.5)
The errors of the numerical solutions obtained by Algorithms 3.2 and 3.3 and their combination are presented in Tables 5 and 6. Numerical results still show that the combination method has obviously higher convergence rate than the two algorithms, as well as better maximum error ratio, which is greater than 2^2+(−3/5)=2.639.

Table 5 

The errors and the a posterior error estimate of 𝑥(𝑠), that is, 𝑢(𝑠) in Example 6.2 at mesh point (ℎ=1/10, 𝛼=−3/5 in combination coefficient).

Table 6 

The errors and the a posterior error estimate of 𝑥(𝑠), that is, 𝑢(𝑠) in Example 6.2 at mesh point (ℎ=1/20, 𝛼=−3/5 in combination coefficient).

7. Conclusions

In this paper we use the combination method to solve the systems of integral and integro-differential equations with weakly singular kernels of the second kind, which are important to many applications but have few results. High accuracy and high parallelism are two features of this method. Our numerical results also confirm the theoretical conclusions.

Acknowledgment

This work was supported by the National Natural Science Foundation of China (10671136).