A Block by Block Method for Solving System of Volterra Integral Equations with Continuous and Abel Kernels

The aim of the present paper is to introduce a block by block method for solving system of nonlinear Volterra integral equations with continuous kernels and system of Abel integral equations. We prove convergence of the method and show that its convergence order is at least six. To illustrate performance of the method, numerical experiments are presented and they are compared with HPM (Homotopy Perturbation Method) and RBFN (Radial Basis Function Network) method. The given results demonstrate remarkable ability of the proposed method.


Introduction
Consider a system of Volterra integral equations (VIEs) of the form f (x) = g(x) + There are many physical and biological problems which are modeled as VIEs with weakly singular kernels such as reaction-diffusion problems [9], the behavior of viscoelastic materials in mechanics, superfluidity problems, propagation of a flame, soft tissues like mitral valves of the aorta in the human heart (see [6]). A system of Abel-VIEs may arise either from semi-discretization along space of partial differential equations of fractional order, which models for example some anomalous diffusion and sub-diffusion processes [15] or arises from semi-discretization of Volterra-Fredholm integral equations with singular kernels, some of which occurs in the modeling of coding mechanism in the transmission of nervous signals among neurons [7].
There are also many numerical methods for solving system of VIEs, which some of them are briefly described as follows: The Chebyshev collocation method has been presented in [2] to solve a system of linear integral equations in terms of Chebyshev polynomials. This method transforms the integral system to the corresponding matrix equation by use of Chebyshev collocation points which is solved for the Chebyshev coefficients. Golbabai et al. [8] proposed a novel learning strategy for radial basis function networks (RBFN). By adjusting the parameters of the hidden layer, including the RBF centers and widths, the weights of the output layer are adapted by local optimization methods. In [18] continuous approximations to the solution of the system of Volterra integral equations of the first and second kinds are sought by methods using spline functions of degree m and deficiency−(k − 1), i.e. in C m−k . The resulting method is called an (m, k)method. An expansion method in [16] reduces a system of integral equations to the corresponding linear system of ordinary differential equations. After constructing boundary conditions, this system reduces to a system of algebraic equations that can be solved by a direct or iterative method. Maleknejad et al. [13] used operational matrices of piecewise constant orthogonal functions on the interval [0, 1) to solve system of singular VIEs of convolution type. They use the Laplace transform and then found numerical inversion of Laplace transform by operational matrices. In this field some new versions of the Adomian Decomposition Method (ADM), Homotopy Perturbation Method (HPM) and Runge Kutta methods were also introduced (see for example [3,12,20]).
A block by block method is essentially an extrapolation procedure described by Young [19] for the first time. It has the advantages of being self starting and producing a block of values simultaneously. It is also effective for the problems defined in a large interval. In what follows, we propose a new extension of the block by block method for solving system of Volterra integral equations with continuous and Abel kernels. Of course for Abel equations, to note that the results of this paper are valid if f (x) has sufficient continuity, otherwise it would be necessary to develop special starting formula, for this case see [1,4,5].
Our method has the following extra advantages: 3. In the first step of Romberg quadrature rule, the Simpson rule can be used instead of trapezoidal rule, then the order of convergence will be at least O(h 8 ) for a 4 blocks method.
The rest of the paper is organized as follows. In sections 2 and 3, the method will be described for the system of VIEs with continuous and Abel kernels respectively and then a convergence result will be proven. Finally, the method will be illustrated by numerical results in section 4.
2 System of VIEs with continuous kernels Consider a system of Volterra integral equations of the form

Existence and uniqueness of solution
Sufficient conditions for the existence and uniqueness of solution for equation (2.1) are described by the Theorem 1.

The quadrature rule
In this section, we state basic formulation of the Romberg quadrature rule. To recall this, we use the trapezoidal rule for u v k i,j (x, s, f (s))ds (i, j = 1, 2, ..., n) and define , denotes the trapezoidal rule with 2 l subintervals and .., u. Therefore two and three stages Romberg quadrature rules can be respectively written as

Description of the method
Then put x = x m,l in eq. (2.1) and write it in the form .., F n,j,l ] T (j = 0, 1, ..., m − 1) are known, then the first integral in equation (2.5) can be approximated by standard quadrature rules and the second one is estimated by Romberg quadrature rule. Thus a system of equations will be obtained. In what follows, we describe a method with four blocks which has at least six order of convergence. This can be generalized by increasing the number of blocks to obtain a method with higher order of convergence. Thus let p = 4 and u l = l/4 for l = 0, 1, ..., 4. At first, we describe the method for a single VIE, then we generalize it to the system of VIEs. Thus, we put n = 1 in eq. (2.5) and obtain Then we use eq. (2.3) for the second integral and obtain x m,l xm,0 make a difficulty in computing f 1 (x m,l ). For such a case, we use the Lagrange interpolation at the points x m,0 , x m,1 , x m,2 , x m,3 , x m,4 and obtain where L l (x) is the fundamental Lagrange polynomial. When x changes in a large interval or the step size is very small, we need more accurate approximation for the first integral in eq. (2.5) to reduce the effect of accumulated errors in evaluating the unknown function at the points near to the end of interval. Therefore, it is better the three-stage Romberg rule, B i,j (u, v), be used to approximate this integral. To do this, we consider the following cases.

Case(I):
If m is even, then and substituting it in eq. (2.6) yields Therefore, in each case a system of 4 equations is obtained which is solved for the unknowns F 1,m,1 , F 1,m,2 , F 1,m,3 and F 1,m,4 .

The general process
Now, let us assume that f , g and H are vectors with n component (where H j (t, s) = 1, j = 1, 2, ..., n) and K is an n × n matrix and the grid points are defined similar to the case n = 1. Then eq.(2.5) can be written as for m = 0, 1, 2, ..., N − 1 and l = 1, ..., 4. Consequently, at each step, we get a system of 4n equations with the unknowns F 1,m,1 , F 1,m,2 , ..., F 1,m,4 , F 2,m,1 , ..., F 2,m,4 , ..., F n,m,1 , ..., F n,m,4 , which will be linear and nonlinear respectively for linear and nonlinear integral equations. For the linear case it is solved via a direct method but for the nonlinear case, the system may be solved by using an iterative method or by using a suitable software package such as Maple.

Convergence analysis
Theorem 1. The approximation method given by the system (2.11) is convergent and its order of convergence is at least 6 for the functions k i,j and f i (i, j = 1, 2, ..., n) with at least six-order continuous derivations.
Proof. Define e m,l := F m,l − f (x m,l ), where f (x m,l ) ∈ R n and F m,l ∈ R n denote respectively the exact and approximate solutions of (1.1) at the point x = x m,l . From (2.8) and (2.9) we have (2.12) Then subtracting (2.5) from (2.7) and (2.12) (for l = 1) gives where P(x, s) = [P 1 (x, s), ..., P n (x, s)]. By adding and diminishing the terms and using the conditions of Theorem 1, one can write where R 1 and R 2 are the errors of numerical integrations. Since Q(t, s) is continuous, one can write |Q(x m,1 , x i,j )| ≤ C ij , thus 3 System of Abel VIEs Consider a system of nonlinear Volterra integral equations of the Abel types, i.e.

Existence and uniqueness of solution
Theorem 2. Let: is continuous (i.e. every component is continuous), (ii) K(x,s,f(s)) satisfies Lipshitz condition, i.e. for fixed x, s with 0 ≤ s ≤ x ≤ X there is a positive constant l independent of s and x, such that then (3.1) has a unique solution.
Proof. See [17]. The method of the previous section can't be used for this system, because of the singularity of K i,j (x m,l , s, f (s))H j (x m,l , s) at s = x m,l . So we use the following product integration method instead of a usual integration method. We write and use the Lagrange interpolation, to write

Description of the method
.

Convergence analysis
Theorem 3. The block by block method for Abel type integral equations (introduced in previous subsection) is convergent and its order of convergence is at least 5 for the functions k i,j and f i , (i, j = 1, 2, ..., n) with at least 6 order continuous derivatives.
Let again e m,1 = max l=1,2,3,4 e m,l , then  Table 1 show the superiority of the block by block method in comparing its results (for h = 1/5) with the results of RBFN-MshA (a modified version of Shi's algorithm) [8] and HPM [20], where the results of RBFN-MshA obtained with 6 hidden nodes and the results of HPM obtained with 4 iterations. Moreover, • The time of computation in the block by block method is less than that in the HPM, whenever programming of both methods is done using Maple package. Also, according to the structure of HPM, increasing the number of iterations do not affect on the precision of errors.
• The values of the RBF widths affect significantly on the accuracy of results and determination of them is still a challenging problem, whereas the block by block method dose not need any starting values.
• At each step of the RBFN method, the weights are updated by using an optimization method, but the block by block method is independent of using any other method.Hence the RBFN method is more complicated than the block by block method. Tables 2 and 3 show the numerical results for the Abel systems. Most of the available methods for solving VIEs are based on expansion of solution, for example the Taylor and Chebyshev expansion methods, Tau method, Adomian and homotopy methods and so on. These methods are efficient only for the intervals with small length (say [0, 1] or [-1, 1]) and they are useless for the large intervals, whereas the block by block method is one of the most suitable methods for the large intervals (see table 3 for X = 4).