Multistage Numerical Picard Iteration Methods for Nonlinear Volterra Integral Equations of the Second Kind

Using the Picard iteration method and treating the involved integration by numerical quadrature formulas, we propose a numerical scheme for the second kind nonlinear Volterra integral equations. For enlarging the convergence region of the Picard iteration method, multistage algorithm is devised. We also introduce an algorithm for problems with some singularities at the limits of integration including fractional integral equations. Numerical tests verify the validity of the proposed schemes.


Introduction
The Volterra integral equations arise in many scientific and engineering fields such as the population dynamics, spread of epidemics, semi-conductor devices, vehicular traffic, the theory of optimal control, the kinetic theory of gases and economics [1]- [7].The initial or boundary value problems for ordinary differential equations and some fractional differential equations can be equivalently expressed by the second-kind Volterra integral equation [6]- [9].
In this work, we consider the general nonlinear Volterra integral equation of the second kind where it permits weak singularity at the limits of integration.The specific conditions under which a solution exists for the nonlinear Volterra integral equation are considered in [1]- [4] [7].Many analytical and numerical methods have been proposed for solving this type of equations, such as the linearization and collocation method [10]- [14], the trapezoidal numerical integration and implicit scheme method [15], the implicit multistep collocation methods [16], the reproducing kernel method [17], the wavelet method [18] [19], the Adomian decomposition method [6] [7] [20] and the methods by using function approximation [21]- [23].
The Picard iteration method, or the successive approximations method, is a direct and convenient technique for the resolution of differential equations.This method solves any problem by finding successive approximations to the solution by starting with the zeroth approximation.The symbolic computation applied to the Picard iteration is considered in [24] [25], and the Picard iteration can be used to generate the Taylor series solution for an ordinary differential equation [25].
In this work, we concern on the numerical Picard iteration methods for nonlinear Volterra integral Equation (1).By using the proposed methods, we treat the involved integrals numerically and enlarge the effective region of convergence of the Picard iteration.The rest of the paper is organized as follows.In Section 2, the scheme in a single interval is considered, and the validity of the method is verified by some numerical tests.Basing on the scheme proposed in Section 2, we devise a multistage algorithm in Section 3 for enlarging the convergence region.In Section 4, an algorithm is introduced for problems with some singularity.To show the effectiveness of the proposed algorithms, we perform some numerical results.

Numerical Picard Iteration Method for Integral Equations
The Picard iteration scheme for the considered Equation ( 1) reads [7] [26] ( ) ( ) t The Picard iteration scheme has been applied in almost each textbook on differential equations to mainly prove the existence and uniqueness of solutions.It is direct and easily learned for numerical calculation.
Assume the recursion scheme is convergent for Treating the integral involved in (4) by numerical quadrature formulas, we have the numerical Picard iteration scheme for (1) over [ ] where 0 n ≥ and , i j ω are the corresponding weights.Considering the compound trapezoidal formula in (6), the weights are 0 0 ,0 , , Numerical results are given to validate the proposed scheme.Let us start with an example in which the integrand ( ) ( ) , , f t s u s is independent with t.Example 1 Consider the initial value problem (IVP) for the nonlinear differential equation ( ) This IVP has the exact solution The equivalent integral equation of the IVP is ( ) ( ) ϕ the result after n iterations when discretization parameter N is taken.Take T = 10, N = 20.
The exact solution can be expressed in terms of the Jacobi elliptic function Integrating the differential equation in (7) yields Similar behavior of errors as in Figure 1 can be observed from Figure 3 which shows  the results of the first 5 iterations and the errors at T for each iteration.It confirms the validity of the scheme ( 5), (6) for equations with general integrand f.What's different from Example 1 is that, at T, the results of the second and the third iterations are even worse than the first one.However, it can be noticed that, in the interval closer to t = 0, for example [ ] 0,3 , the errors decrease as n increases all the same.So the underlying numerical iteration method can be viewed as a point-bypoint correction process.

Multistage Scheme
It's well-known that the convergence of the Picard iteration is constrained in some interval.Then how can we get the numerical solution to the integral Equation (1) when t is outside the interval of convergence?We will take advantage of the multistage method and design a scheme by which the considered problem can be solved interval by interval.For example, the Equation ( 1) is considered on [ ] 0 , t T , however, assume that the single- stage-scheme designed in the previous section is convergent only on [ ] 0 1 , t t , where t 1 < T. For achieving the numerical result at T, we can regard the problem on [ ] 1 , t T as a new one, in which we take the numerical result at 1 t as the initial value.Now we begin to design the multistage scheme in detail.Denote the time interval considered for (1) by [ ] 0 , I t T = .For a given positive integer K, we break Suppose the equation has been solved on 0 , k t t denote the times of iteration by n κ and the iterative solutions by { } . Now we consider the solution on ( ) the right hand side of which will be analyzed below.
• An approximation k N k ϕ  of the first term ( ) k u t has been gotten in previous resolution.• The second part, with the approximations of u on nodes in 0 , k t t     having been gained, can also be ap- proximated ) where the corresponding weights for numerical integration on I κ are ( ) Using ( 5), ( 6) over 1 ,  , numerical solution to (12) can be obtained.We conclude the previous analysis as an algorithm.
Example 3 Consider the Lane-Emden equation ( ) ( ) The equivalent integral form of the Lane-Emden equation is [20] ( ) ( ) First, taking T = 4, 20 N = , we solve the current problem by ( 5), (6).The numerical solutions of the first 5 iterations and the errors at T are shown in Figure 4 from which the convergence can be observed.Unfortunately, the scheme is not convergent for T = 6.
Consider the underlying problem for larger T by Algorithm 1.The time interval [ ] 0,T is uniformly divided into K subintervals, in which the same discretization parameter, denoted by N, is taken.Take 3, 4, 6,12 K = and 10, 20,30, 40 < 10 , where , N K n ϕ denotes the result after n iterations when discretization parameters N and K are taken.Errors and convergence rates respect to N at t = 12 are reported in Table 1, from which one can see that the underlying scheme is of order   In fact, from the errors reported in the table, the convergence order 2 K − can also be obtained.So the scheme is of order ( ) It's an interesting phenomenon observed from Table 1 that almost the same results are obtained for same NK.For example, when NK = 120, the errors are all 9.040e 5 − .This may be because "enough" iteration numbers are taken for all subintervals in the sense of (13).Setting the maximal iteration number allowed for each subinterval to 3 and taking NK = 120, we recalculate the current example up to T = 12 for K = 3, 4, 5, 6, 8, 10, 12, 15.The errors at T are presented in Figure 5(b) which shows the decrement of the errors respect to K.

Problem with Singular Integrand
In recent years, the fractional differential or integral equations are much involved.In fact, fractional integral is a class of integration with weak singular kernel.So many fractional differential and integral equations can be equivalently expressed by the singular Volterra integral equation of the second kind.Let us consider such an integral equation with some singularity.
Example 4 Consider the singular Volterra integral equation [14] ( ) ( ) . Note that in the integrand there has 1 t s − , which is infinity at s t = .In such case, the numerical scheme ( 5), ( 6) and corresponding multistage scheme (Algorithm 1) are not valid any more.
A simple idea is to avoid the value of the integrand at s = t in the numerical integration, so an alternative is to integrate with compound rectangular formula.The only things we need to do are changing the nodes of numerical integration and generating approximations for the values of ϕ on these points since only the values on the nodes { } ) and the corresponding weight by ) ( ) ( ) in which ( ) We present the following algorithm.
Now, we come back to Example 4. Taking 0.5, 0.8 t = to subdivide the time interval [ ] 0,1 and N = 5, 10, 20, 40. Figure 6 presents the dependence of the error on t for each N and that on N at t = 1.The results verify the validity of Algorithm 2 in solving problems with some singularity at the limits of integration.However, the method is of order about only 0.5 N − for this example.Remark 1. Algorithm 2 is devised not especially for singular problems.It's also valid for regular problems.For instance, we recalculate Example 1 with K = 2 and N = 5, 10, 20, 40, 80, 160.Errors and convergence rates respect to N are reported in Table 2, from which we can find the order is 2 N − .

Conclusions
In this work, Picard iteration methods with numerical integration are devised for the second kind nonlinear Volterra integral equations.The Picard iteration method solves the considered nonlinear equation explicitly, while the multistage scheme solves it interval by interval and enlarges the convergence region of the Picard iteration method.Numerical results validate the proposed schemes and algorithms and reveal that the schemes are of order ( )   What should be noticed is that the errors reported in the numerical results decrease exponentially respect to times of iteration n (for example, through simple calculation, we can observe from

Figure 1 ( 10 − 2
Figure 1(a) and Figure 1(b) show the results of the first 5 iterations and the errors at T for each iteration respectively.It's shown in the figure that, the iterative solution converges exponentially respect to iteration number n.The relative errors

( 9 )
leads to a new equation, which is similar to the considered problem (1),

Figure 5 .
Figure 5. Example 3 is simulated by Algorithm 1.(a) Dependence of the error get the initial value of iteration:

Figure 6 .
Figure 6.Example 4 is simulated by Algorithm 2. (a) Dependence of the error

Figure 3 (
b) and Figure 4(b) that the convergence rates are about 4 n − for Examples 2 and 3) and are of order 2 − respect to discretization parameter NK.Future work may concern on enhancing the rate of convergence respect to NK.

Table 1 .
The error −and convergence rate at t = 12 (Example 3 is simulated by Algorithm 1 with various discretization parameter N and number of subintervals K).

Table 2 .
The error