Convergence Analysis of Legendre Pseudospectral Scheme for Solving Nonlinear Systems of Volterra Integral Equations

We are concerned with the extension of a Legendre spectral method to the numerical solution of nonlinear systems of Volterra integral equations of the second kind. It is proved theoretically that the proposed method converges exponentially provided that the solution is sufficiently smooth. Also, three biological systems which are known as the systems of Lotka-Volterra equations are approximately solved by the presented method. Numerical results confirm the theoretical prediction of the exponential rate of convergence.


Introduction
Volterra-type integral equations (VIEs) are the mathematical model of many evolutionary problems with memory arising from biology, chemistry, physics, and engineering.For instance, in several heat transfer problems in physics, the equations are usually replaced by systems of Volterra integral equations (SVIEs).Since just few of these equations (i.e., VIEs and SVIEs) can be solved analytically, it is often necessary to apply appropriate numerical techniques.
The remainder of this paper is organized as follows.The LSM is introduced in Section 2. Convergence analysis of the proposed method is discussed in Section 3. Section 4 states three applications of the desired equation in the biological systems.In Section 5, four types of biological models that are known as Lotka-Volterra system of equations are solved by the LSM to show the efficiency of the presented method and to verify the theoretical results obtained in Section 3. Also some comparisons are made with the existing results.Finally, Section 6 includes some concluding remarks.

Implementation of the Legendre Spectral Method
In this section, we apply the basic idea of Tang et al. [11], which was previously used by the authors in [10,12], for discretizing the nonlinear SVIEs (1).In the procedure of approximation, Legendre Gauss quadrature rule together with Lagrange interpolation is used.
In order to use the spectral method, we consider the collocation points {  }  =0 as the set of (+1) Legendre-Gauss points (i.e., the roots of  +1 () = 0, where  +1 () is the ( + 1)th Legendre polynomial).Assume that the system (1) holds at   : Gauss quadrature rules can be used to compute approximately the integral terms in (2).To this end, we make the change of variable and obtain By applying the ( + 1)-point Gauss quadrature formula associated with the Legendre weights {  }, we get where    = (  ,   ) and the points {  }  =0 coincide with the collocation points {  }  =0 .
Let (  ) =   and V(  ) = V  , for  = 0, 1, . . ., , and represent the Lagrange interpolation polynomials of  and V, in which   is the th Lagrange basis function.The use of these interpolation polynomials for representing (   ) and V(   ) in terms of   and V  implies that Let   and   denote the approximation of   and V  , respectively.Then, from (7), we obtain the discrete problem which is a system of 2( + 1) nonlinear algebraic equations and 2( + 1) unknown coefficients {  }  =0 and {  }  =0 .The nonlinear system (8) can be solved by an appropriate numerical method and the Lagrange interpolation of the solutions can be then obtained from (6).The numerical experiments show that the nonlinear system (8) can be solved easily by fsolve command in the Maple software.

Convergence Analysis
In this section, convergence analysis of the proposed method for the system of Volterra integral equations (1) will be provided.We show that the rate of convergence is exponential.
For convenience, we need some definitions and lemmas for providing the proof of the main Theorem.These lemmas include integral error of the Gauss quadrature rules, estimates of the interpolation error, Lebesgue constant of the Legendre series, and finally the Gronwall inequality.Definition 1.Let  be a bounded interval of R, and let 1 ≤  < +∞.One denotes by   () the space of the measurable functions  :  → R such that ∫   |()|   < +∞.Endowed with the norm it is a Banach space.
Definition 2. Let  be a bounded interval of R, and let  ≥ 0 be an integer.One defines   () to be the vector space of the functions V ∈  2 () such that all the distributional derivatives of  of order up to  can be represented by functions in  2 ().
In short, () is endowed with the inner product for which   () is a Hilbert space.The associated norm is Lemma 3 (integration error from Gauss quadrature [14, page 290]).Assume that a ( + 1)-point Gauss or Gauss-Radau or Gauss-Lobatto quadrature formula relative to the Legendre weight is used to integrate the product , where  ∈   () with  := (−1, 1) for some  ≥ 1 and  ∈ P  .Then, there exists a constant  independent of  such that          ∫ where Lemma 4 (estimates for the interpolation error [14, page 289]).Assume that  ∈   () and denote   () as its interpolation polynomial associated with the ( + 1)-point Gauss or Gauss-Radau or Gauss-Lobatto points {  }  =0 , namely, Then, the following estimates hold: Lemma 5 (Lebesgue constant for the Legendre series [15]).
Assume that   () is the th Lagrange interpolation polynomials associated with the Gauss or Gauss-Radau or Gauss-Lobatto points.Then, where  0 is a bounded constant.

Lemma 6 (Gronwall inequality). If a nonnegative integrable function 𝐸(𝑡) satisfies
where () is an integrable function, then Here, we assume that the kernel (, , ()) has the two following properties which are required for the proof of the convergence analysis: where In the following, we will provide the main theorem of this section in  2 .A similar technique could be designed in  ∞ by using an extrapolation between  2 and  1 [11].

Advances in Mathematical Physics
Proof.According to notation (15), if then the numerical schemes (5) can be written as which gives where From Lemma 3 and the assumption  ≤ 1 (or On the other hand, ( 27) can be rewritten as follows: Multiplying   () on both sides of (30) and summing up from which can be restated in the following form: where   is defined by (23) and the interpolation operator   is defined by (16).It follows from ( 32) and ( 1 Consequently, According to the Lipschitz property of the kernel , we have By letting  = 1 in (17), we have The above estimates, together with (38), yield which leads to (24) provided that  is sufficiently large.This completes the proof.

Applications
The Lotka-Volterra equations model the dynamic behavior of an arbitrary number of competitors [16].Though originally formulated to describe the time history of a biological system, these equations find their application in a number of engineering fields such as nonlinear control.The accurate solutions of the Lotka-Volterra equations may become a difficult task either if the equations are stiff (even with a small number of species) or when the number of species is large [17].Therefore, it is necessary to apply the robust numerical techniques to achieve the best approximations.We refer the interested reader to [18][19][20][21][22] for more information on the biological models and the Lotka-Volterra equations.
First, we consider the prey-predator model: Lotka-Volterra system as an interacting species model to be governed by [23,24] where , , , ,  1 , and  2 are appropriate constants.Here,  1 =  1 () is the prey (e.g., rabbits) population and  2 =  2 () is the predator (e.g., foxes) population at time .As the second system, we consider the simple 2-species Lotka-Volterra competition model with each species  1 and  2 having logistic growth in the absence of the other [23,24]: where  1 , the system (44) would be changed into the following system: where  1 =  1 / 1 and  2 =  2 / 2 are the new initial values.
As the final system, we consider the following version of the Lotka-Volterra equations [23]: where , ,  1 ,  2 , and  3 are constants.These models can be easily transformed into their associated systems of Volterra integral equations.This idea (i.e., changing the IVPs into their associated Volterra integral forms) has been done by many authors like [25,26].This approach of transforming the IVPs into their Volterra form has many interesting advantages such as imposing the initial conditions to the new equations and ease of applying high order Gauss quadrature rules for getting highly accurate approximations.
For instance, the first model ( 43) can be transformed into a system of Volterra integral equations as follows: For the ease of applying the spectral method, as [11], we make the change of variable  = ((1 + )/2) and observe the problem where (50) Similar procedures can be applied to restate the other models ( 46)-( 47) as Volterra integral equations system in [−1, 1].

Numerical Examples
In this section, some numerical examples are considered to illustrate the efficiency and accuracy of the proposed method.
In all examples, we consider {  }  =0 as the Legendre-Gauss points with the corresponding weights where  +1 () is the (+1)th Legendre polynomial.Also, the nonlinear algebraic systems are solved directly by using fsolve command in Maple 13 software with the Digits environment variable assigned to be 30.All calculations are run on a Pentium 4 PC with 2.70 GHz CPU and 4 GB RAM.In order to show the efficiency of the LSM, we compare our results with those of the other methods that were proposed recently in the literature like Bessel collocation method [24] (BCM), HPM [27], and VIM [28].It should be noted that, in all tables, the absolute values of residuals are provided in the uniform nodes  = 0,  = 0.2,  = 0.4,  = 0.6,  = 0.8, and  = 1.
For comparing the LSM and BCM [24], for  = 3, 6, 9, the absolute values of residual functions for the approximate solutions obtained by the LSM and BCM are provided in Table 5.These comparisons are also depicted in Figure 1.Moreover, Figure 2 displays the residual functions  , () =  , (2 − 1) ( = 1, 2) which are obtained by our method for  = 15.From Table 1 and Figure 1, we observe that the presented method is very effective and the obtained results are better than those of the BCM.In [24], one can see that the BCM is better than the ADM [29] and HPM [27]; thus, the current method is more effective than these methods too.From Figure 2, one can conclude that numerical solution with high accuracy is furnished by the presented method.Example 2. We now consider the following problem for model (46): where  =  =  1 =  2 = 1 and  = 0.8.
We apply the presented method to find the approximate solutions of the equivalent system of Volterra integral equations.functions  , ()( = 1, 2) obtained by LSM for  = 5, 10, 15.In order to compare the results of our method for  = 5 with the five-term HPM solutions [27], the residual functions  , () ( = 1, 2) of these methods are depicted in Figures 3 and 4. From Table 2, one can conclude that the LSM provides the numerical solutions with high accuracy.Also, from Figures 3 and 4, we see that the results obtained by the presented method are better than those obtained by the HPM.
Example 3. At this stage, we solve a typical system in model (47) using the LSM with  =  = 0.1,  1 = 0.2,  2 = 0.3, and  3 = 0.5: After transforming this system to the equivalent Volterra system, we apply the LSM for different values of .Numerical values of the residual functions obtained by LSM are provided in Table 3.Moreover, for comparing with other methods, the numerical results of the variational iteration method (VIM) [28], the 4th-order Runge-Kutta method (RK4), and the LSM for  = 3 are given in Table 4. Also, Figure 5 displays the residual functions  , () ( = 1, 2, 3), for  = 10 which are furnished by the LSM.From these tables and figures, we observe that the results of our method are very accurate and they are better than those obtained by the other methods.
Example 4. As the final example, we consider the following problem that is selected from [30]: for the mentioned values of  and report them in Table 5.

Conclusions
This paper deals with the LSM for computing the approximate solution of the systems of nonlinear Volterra integral equations by using the Lagrange interpolations and Gauss quadrature rules.We demonstrated that the errors of the spectral approximations decay exponentially in the nonlinear case.The numerical results obtained for the solutions of the systems of the Lotka-Volterra equations confirm the spectral accuracy of the LSM.In addition, the comparisons of the residual functions obtained by our scheme with those obtained by other methods show that the LSM is more effective than the other methods.

Figure 1 :
Figure 1: Comparison of the residuals at the selected points of Example 1, for  = 3, 6, and 9.

Table 1 :
Comparison of the residual functions of (53) for the   values.

Table 2 :
The residual functions  , () for the selected nodes of Example 2.

Table 4 :
Comparisons of the numerical solutions for Example 3.