High Order Boundary Conditions Technique for the Computation of Global Homoclinic Bifurcation

In the present paper, a custom algorithm based on the method of orthogonal collocation on finite elements is presented and used for the location of global homoclinic point-to-point asymptotic connecting orbits. This kind of global bifurcation occurs in a large variety of problems in Applied Sciences, being associated to specific, significant physical aspects of the problem under consideration. In order to confront the difficulties faced when the location of such orbits is attempted, high order boundary conditions are constructed through scale order approximations, and used instead of the more common first order ones. The effectiveness of the implemented algorithm is justified by means of the specific applications and the figures presented.


Introduction
In recent years, the improvement of hardware capabilities, such as operational CPU frequency, the increasing amount of RAM equipped, the use of parallelization and the improvement of symbolic mathematical software has enabled researchers to numerically compute global asymptotic orbits more easily, since their computation constitutes a computationally demanding task, even in the case of low dimensional systems.Homoclinic point-to-point connecting orbits arise in various occasions where hysteresis and saturation phenomena are encountered.Also, these orbits act as separatrices for the nonlinear state space in 2D conservative ODEs, since they divide the phase space into regions of periodic

Description of Algorithm
The well-known algorithm of orthogonal collocation on finite elements (widely used in the famous software AUTO86 (see [1] and MATCONT (see [2]) has been implemented and Legendre orthogonal polynomials of maximal degree m , defined by user, have been used for the approximation of the orbits of interest.
Regarding the numerical approximation of the homoclinic orbit of the dynamical differential system ( ) 1 ; , , , : x f x a x a f is transformed to (see [3]) ; , y T T f y α y y τ together with an integral phase condition (see below).Here we have chosen a symmetrically truncated time interval, that is T T T Then, through the appropriate normalization, the independent variable of time is scaled to [0, 1] and the system is further reduced to ( ) Thus after defining the maximal degree of the orthogonal basis polynomials (the Legendre polynomials here), the time interval [ ] , being m internal points in each time subinterval (See Figure 1), are defined as the translated roots of the original Legendre polynomial of degree m .Then, by substituting (4) into (3) we obtain the discretized system of differential equations, that is a system of nonlinear algebraic collocation equations, which are required to be exact at the collocation points.Thus the following ( ) .. , then the requirement that the solution be continuous within the whole time interval, leads to the associated ( ) . Also, the discrete counterparts of the BC that are used for the location of limit cycles are the n equations A technique for the determination of high order BC in the case of homoclinic orbits will be described below.
Finally, the discrete counterpart of an integral type phase condition is utilized for both limit cycles and homoclinic orbits; the continuous form of this condition is (See [1]) where the approximation ( ) ( ) ( ) is used.Note that for the location of limit cycles, by the time normalization T the fundamental period of the limit cycle, the original system (1) is transformed to Figure 1.Meshing of time interval, collocation points.

( ) ( )
; , Also, (5) becomes ( ) So, the period appears explicitly as a system parameter to benefit the numerical continuation performed.Now, by using the Gauss-Legendre quadrature, the discrete counterpart of the integral phase condition (8) takes the form where , i j ω denote the Gauss-Legendre quadrature coefficients-weights, , , , i j i j v y  are the values of the derivatives and the points of the computed orbit itself at a previous step, respectively, and , i j y are the points to be determined.This con- dition has certain advantages over the Poincaré type, scalar ones, especially when continuation of the orbit of interest is carried out.So, counting up the number of the unknowns and the number of the equations for the case of limit cycle location, we deduce that the problem is well-posed, since these numbers are equal.

High Order Boundary Conditions
Both the well-known and widely used techniques of projection BC and the method of eigenvectors [4] are of first order.Thus quite good initial data is required for the successful computation of the orbits of interest and this is becoming harder and harder to achieve as the number of state variables together with the number of active parameters increase.Here, by appropriately combining the multiple scales approximation method with that of successive approximations, we construct a technique for the determination of high order BC, in order to locate numerically homoclinic P2P orbits.The idea for the above mentioned combination comes from Deprit and Henrard [5], Bennett [6] and the relevant references therein.Also, Hassard [7] presented the idea to use high order BC instead of the projection ones.Let us give a quick description of this technique.
We consider a dynamical differential system of the form (1) which possesses a number of equilibrium points, and let ( ) be the fixed point of interest (i.e. the one associated with the homoclinic P2P orbit).Then by setting 0 , 1,..., we expand the right-hand side of (1) in a Taylor series and keep terms up to the fourth order, that is ( ) 1 ,..., ; , 1,..., with i h polynomial functions of 1 ,..., n ξ ξ , of order less or equal to four.
Moreover, assuming ( ), 1,..., i t i n = ξ , can be approximated up to the order of interest in positive integerpowers of a small amplitude orbital parameter, let it be ε , as where k is the desired order of approximation.Then, substituting (14) in the expanded equations (13) and equating the terms of the same order, we derive k respective linear (with regard to the variables ( ) , , 1,..., to the j − order, 1,..., j k = ) scale order systems.Further assuming the fixed point of interest is hyperbolic (i.e.no eigenvalue associated with it, is trivial), the respective systems are solved successively.Then, in order to substitute the solutions associated with a specific order to the higher order systems, all the integration constants which are not associated with the manifold of interest (the high order approximation of which is sought), as well as the integration constants corresponding to the homogeneous part of the respective 2 nd and higher order systems are set equal to zero.Finally, the total solution up to the desired order, k , is given by the n sums of (14).
More precisely, regarding the scale order approximations of the outgoing (incoming) solution vector, associated with the unstable (stable) manifold of the equilibrium (to which the orbit under consideration is homoclinic), by considering the solutions of the respective systems, we remove the homogeneous part of the solution ( 2, 1,..., j j k ≥ = ) and set the integration constants correspond- ing to the eigenvalues with negative (positive) real part equal to zero.Then, we define the high order BC for the computation of the homoclinic orbit, in the form ( ) ( ) where the upper index in 1 y and 1 N y − denotes the i -coordinate of these vector functions, Equation ( 15) are differentiable, so that the standard iterative correcting methods of Jacobian based solvers can be used during the solution procedure.By employing the aforementioned high order BC a computation of the homoclinic orbit of interest is achieved and this type of BC can be proved useful in cases like those mentioned in the introduction.Ideally (that is when (15) is zero for all ), both the start and endpoints are equal to the ones defined through the high order BC presented above and the corresponding relations are differentiable, so that standard iterative correcting, Jacobian based solvers can be used during the solution procedure.All the aforementioned coefficients have been set equal to 1 for the applications under consideration in this work.
In general, considering the system (1) (with is the connecting orbit pair.More specifically, for the a − and ω − limit sets it holds that a( ) , ω( ) where ( ) we have that 0 c c n n ( ) s M x + intersect in at least one point, they intersect in at least one orbit by uniqueness of solutions to the initial value problem and , where , Then we expect the orbit γ to be an isolated connecting orbit if for the tangent spaces at ( ) z t .Furthermore, the connecting orbit is persistent in parametric p − systems as long as the intersection is transversal [3], that is if So, based on (20) and ( 21), the sum of dimensions results in or (since 0 c c n n A more thorough explanation lies in the field of differential topology and more specifically according to transversality theorem of the homonymous theory, when two manifolds of dimensions k and l intersect in an n − dimensional space, then generally their intersection will be a manifold of dimension k l n + − (See [8]).Let us also note that the non-transversal intersections can be perturbed to transversal ones, whereas the transversal ones retain their topology under because the orbit sought is a homoclinic point-to-point one, so according to (23), 1 p = and thus we have chosen one active parameter in each case.

Application to the Lorenz System
The first application of the methodology presented in this paper is on the wellknown Lorenz model [9], which not only serves as a pure mathematical model, but it is applied in various fields of applied physics (thermosyphons, dynamos, meteorology, convection, etc.).This example serves as a validation of the integrity, accuracy and precision of the implemented algorithm.

Equilibrium Analysis
As it is known, the system possesses three equilibria, one at the origin, ( ) 0 0, 0, 0 E and two nonzero ones (see [10])

Hopf Bifurcation-Limit Cycles Continuation
By using the Liu criterion [11], regarding the equilibrium point E + , in order to investigate the occurrence of the Hopf bifurcation, we first evaluate the Jacobian Then by writing the characteristic polynomial  8).(Also, the numerical continuation can be performed by means of the well-known numerical toolbox MATCONT [12]).The Jacobian matrix evaluated at 0 E has two negative eigenvalues, as well as a posi- tive one, namely [6] ( ) The continued limit cycles are presented in Figure 2. As soon as the numerically continued cycles are close enough to the fixed point of interest and the free parameter remains practically unchanged, then the main computation of the homoclinic connecting orbit can be initiated, as the former (the largest cycle) can be considered as a good initial approximation of the latter.

Application of High Order Boundary Conditions
Let us describe the definition and application of high order BC.Assuming the solutions of the dynamical differential system of interest can be approximated by , , where ε denotes the orbital parameter and k is the desired order of approximation.Then, by substituting (29) into (24) and equating the terms of the same order, we get the respective linear (as regards the variables ( , , ) x y z corresponding to the j -order, 1,..., j k = ) systems.In terms of the present analysis the desired order of approximation has been chosen as 4 k = , so that we obtain 1 st order approximation By means of the procedure described in Section 3, we arrive at the approximations of both the outgoing (locally asymptotically unstable) vector solution and the incoming (locally asymptotically stable) one.Then, by using (30), the total solution up to the fourth order is given by x y z , respectively.
Moreover, since all the demanded calculations can be lengthy even for low dimensional systems, the above approximations can be obtained via a symbolic computational package, such as Mathematica or Maplesoft Maple (which offers direct integration with MathworksMatlab).We present below the solutions associated with the outgoing vectors: where 1 c can be user defined and ( Similarly, the expressions of the incoming vector can be set, as well.There, the integration constants associated with the unstable eigenvalues must be set equal to zero.Via the method of orthogonal collocation on finite elements and 4 th or-derBC as described above, the homoclinic connecting orbit of interest has been located inside the truncated time interval [ 7.2018, 7.2018] − , which has been determined by means of the well-known Beyn's method [1].The trajectory of the homoclinic orbit with 13.92655740.. r  is presented in Figure 3 together with the orbit obtained by use of the standard predictor-corrector method Adams-Bashforth-Moulton (ode113 of MathworksMatlab).The improvement achieved by the use of high order BC as compared to the traditional first order ones, is shown in Figure 4, Figure 5

Application to a Model of Three-Species Food Chain with Group Defence Ecosystem
The model used to describe a food-chain with group defense ecosystem is an instantaneous one (i.e.no time delays appear), expressed by the system of autonomous ordinary differential equations (See [13]) with ( ) ( ) ( ) Here ( ) x t denotes the prey population, Journal of Applied Mathematics and Physics   ( ) y t denotes the intermediate population which feeds upon x and ( ) z t de- notes the top predator population that feeds upon y .Additionally, , , ,

K r c d
are positive constants and , , g p q are analytic functions.More specifically, K denotes the carrying capacity of the environment and will be used as the bi- furcation parameter (active parameter).The function ( , ) g x K denotes the spe- cific growth rate of the prey in the absence of predation, ( ) p x denotes the predator functional response and ( ) q y a predator functional response of z on y .In the present analysis logistic growth rate is considered for g , that is , ( ) x p x xe − = and ( ) q y y = .These choices satisfy certain conditions mentioned in [13], which the reader can refer to if further information is sought.Group defence is a phenomenon whereby there is a decrease (or even total prevention) in predation when the numbers of the prey are large, due to the ability of the prey to better defend or disguise themselves.So, the system takes the form

Fixed Point Analysis
The system (42) possesses four types of fixed points.More precisely we have ( The latter is an internal equilibrium (i.e. in general it might or might not be an equilibrium of (42)), where Note that there is a second solution of the system of the third equation of (28) together with (44), be it − .This solution is rejected, as the fixed point associated with it possesses two purely imaginary eigenvalues together with a trivial one, that is the first and second inequality of (28) do not hold.
Additionally, the transversality condition ( ) The sought homoclinic orbit, associated with the equilibrium h E , will constitute a structurally stable phenomenon for the system of interest for 1 p = according to (23).The aforementioned numerically continued limit cycles are presented in Figure 6.

Application of High Order Boundary Conditions
We briefly present the corresponding scale order systems.Thus regarding a system fixed point terms up to the fourth order (hence, finally 4 th order BC are extracted) we have where Then, for ( ) ( ) u v w j = denoting the successive approximations of the state variables, we take the following scale order systems: 1 st order approximation

Conclusion-Discussion
Firstly, an efficient custom algorithm of orthogonal collocation on finite elements implemented in MathworksMatlab has been presented together with two applications in different fields of Applied Mathematics, be them the well-known Lorenz system and a model of a three-species food chain with group defence ecosystem.As a result, global homoclinic asymptotic point-to-point connecting orbits have been obtained numerically, regarding the specific applications.In both cases an initial approximation of the homoclinic connecting orbits of interest has been acquired by continuing limit cycles, which have emerged from a Hopf bifurcation, numerically up to a high value of the fundamental period of them.
The efficiency of the algorithm also lies in the fact that all the required equations (that is, the collocation equations, the continuity equations, the BC and the phase condition) are converted to Matlab functions automatically, so that integrated, sophisticated Matlab routines used for solving systems of nonlinear algebraic equations, as well as optimization routines or any other relevant, userdefined routines can be applied directly for the solution of the aforementioned system of nonlinear algebraic equations.Furthermore, the high order BC defined and used herein can be useful when ordinary PCs of low to moderate computational power are used for the location of homoclinic orbits, as they did not require the increase of the length of the truncation interval in order to achieve the precision sought for the computation.
Finally, regarding the equilibrium point * E of the ecosystem model, the physical meaning of the homoclinic orbit is that if the carrying capacity K is increased (i.e. if enrichment is attempted) above the critical value cr K (leading to a supercritical Hopf bifurcation), then it approaches a limiting value, that of the homoclinic orbit, where the top predator is extinct (i.e. it eventually dies out).So, enrichment needs to be carried out with caution and occasions like that have to be taken seriously into account.Moreover, the homoclinic orbit of the Lorenz system can also be seen as a pure mathematical result as well, while also it serves as validation study of the implemented algorithm and the whole methodology presented in this paper.
state variables and a a parameter of the system, the infinite time horizon ( ), −∞ +∞ is generally truncated to [ ]

 ( 4 )
, the solution of the differential system is approximated by a weighted sum of the basis polynomials Journal of Applied Mathematics and Physics in each time subinterval, where ( ) , i l P τ are the Legendre polynomials of degree l and the coefficients , i l c must be determined.The positions of the colloca- tion points, weighting coefficients (set equal to 1 for the applications under consideration in this paper)and ( , )

pa
∈  ), we should determine the number of control parameters, p , for which the connecting orbits are iso- lated and structurally stable phenomena.The fixed points , x x − + associated with the global connecting orbit are in fact two invariant sets of the system, ,  , respectively, and the orbit { } ( ( ), ) : the connecting orbit is called homoclinic.It is called heteroclinic, otherwise.Assuming that M ± are adequately smooth invariant bifurcation takes place, which is subcritical, since the first Lyapunov coefficient is evaluated ( ) both applications presented in this work has been carried out by means of a custom algorithm of sequential numerical continuation based on the method of orthogonal collocation on finite elements and the integral phase condition (
denote the integration constants (from now on , 1, 2,3, 4,... i c i = will denote integration constants, unless stated otherwise).By setting 2 the next orders of approximation, only the outgoing solution vectors are presented, since the formulae of the respective full solutions are quite lengthy.However, the symbolic engine of Maplesoft Maple proved excellent at P. S. Douris, M. P. Markakis DOI: 10.4236/jamp.

Figure 3 .
Figure 3. Point-to-point homoclinic connecting orbit at
conclude with the occurrence of a supercritical Hopf bifurcation.Thus, a stable limit cycle bifurcates from the equilibrium and it is numerically continued towards the direction of increasing period.

Figure 7 .
Figure 7. Point-to-point homoclinic connecting orbit at a are compact invariant sets of the sys- ±are hyperbolic equilibria of (1) (assumed throughout the present paper unless stated otherwise), that is The above dimension formula can also be expressed in terms of codimension.The codimension of an l -dimensional submanifold of n -space is n l − .Now, if case (22) or equivalently (23) does not hold, then either the two manifolds intersect at more than one orbit, so that extra conditions (restrictions) are necessary for the parameterization of a unique orbit, or the connecting orbit is not a structurally stable phenomenon, so a number of extra parameters need to be considered.For the systems under consideration in this paper, P. S. Douris, M. P. Markakis DOI: 10.4236/jamp.2018.63049560 Journal of Applied Mathematics and Physics perturbation.