Numerical analysis of the Volterra differential equations via a combination of evolutionary algorithms with moving least squares and ﬁnite difference methods

In this paper, the Volterra differential equations governing on the predator-prey dynamics is numerically solved through the combination of evolutionary algorithms with Moving Least Squares as well as the Finite Difference Method. Among the numerous evolutionary algorithms, the Genetic Algorithms and Particle Swarm Optimization, as two powerful and effective approaches, are selected. The penalty method is chosen to impose the initial conditions on the objective function, while the Moving Least Squares is implemented to interpolate the derivatives in the differential equations. In order to illustrate the effectiveness of the Moving Least Squares, the results are compared with those obtained via the Finite Difference Method based evolutionary algorithms. Moreover, the results are associated with those of a numerical method and an analytical scheme to evaluate the obtained populations of the prey and the predator. These comparisons clearly illustrate the reliability and efﬁciency of the proposed methodology to solve nonlinear differential equations.


INTRODUCTION
Most governing equations of real world phenomena are enormously nonlinear; and it is challenging to find their solutions [1][2][3][4][5]. Specially, Volterra differential equations governed on the predator-prey dynamics belong to these drawbacks that many researchers have been attracted to solve them [6][7][8][9][10][11][12][13][14][15][16]. Some notable studies on Volterra differential equations are as follows. The stability analysis and convergence studies of Volterra integro-differential equations were conducted based on the convolution test equation and theoretical results and using the corresponding recurrence relation and stability matrix [17]. The nonlinear stability properties of the discontinuous Galerkin method were investigated for stiff Volterra functional differential equations [18]. Spectral collocation methods were developed for weakly singular Volterra integro-differential equations, and convergence analysis results illustrated that the global convergence order depends on the regularities of the kernels functions [19]. solutions of the nonlinear Volterra integro-differential equations [20]. The multi-wavelets Galerkin method was used for linear and nonlinear second order Volterra integro-differential equations, and the problem was reduced to a set of algebraic equations using the operational matrices of integration and the wavelet transform matrix [21]. On the other hand, the moving least squares approach is implemented to reconstruct continuous functions from a set of unorganized point samples [22][23][24][25][26]. As another discretization scheme, the finite difference method converts linear/nonlinear differential equations into a system of linear/non-linear algebraic equations [27][28][29][30][31][32][33][34][35].
Besides, evolutionary algorithms inspired by biological evolution are primarily suited for numerical optimization problems because they often found well approximating solutions for all types of challenges. In the recent years, numerous types of these schemes have been introduced by scientists and researchers, to name but a few, genetic algorithm [39], particle swarm optimization [40], ant colony [36], artificial bee colony [37], bat algorithm [38], firefly algorithm [46], grey wolf optimizer [47], invasive weed optimization [48], shuffled frog leaping algorithm [49], spider algorithm [50], whale optimization algorithm [51] and team game algorithm [52].
The motivation of this paper is to extend two evolutionary algorithms, namely the genetic algorithm [39] and particle swarm optimization [40], to solve the predator-prey dynamics. To reach this goal, the penalty method is used to impose the boundary conditions, and the MLS and FDM are implemented to approximate the derivatives appeared in the equation. The success and effectiveness of the introduced methodology are shown via comparing its results with reported ones in literature.
The rest of the paper is organized as follows. The governing equations are briefly represented in Section 2. The moving least squares approach is formulated in Section 3. The PSO and GA are described Sections 4 and 5, respectively. Discretization of the governing equations to provide a proper objective function is explained in Section 6. Simulation results and discussions are given in Section 7. Finally, Section 8 concludes the findings of this paper.

GOVERNING EQUATIONS
There are some rabbits and foxes living together that foxes eat the rabbits and rabbits eat clover, and there is an increase and decrease in their number. The governing relations to this problem known as Volterra differential equation could be stated as follows [41].
where, ( ) and ( ) respectively represent the populations of the rabbits and foxes at time . Further, , , and are constant parameters dependent to the environment conditions.

MOVING LEAST SQUARES
The moving least squares approximation is formulized in this section. At first, the shape function ( i ) for a selected point where, denotes the monomial basis vector, and m is its order. Matrices Γ( ) and E( ) are calculated by the following relations.

( )
= Ω T and E ( ) where, = [ ( 1 ), ( 2 ), … , ( n )] T signifies the monomial basis matrix; n represents to the number of nodes in the supporting domain Ψ of the global domain Ψ; and Λ defines the diagonal weight function matrix whose jth diagonal member is calculated by the following equation.
where, j = 1, 2, … , n;̄i j denotes the distance between points i and j ; and represents the size of the supporting domain. After calculation of the shape function ( i ), the unknown value F h ( i ) is approximated from the known nodal valuesF ( ) as bellow [42].

PARTICLE SWARM OPTIMIZATION
In the particle swarm optimization algorithm, the position of each particle is enhanced with respect to its previous best position and the position of the best particle among population. The mathematical equations of the PSO are generally represented as follows: where, ⃗ i (it ) and ⃗ V i (it ) describe the position and velocity of particle i at certain iteration it , respectively. Moreover, c 1 and c 2 are coefficients that express the tendency of the particle to the personal and social successes, respectively. ⃖⃗ R 1 and ⃖⃗ R 2 are random vectors in interval [0,1], and G namely inertia weight illustrates the effect of the particle velocity at the previous iteration on the next one. Furthermore, ⃗ pbest i (it ) and ⃗ gbest (it ) are personal and global best positions of particle i and the swarm, respectively [43]. In this paper, the tendency coefficients are regarded as c 1 = c 2 = 3; and the inertia weight is calculated through the following adaptive relation [44].
where, r presents the evolutionary factor and could be calculated by the following equation.
where, s i is the mean distance of particle i to all other particles. N and D are the population size and the number of dimensions, respectively. s i for the global best particle is shown as s global . It is clear that by comparing s i , the maximum and minimum distances (s max and s min ) would be determined.

GENETIC ALGORITHM
The genetic algorithm is considered as an approach for solving the optimization problems based on the biological evolution via modification of a population of individual solutions, repeatedly. At each iteration, individuals are chosen randomly from the current population (as parents) and employed to produce the children for the next generation. The main operators of the real coded genetic algorithm utilized in this work could be described as follows.
• Crossover: This operator selects two chromosomes from the mating pool, randomly. Then, it creates two offspring from the selected parents, and the generated offspring are replaced by the patents. Let ⃗ i (it ) and ⃗ j (it ) represent two random selected chromosomes, and ⃗ i (it ) has the smaller fitness value, the crossover formula could be written as follows: where, ⃖⃗ R 1 and ⃖⃗ R 2 ∈ [0, 1] are random vectors. After calculation of Equations (11) and (12), the superior offspring between ⃗ (it ) and ⃗ (it + 1) should be selected.
• Mutation: The mutation operation causes some varieties in the population and decreases the probability of convergence towards local optima. Let ⃗ i (it ) represents a randomly selected chromosome, then the mutation formula used in this algorithm is defined as: where, ⃗ max and ⃗ min are upper and lower bounds related to the search domain, respectively. Moreover, ⃗ R ∈ [0, 1] represents a random vector [39].

DISCRETIZATION OF THE GOVERNING EQUATIONS AS THE OBJECTIVE FUNCTION
In this section, the objective function of the problem is identified for minimization by the detailed PSO algorithm. Suppose that i denotes nodal position i on the problem domain, and n ( i ) and n ( i ) respectively indicate the numerical values of the populations of the rabbits and the foxes at i computed via the regarded approaches. Using Equation (6), Equations (1) and (2) could be rewritten for node i as follows.
that, (0) = 0 and (0) = 0 . The summation of absolute values of these equations for all nodal points generates the first part of the regarded objective function.

NUMERICAL RESULTS AND COMPARISONS
In this section, the introduced strategy based on the PSO/GA and MLS/FDM is applied to analyse the nonlinear predatorprey dynamics in several cases. The initial population size for the PSO and GA is considered as 100; their maximum number of iterations is set as 10,000; the radius of the supporting domain is regarded as = ( max − min ) + 0.1; 6 nodes are regularly distributed on the defined domains; and the penalty parameters are set at Π 1 = Π 2 = 100.
The global domain is considered as a line geometry from min = 0 to max = 1. The convergence trajectories of the MLS based PSO, MLS based GA, FDM based PSO and FDM based GA methods are depicted in Figure 1. For this case, the found numerical solutions by the four mentioned methods are compared together and with those of the Runge-Kutta Method (RKM) [45] and the exact solutions in Figure 2. These comparisons obviously demonstrate that the MLS based PSO and

Example 2
In this case, in order to have a reasonable comparison with the method presented in Ref. [8], the fixed values in the governing equations are regarded similar to this reference as = = = 1 and = 0.1; and the initial conditions are appreciated as 0 = 14 and 0 = 18. In Ref. [8], the problem has been solved through Homotopy Perturbation Method (HPM) and the global domain is considered as a line geometry from min = 0 to max = 0.125. The convergence diagrams of the evolutionary algorithms are depicted in Figure 3. The numerical results obtained by these methods are compared together and with those of the RKM [45] and HPM [8] in Figure 4. This figure shows that MLS based PSO has most agreement with RKM for the prey population, while FDM based PSO has less error for the predator population. These comparisons obviously display the ability and competence of the methods based PSO algorithm to solve the considered nonlinear differential equation.  Figure 5. The explanations obtained by these methods are compared based schemes with those of the RKM [45] and HPM [8] in Figure 6 that clearly verifies the ability of the PSO algorithm to find the answers of the prey and predator population problem.

Example 4
For this example, the constant parameters and the initial conditions are regarded as = = = 1, = 0.1, 0 = 14 and  Figure 7. The solutions obtained by these methods are associated with those of the RKM [45] and HPM [8] in Figure 8. These assessments visibly confirm the effectiveness and efficiency of the MLS based PSO method to solve the differential equation governed to the considered nonlinear problem.

Example 5
In this case, in order to prepare another fair comparison, the parameters in the Volterra equations are considered similar to Ref. [8] as = = = 1 and = 0.1 with the initial conditions as 0 = 14 and 0 = 18. As it was mentioned earlier, in Ref. [8],  Figure 9. The populations of the preys and predators are compared with those of the RKM [45] and HPM [8] in Figure 10. These comparisons demonstrate that the MLS based PSO method could find the nearest solutions to the RKM for the prey population while the FDM based PSO approach has found the closest solutions to the RKM for the predator population.

Example 6
In this instance, according to Ref. [8], the constants, initial conditions and problem domain are fixed at = = = = 1, 0 = 3, 0 = 2, min = 0, max = 1. The convergences of the GA and PSO based methodologies are signified in Figure 11. The numerical findings achieved by these methods are contrasted together and with those of the RKM [45] and HPM [8]  in Figure 12. These comparisons obviously verify the success of the MLS based PSO method to solve the considered Volterra differential equation.

CONCLUSIONS
In this paper, new combinations of the evolutionary algorithms with the approximation methods were described to analyse the nonlinear Volterra differential equation based predator-prey dynamics. The mathematical modelling of the objective functions were displayed by applying both the FDM and MLS approximation to interpolate the derivatives and the penalty method to impose the initial conditions. Simulation results were illuminated for several cases with different values for the constants and initial conditions of the problem to challenge the introduced strategies. The convergence trajectories of the MLS based PSO, MLS based GA, FDM based PSO and FDM based GA methods were plotted for each case. The diagrams of the prey and predator populations evidently expressed that the MLS based PSO method is competent for solving this class of nonlinear differential equations with an acceptable convergent speed and high accuracy without any restrictive assumptions.