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 and non-periodic solutions, respectively, as it is known. In that sense, we can consider homoclinic orbits as the limit of periodic solutions, that is a periodic orbit with infinite fundamental period, while it remains bounded.
In the present paper, we present an algorithm for the numerical computation of global homoclinic asymptotic point-to-point connecting orbits (homoclinic P2P orbits for short from now on), where the evaluation of high order boundary conditions (BC from now on) is involved. It is well known that the projection method converges exponentially with the increase of the truncation interval, which in turn, however, can increase the computational time (mainly in ordinary PCs with low to moderate CPU power such as an Intel Core i7 870). So, the use of high order BC can be proved useful in cases like that. In Section 2, a brief description of the algorithm is given and in Section 3 the procedure of the extraction of high order BC is presented. The algorithm is applied to the Lorenz system (Section 4), as well as to a three-species food chain model with group defence ecosystem (Section 5), both being three dimensional cases of dynamical systems. The analysis is carried out with the aid of MathworksMatlab and the symbolic engine of Maplesoft Maple, by means of which we also obtain the respective graphs.
2. Description of Algorithm
The well-known algorithm of orthogonal collocation on finite elements (widely used in the famous software AUTO86 (see  and MATCONT (see  ) has been implemented and Legendre orthogonal polynomials of maximal degree , 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
With the vector of state variables and a parameter of the system, the infinite time horizon is generally truncated to , so that by setting , (1) is transformed to (see  )
together with an integral phase condition (see below). Here we have chosen a symmetrically truncated time interval, that is . Then, through the appropriate normalization, the independent variable of time is scaled to 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 is divided into subintervals-elements, . Moreover, the solution of the differential system is approximated by a weighted sum of the basis polynomials
in each time subinterval, where are the Legendre polynomials of degree and the coefficients must be determined. The positions of the collocation points, , being internal points in each time subinterval (See Figure 1), are defined as the translated roots of the original Legendre polynomial of degree . 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 equations ( is the number of state variables) should hold:
for and the scaled independent variable of time. By setting
then the requirement that the solution be continuous within the whole time interval, leads to the associated continuation (matching) conditions
where . Also, the discrete counterparts of the BC that are used for the location of limit cycles are the 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  )
where the approximation is used. Note that for the location of limit cycles, by the time normalization , mapping to , with 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 denote the Gauss-Legendre quadrature coefficients-weights, are the values of the derivatives and the points of the computed orbit itself at a previous step, respectively, and are the points to be determined. This condition 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.
3. High Order Boundary Conditions
Both the well-known and widely used techniques of projection BC and the method of eigenvectors  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  , Bennett  and the relevant references therein. Also, Hassard  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
we expand the right-hand side of (1) in a Taylor series and keep terms up to the fourth order, that is
with polynomial functions of , of order less or equal to four. Moreover, assuming , can be approximated up to the order of interest in positive integerpowers of a small amplitude orbital parameter, let it be , as
where is the desired order of approximation. Then, substituting (14) in the expanded equations (13) and equating the terms of the same order, we derive respective linear (with regard to the variables corresponding to the order, ) 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 2nd and higher order systems are set equal to zero. Finally, the total solution up to the desired order, , is given by the 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 ( ) and set the integration constants corresponding 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 and denotes the -coordinate of these vector functions, are some positive weighting coefficients (set equal to 1 for the applications under consideration in this paper)and have been evaluated via (14).
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 ), we should determine the number of control parameters, , for which the connecting orbits are isolated and structurally stable phenomena. The fixed points associated with the global connecting orbit are in fact two invariant sets of the system, , respectively, and the orbit is called a connecting orbit from to if
where is the connecting orbit pair. More specifically, for the and limit sets it holds that
So, if , the connecting orbit is called homoclinic. It is called heteroclinic, otherwise. Assuming that are adequately smooth invariant manifolds of dimensions , where denote the dimensions of the manifolds in the following decomposition (See  )
where is compact and are compact invariant sets of the system. If we further assume that , have unstable manifolds and stable manifolds with dimensions and , respectively, and also have a hyperbolic structure (so that these are independent of ), then . However, for the steady case, since are hyperbolic equilibria of (1) (assumed throughout the present paper unless stated otherwise), that is
we have that , so . Whenever and intersect in at least one point, they intersect in at least one orbit by uniqueness of solutions to the initial value problem and
Then we expect the orbit to be an isolated connecting orbit if
for all (20)
for the tangent spaces at . Furthermore, the connecting orbit is persistent in systems as long as the intersection is transversal  , that is if
So, based on (20) and (21), the sum of dimensions results in
or (since )
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 and intersect in an dimensional space, then generally their intersection will be a manifold of dimension (See  ). Let us also note that the non-transversal intersections can be perturbed to transversal ones, whereas the transversal ones retain their topology under perturbation. The above dimension formula can also be expressed in terms of codimension. The codimension of an -dimensional submanifold of -space is . 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, , because the orbit sought is a homoclinic point-to-point one, so according to (23), and thus we have chosen one active parameter in each case.
4. Application to the Lorenz System
The first application of the methodology presented in this paper is on the well-known Lorenz model  , 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. The system equations are
We choose the standard parameter values and use as the bifurcation parameter.
4.1. Equilibrium Analysis
As it is known, the system possesses three equilibria, one at the origin, and two nonzero ones (see  )
4.2. Hopf Bifurcation-Limit Cycles Continuation
By using the Liu criterion  , regarding the equilibrium point , in order to investigate the occurrence of the Hopf bifurcation, we first evaluate the Jacobian
Then by writing the characteristic polynomial
where , the above mentioned criterion yields
hence a Hopf bifurcation takes place, which is subcritical, since the first Lyapunov coefficient is evaluated . The limit cycles are continued up to one close enough to the saddle equilibrium . The numerical continuation for 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 (8). (Also, the numerical continuation can be performed by means of the well-known numerical toolbox MATCONT  ). The Jacobian matrix evaluated at has two negative eigenvalues, as well as a positive one, namely 
The continued limit cycles are presented in Figure 2.
Figure 2. Continued limit cycles for the Lorenz system with .
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.
4.3. 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 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 corresponding to the -order, ) systems. In terms of the present analysis the desired order of approximation has been chosen as , so that we obtain
1st order approximation
2nd order approximation
3rd order approximation
4th 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
where represent the outgoing or the incoming solution and stand for or , 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:
1st order approximation
where denote the integration constants (from now on will denote integration constants, unless stated otherwise). By setting we obtain the first order approximation of the outgoing solution vector
2nd order approximation
and setting we get the second order approximation
Regarding 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 handling them.
3rd order approximation
4th order approximation
where 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 4th orderBC as described above, the homoclinic connecting orbit of interest has been located inside the truncated time interval , which has been determined by means of the well-known Beyn’s method  . The trajectory of the homoclinic orbit with 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(a) and Figure 5(b).
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  )
with . Here denotes the prey population,
Figure 3. Point-to-point homoclinic connecting orbit at for the Lorenz system with the method of orthogonal collocation on finite elements and .
Figure 4. Comparison between 1st, 2nd and 4th order BC for the Lorenz system (outgoing vector) with .
Figure 5. Comparison between (a) 1st and 4th order BC and (b) 2nd and 4th order BC, for the Lorenz system (incoming vector) with .
denotes the intermediate population which feeds upon and denotes the top predator population that feeds upon . Additionally, are positive constants and are analytic functions. More specifically, denotes the carrying capacity of the environment and will be used as the bifurcation parameter (active parameter). The function denotes the specific growth rate of the prey in the absence of predation, denotes the predator functional response and a predator functional response of on . In the present analysis logistic growth rate is considered for , that is , and . These choices satisfy certain conditions mentioned in  , 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
5.1. 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
with . Thus regarding the fixed point , by evaluating the Jacobianmatrix , then the third equation of (28) together with (44) yield the critical equilibrium, as well as the critical value of the active parameter . In particular, by setting , we get
Also, the corresponding critical eigenvalues are
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
is verified numerically. Then by computing the first Lyapunov coefficient at , , we 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. The sought homoclinic orbit, associated with the equilibrium , will constitute a structurally stable phenomenon for the system of interest for according to (23). The aforementioned numerically continued limit cycles are presented in Figure 6.
5.2. Application of High Order Boundary Conditions
We briefly present the corresponding scale order systems. Thus regarding a system fixed point , that is or , by setting
expanding the right-hand sides of (42) in Taylor series around and keeping
Figure 6. Numerically continued limit cycles for the model of the three-species food chain with group defence ecosystem, for .
terms up to the fourth order (hence, finally 4th order BC are extracted) we have
Then, for , where are given by (43), by substituting (14) in (48) with , and denoting the successive approximations of the state variables, we take the following scale order systems:
1st order approximation
2nd order approximation
3rd order approximation
4th order approximation
For , the homoclinic orbit of interest is numerically computed within the truncated time interval with the value of the active parameter being equal to . Then, we have and the respective homoclinic orbit is presented in Figure 7.
Firstly, an efficient custom algorithm of orthogonal collocation on finite elements implemented in MathworksMatlab has been presented together with two
Figure 7. Point-to-point homoclinic connecting orbit at for the model of the three-species food chain with group defence ecosystem, for .
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, user-defined 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 of the ecosystem model, the physical meaning of the homoclinic orbit is that if the carrying capacity is increased (i.e. if enrichment is attempted) above the critical value (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.
P.S. Douris is pleased to acknowledge his financial support from “Andreas Mentzelopoulos Scholarships for the University of Patras”.
The authors declare that there is no conflict of interest regarding the publication of this paper.
 De Witte, V., Govaerts, W., Kuznetsov, Y.A. and Friedman, M. (2012) Interactive Initialization and Continuation of Homoclinic and Heteroclinic Orbitsin MATLAB. ACM Transactions on Mathematical Software (TOMS), 38, 18.
 Liu, Y. and Liu, L. and Tang,T. (1994) The Numerical Computation of Connecting Orbits in Dynamical Systems: A Rational Spectral Approach. Journal of Computational Physics, 111, 373-380.