ABSTRACT

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.

Keywords:

Homoclinic Connections, Orthogonal Collocation on Finite Elements, Limit Cycles

1. 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 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 [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

$\stackrel{\dot{}}{x}=f\left(x;a\right),x\in {ℝ}^n,a\in ℝ,f:{ℝ}^{n+1}\to {ℝ}^n$ (1)

With $x={\left(x_1,\mathrm{...},x_n\right)}^T$ the vector of state variables and $a$ a parameter of the system, the infinite time horizon $\left(-\infty ,+\infty \right)$ is generally truncated to $\left[T_{-},T_{+}\right]$ , so that by setting $\tau =t/\left(T_{+}-T_{-}\right)$ , (1) is transformed to (see [3] )

$\stackrel{\dot{}}{y}=\left(T_{+}-T_{-}\right)f\left(y;\alpha \right),y=y\left(\tau \right)$ (2)

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

$\stackrel{\dot{}}{y}=2Tf\left(y;\alpha \right)$ (3)

Thus after defining the maximal degree of the orthogonal basis polynomials (the Legendre polynomials here), the time interval $\left[0,1\right]$ is divided into $N-1$ subintervals-elements, $\left[{\tau }_i,{\tau }_{i+1}\right],i=1,\mathrm{..},N-1$ . Moreover, the solution of the differential system is approximated by a weighted sum of the basis polynomials

$y_i\left(\tau \right)={\displaystyle \underset{l=0}{\overset{m}{\sum }}{c}_{i,l}{P}_{i,l}\left(\tau \right)},i=1,\mathrm{..},N-1,y_i\in {ℝ}^n$ (4)

in each time subinterval, where $P_{i,l}\left(\tau \right)$ are the Legendre polynomials of degree $l$ and the coefficients $c_{i,l}$ must be determined. The positions of the collocation points, ${\tau }_{i,j},j=1,\mathrm{...},m$ , 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 $nm\left(N-1\right)$ equations ( $n$ is the number of state variables) should hold:

${\frac{d{y}_i}{d\tau }|}_{{\tau }_{i,j}}=2Tf\left[y_i\left({\tau }_{i,j}\right)\right]$ (5)

for $i=1,\mathrm{..},N-1\text{,}j=1,\mathrm{..},m$ and $\tau$ the scaled independent variable of time. By setting

$y_{i,j}=y_i\left({\tau }_{i,j}\right),i=1,\mathrm{..},N-1\text{,}j=1,\mathrm{..},m$

then the requirement that the solution be continuous within the whole time interval, leads to the associated $n\left(N-2\right)$ continuation (matching) conditions

$y_{i-1,m}=y_{i,0},i=2,\mathrm{..},N-1$ (6)

where $y_{i,0}=y_i\left({\tau }_i\right),i=1,\mathrm{..},N-1$ . Also, the discrete counterparts of the BC that are used for the location of limit cycles are the $n$ equations

$y_{1,0}=y_{N-1,m}$ (7)

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] )

${\displaystyle {\int }_0^1\stackrel{\dot{}}{\hat{y}}{\left(\tau \right)}^{\text{T}}}\left[y\left(\tau \right)-\hat{y}\left(\tau \right)\right]\text{d}\tau =0$ (8)

where the approximation $\stackrel{\dot{}}{\hat{y}}\left(\tau \right)\approx f\left(y\left(\tau \right)\right)$ is used. Note that for the location of limit cycles, by the time normalization $\tau =t/T_p$ , mapping $\left[0,T_p\right]$ to $\left[0,1\right]$ , with $T_p$ the fundamental period of the limit cycle, the original system (1) is transformed to

$\stackrel{\dot{}}{u}=T_{p}f\left(u;\alpha \right),u=u\left(\tau \right)$ (9)

Also, (5) becomes

${\frac{d{y}_i}{d\tau }|}_{{\tau }_{i,j}}=T_{p}f\left[y_i\left({\tau }_{i,j}\right)\right]$ (10)

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

${\displaystyle \underset{i=1}{\overset{N-1}{\sum }}{\displaystyle \underset{j=0}{\overset{m}{\sum }}{\omega }_{i,j}\langle y_{i,j}-{\hat{y}}_{i,j},{\stackrel{\dot{}}{v}}_{i,j}\rangle }}=0$ (11)

where ${\omega }_{i,j}$ denote the Gauss-Legendre quadrature coefficients-weights, ${\stackrel{\dot{}}{v}}_{i,j},{\hat{y}}_{i,j}$ are the values of the derivatives and the points of the computed orbit itself at a previous step, respectively, and $y_{i,j}$ 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 [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 $x^0={\left(x_1^0,x_2^0,\mathrm{..},x_n^0\right)}^T$ be the fixed point of interest (i.e. the one associated with the homoclinic P2P orbit). Then by setting

${\xi }_i=x_i-x_i^0,i=1,\mathrm{...},n$ (12)

we expand the right-hand side of (1) in a Taylor series and keep terms up to the fourth order, that is

${\stackrel{\dot{}}{\xi }}_i=h_i\left({\xi }_1,\mathrm{...},{\xi }_n;a\right),i=1,\mathrm{...},n$ (13)

with $h_i$ polynomial functions of ${\xi }_1,\mathrm{...},{\xi }_n$ , of order less or equal to four. Moreover, assuming ${\xi }_i(t),i=1,\mathrm{...},n$ , can be approximated up to the order of interest in positive integerpowers of a small amplitude orbital parameter, let it be $\epsilon$ , as

${\xi }_i\left(t\right)\approx {\displaystyle \underset{j=1}{\overset{k}{\sum }}{\epsilon }^j{\xi }_{i,j}}\left(t\right)$ (14)

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 $\left({\xi }_{i,j},i=1,\mathrm{...},n\right)$ corresponding to the $j-$ order, $j=1,\mathrm{...},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 ( $j\ge 2,j=1,\mathrm{...},k$ ) 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

$C_{i,0}{\left[y_1^i\left(\tau =0\right)-{\xi }_i^{out}\right]}^2+C_{i,1}{\left[y_{N-1}^i\left(\tau =1\right)-{\xi }_i^{in}\right]}^2=0,i=1,\mathrm{...},n$ (15)

where the upper index in $y_1$ and $y_{N-1}$ denotes the $i$ -coordinate of these vector functions, $(C_{i,0},C_{i,1})$ are some positive weighting coefficients (set equal to 1 for the applications under consideration in this paper)and $({\xi }_i^{out},{\xi }_i^{in})$ 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 $i=1,\mathrm{...},n$ ), 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 $a\in {ℝ}^p$ ), we should determine the number of control parameters, $p$ , for which the connecting orbits are isolated and structurally stable phenomena. The fixed points $x_{-},x_{+}$ associated with the global connecting orbit are in fact two invariant sets of the system, $M_{-},M_{+}\subset {ℝ}^{n+p}$ , respectively, and the orbit $\gamma =\left\{(x(t),a):t\in ℝ\right\}$ is called a connecting orbit from $M_{-}$ to $M_{+}$ if

$dist\left(z\left(t\right),M_{\pm }\right)\to 0$ as $t\to \pm \infty$ (16)

where $z=\left(x,a\right)$ is the connecting orbit pair. More specifically, for the $\text{a}-$ and $\omega -$ limit sets it holds that

$\text{a}(\gamma )\subset M_{-},\text{ω}(\gamma )\subset M_{+}$ (17)

So, if $\text{a}\left(\gamma \right)=\omega \left(\gamma \right)$ , the connecting orbit is called homoclinic. It is called heteroclinic, otherwise. Assuming that $M_{\pm }$ are adequately smooth invariant manifolds of dimensions $n_{\pm c}+p$ , where $n_{\pm c}$ denote the dimensions of the manifolds $M_{\pm }\left(a\right)$ in the following decomposition (See [3] )

$M={\displaystyle \underset{a\in \Pi }{\cup }\left(M\left(a\right)\times \left\{a\right\}\right)}$ (18)

where $\Pi \subset {ℝ}^p$ is compact and $M\left(a\right)$ are compact invariant sets of the system. If we further assume that $M_{\pm }\left(a\right)$ , $a\in \Pi$ have unstable manifolds $M_{\pm u}\left(a\right)$ and stable manifolds $M_{\pm s}\left(a\right)$ with dimensions $n_{\pm c}+n_{\pm u}$ and $n_{\pm c}+n_{\pm s}$ , respectively, and $M_{\pm }\left(a\right)$ also have a hyperbolic structure (so that these are independent of $a\in \Pi$ ), then $n=n_{\pm c}+n_{\pm s}+n_{\pm u}$ . However, for the steady case, since $M_{\pm }\left(a\right)$ are hyperbolic equilibria of (1) (assumed throughout the present paper unless stated otherwise), that is

$M_{+}\left(a\right)=x_{+}\left(a\right),M_{-}\left(a\right)=x_{-}\left(a\right)$ (19)

we have that $n_{+c}=n_{-c}=0$ , so $n_{\pm u}=n-n_{\pm s}$ . Whenever $M_u\left(x_{-}\right)$ and $M_s\left(x_{+}\right)$ intersect in at least one point, they intersect in at least one orbit by uniqueness of solutions to the initial value problem and

$\gamma \subset M_{-u}\cap M_{+s},\text{where}{M}_{-u}={\displaystyle \underset{a}{\cup }\left(M_{-u}\left(a\right)\times \left\{a\right\}\right)},M_{+s}={\displaystyle \underset{a}{\cup }\left(M_{+s}\left(a\right)\times \left\{a\right\}\right)}$

Then we expect the orbit $\gamma$ to be an isolated connecting orbit if

$T_{z\left(t\right)}{M}_{-u}\cap T_{z\left(t\right)}{M}_{+s}=T_{z\left(t\right)}\gamma =span\left\{\stackrel{\dot{}}{z}\left(t\right)\right\}$ for all $t\in ℝ$ (20)

for the tangent spaces at $z\left(t\right)$ . Furthermore, the connecting orbit is persistent in $p-\text{parametric}$ systems as long as the intersection is transversal [3] , that is if

$T_{z\left(t\right)}{M}_{-u}+T_{z\left(t\right)}{M}_{+s}={ℝ}^{n+p}$ (21)

So, based on (20) and (21), the sum of dimensions results in

$p+n_{-u}+n_{-c}+n_{+s}+n_{+c}+p-1=n+p$ (22)

or (since $n_{+c}=n_{-c}=0$ )

$p=n_{+u}-n_{-u}+1.$ (23)

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 perturbation. 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, $n_{+u}=n_{-u}$ , because the orbit sought is a homoclinic point-to-point one, so according to (23), $p=1$ 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 [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. The system equations are

$\begin{array}{l}\stackrel{\dot{}}{x}=-\sigma x+\sigma y\\ \stackrel{\dot{}}{y}=rx-y-xz\\ \stackrel{\dot{}}{z}=-bz+xy\end{array}$ (24)

We choose the standard parameter values $\sigma =10,b=8/3$ and use $r$ as the bifurcation parameter.

4.1. Equilibrium Analysis

As it is known, the system possesses three equilibria, one at the origin, $E_0\left(0,0,0\right)$ and two nonzero ones (see [10] )

$E_{+}\left(\sqrt{b\left(r-1\right)},\sqrt{b\left(r-1\right)},r-1\right),E_{-}\left(-\sqrt{b\left(r-1\right)},-\sqrt{b\left(r-1\right)},r-1\right)$

4.2. 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

$J=\left(\begin{array}{ccc}-\sigma & \sigma & 0\\ r-z& -1& -x\\ y& x& -b\end{array}\right)$ (25)

so that

$J_{E_{+}}=\left(\begin{array}{ccc}-\sigma & \sigma & 0\\ 1& -1& -\sqrt{b\left(r-1\right)}\\ \sqrt{b\left(r-1\right)}& \sqrt{b\left(r-1\right)}& -b\end{array}\right)$ (26)

Then by writing the characteristic polynomial

${\lambda }^3+p_2{\lambda }^2+p_1\lambda +p_0=0$ (27)

where $p_0=2b\sigma (r-1),p_1=b(\sigma +r),p_2=\sigma +b+1$ , the above mentioned criterion yields

$\begin{array}{l}{p}_0>0\to r>1\\ p_1>0\to b\left(\sigma +r\right)>0\\ D_2\triangleq p_1{p}_2-p_0=0\to r=r_{cr}=\frac{2\sigma +\left(\sigma +b+1\right)\sigma }{\sigma -1-b},\end{array}$ (28)

hence a Hopf bifurcation takes place, which is subcritical, since the first Lyapunov coefficient is evaluated $l_1\left(r_{cr}\right)\simeq 2.55\times {10}^{-3}>0$ . The limit cycles are continued up to one close enough to the saddle equilibrium $E_0\left(0,0,0\right)$ . 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 [12] ). The Jacobian matrix evaluated at $E_0$ has two negative eigenvalues, as well as a positive one, namely [6]

${\lambda }_{1,2}=0.5\times \left[-\sigma -1\pm \sqrt{{\left(\sigma -1\right)}^2+4\sigma r}\right],{\lambda }_3=-b$

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.

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

$x\left(t\right)\approx {\displaystyle \underset{i=1}{\overset{k}{\sum }}{\epsilon }^j{x}_j}\left(t\right),y\left(t\right)\approx {\displaystyle \underset{i=1}{\overset{k}{\sum }}{\epsilon }^j{y}_j}\left(t\right),z\left(t\right)\approx {\displaystyle \underset{i=1}{\overset{k}{\sum }}{\epsilon }^j{z}_j}\left(t\right)$ (29)

where $\epsilon$ 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_j,y_j,z_j)$ corresponding to the $j$ -order, $j=1,\mathrm{...},k$ ) systems. In terms of the present analysis the desired order of approximation has been chosen as $k=4$ , so that we obtain

1^st order approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_1=-\sigma x_1+\sigma y_1\\ {\stackrel{\dot{}}{y}}_1=r{x}_1-y_1\\ {\stackrel{\dot{}}{z}}_1=-b{z}_1\end{array}$ (30)

2^nd order approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_2=-\sigma x_2+\sigma y_2\\ {\stackrel{\dot{}}{y}}_2=r{x}_2-y_2-x_1{z}_1\\ {\stackrel{\dot{}}{z}}_2=-b{z}_2+x_1{y}_1\end{array}$ (31)

3^rd order approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_3=-\sigma x_3+\sigma y_3\\ {\stackrel{\dot{}}{y}}_3=r{x}_3-y_3-x_1{z}_2-x_2{z}_1\\ {\stackrel{\dot{}}{z}}_3=-b{z}_3+x_1{y}_2+x_2{y}_1\end{array}$ (32)

4^th order approximation

$\begin{array}{l}{\stackrel{\dot{}}{x}}_4=-\sigma x_4+\sigma y_4\\ {\stackrel{\dot{}}{y}}_4=r{x}_4-y_4-x_1{z}_3-x_2{z}_2-x_3{z}_1\\ {\stackrel{\dot{}}{z}}_4=-b{z}_4+x_1{y}_3+x_2{y}_2+x_3{y}_1\end{array}$ (33)

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

$\begin{array}{l}x\left(t\right)\approx \epsilon x_1\left(t\right)+{\epsilon }^2{x}_2\left(t\right)+{\epsilon }^3{x}_3\left(t\right)+{\epsilon }^4{x}_4\left(t\right)\\ y\left(t\right)\approx \epsilon y_1\left(t\right)+{\epsilon }^2{y}_2\left(t\right)+{\epsilon }^3{y}_3\left(t\right)+{\epsilon }^4{y}_4\left(t\right)\\ z\left(t\right)\approx \epsilon z_1\left(t\right)+{\epsilon }^2{z}_2\left(t\right)+{\epsilon }^3{z}_3\left(t\right)+{\epsilon }^4{z}_4\left(t\right)\end{array}$ (34)

where$(x,y,z)$ represent the outgoing or the incoming solution and $(x_j,y_j,z_j),j=1,2,3,4$ stand for $(x_j^{out},y_j^{out},z_j^{out})$ or $(x_j^{in},y_j^{in},z_j^{in})$ , 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:

1^st order approximation

$\begin{array}{l}{x}_1\left(t\right)=c_1{e}^{\frac{1}{2}\left(-11+\sqrt{81+40r}\right)t}+c_2{e}^{-\frac{1}{2}\left(11+\sqrt{81+40r}\right)t}\\ y_1\left(t\right)=\left(\frac{9}{20}+\frac{1}{20}\sqrt{81+40r}\right)c_1{e}^{\frac{1}{2}\left(-11+\sqrt{81+40r}\right)t}\\ +\left(\frac{9}{20}-\frac{1}{20}\sqrt{81+40r}\right)c_2{e}^{-\frac{1}{2}\left(11+\sqrt{81+40r}\right)t}\\ z_1\left(t\right)=c_3{e}^{-\frac{8}{3}t}\end{array}$ (35)

where $c_1,c_2,c_3$ denote the integration constants (from now on $c_i,i=1,2,3,4,\mathrm{...}$ will denote integration constants, unless stated otherwise). By setting $c_2=c_3=0$ we obtain the first order approximation of the outgoing solution vector

$\begin{array}{l}{x}_1^{out}=c_1{e}^{\frac{1}{2}\left(-11+\sqrt{81+40r}\right)t}\\ y_1^{out}=\left(\frac{9}{20}+\frac{1}{20}\sqrt{81+40r}\right)c_1{e}^{\frac{1}{2}\left(-11+\sqrt{81+40r}\right)t}\\ z_1^{out}=0\end{array}$ (36)

2^nd order approximation

$\begin{array}{l}{x}_2\left(t\right)=c_4{e}^{\frac{1}{2}\left(-11+\sqrt{81+40r}\right)t}+c_5{e}^{-\frac{1}{2}\left(11+\sqrt{81+40r}\right)t}\\ y_2\left(t\right)=\left(\frac{9}{20}+\frac{1}{20}\sqrt{81+40r}\right)c_4{e}^{\frac{1}{2}\left(-11+\sqrt{81+40r}\right)t}\\ +\left(\frac{9}{20}-\frac{1}{20}\sqrt{81+40r}\right)c_5{e}^{-\frac{1}{2}\left(11+\sqrt{81+40r}\right)t}\\ z_2\left(t\right)=\left(\frac{3}{160}\frac{\left(25+3\sqrt{81+40r}\right)\left(9+\sqrt{81+40r}\right)e^{\frac{1}{3}\left(-25+3\sqrt{81+40r}\right)t}}{13+45r}{c}_1^2+c_6\right)e^{-\frac{8}{3}t}\end{array}$ (37)

and setting $c_4=c_5=c_6=0$ we get the second order approximation

$\begin{array}{l}{x}_2^{out}\left(t\right)=0\\ y_2^{out}\left(t\right)=0\\ z_2^{out}\left(t\right)=\frac{3}{160}\frac{\left(25+3\sqrt{81+40r}\right)\left(9+\sqrt{81+40r}\right)}{13+45r}{c}_1{}^2{e}^{\left(-11+\sqrt{81+40r}\right)t}\end{array}$ (38)

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.

3^rd order approximation

$\begin{array}{l}{x}_3^{out}=C_{x,3}{c}_1{}^3{e}^{\frac{3}{2}\left(-11+\sqrt{81+40r}\right)t}\\ y_3^{out}=C_{y,3}{c}_1{}^3{e}^{\frac{3}{2}\left(-11+\sqrt{81+40r}\right)t}\\ z_3^{out}=0\end{array}$ (39)

where

$C_{x,3}=\frac{1}{16}\frac{-720r^2+\left(-609\sqrt{81+40r}-10503\right)r-2262\sqrt{81+40r}-20358}{\left(r-1\right)\left(160r+203\right)\left(13+45r\right)}$

$C_{y,3}=\frac{1}{80}\frac{\left(-540\sqrt{81+40r}-15930\right)r^2+\left(-5898\sqrt{81+40r}-70722\right)r-7917\sqrt{81+40r}-71253}{\left(r-1\right)\left(160r+203\right)\left(13+45r\right)}$

4^th order approximation

$\begin{array}{l}{x}_4^{out}=0\\ y_4^{out}=0\\ z_4^{out}=C_{z,4}{c}_1{}^4{e}^{\left[\frac{2}{3}\left(3\sqrt{81+40r}-29\right)-\frac{8}{3}\right]t}\end{array}$ (40)

where $c_1$ can be user defined and

$C_{z,4}=-\frac{9}{640}\frac{\left(3\sqrt{81+40r}+29\right)\left(120r^2\sqrt{81+40r}+3940r^2+1649r\sqrt{81+40r}+21551r+316\sqrt{81+40r}+27144\right)}{\left(45r-14\right)\left(r-1\right)\left(160r+203\right)\left(45r+13\right)}$

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 orderBC 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 $r\simeq \mathrm{13.92655740..}$ 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 [13] )

$\begin{array}{l}\stackrel{\dot{}}{x}=xg\left(x,K\right)-yp\left(x\right)\\ \stackrel{\dot{}}{y}=y\left[-r+cp\left(x\right)\right]-zq\left(y\right)\\ \stackrel{\dot{}}{z}=z\left[-s+dq\left(y\right)\right]\end{array}$ (41)

with $x\left(0\right)\ge 0,y\left(0\right)\ge 0,z\left(0\right)\ge 0$ . Here $x(t)$ denotes the prey population,

$y(t)$ denotes the intermediate population which feeds upon $x$ and $z(t)$ denotes 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 bifurcation parameter (active parameter). The function $g(x,K)$ denotes the specific 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 $g(x,K)=1-x/K$ , $p(x)=x{e}^{-x}$ 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

$\begin{array}{l}\stackrel{\dot{}}{x}=x-\frac{1}{K}{x}^2-xy{e}^{-x}\\ \stackrel{\dot{}}{y}=-ry+c\text{}xy{e}^{-x}-yz\\ \stackrel{\dot{}}{z}=-sz+dyz\end{array}$ (42)

5.1. Fixed Point Analysis

The system (42) possesses four types of fixed points. More precisely we have

$\begin{array}{l}{E}_0\left(0,0,0\right),E_K\left(K,0,0\right)\\ E_h\left(x_h,y_h,0\right),x_h{e}^{-x_h}=\frac{r}{c},y_h=\frac{c}{r}{x}_h\left(1-\frac{x_h}{K}\right)\left(\text{oneortwo}\right)\\ E^{*}\left(x^{*},y^{*},z^{*}\right)\left(\text{none,oneortwo}\right)\end{array}$ (43)

The latter is an internal equilibrium (i.e. in general it might or might not be an equilibrium of (42)), where

$1-\frac{x^{*}}{K}-\frac{s}{d}{e}^{-x^{*}}=0,y^{*}=\frac{s}{d},z^{*}=\frac{cd}{s}{x}^{*}\left(1-\frac{x^{*}}{K}\right)-r$ (44)

with $x^{*}>0,y^{*}>0,z^{*}>0$ . Thus regarding the fixed point $E^{*}$ , by evaluating the Jacobianmatrix $J_{E^{*}}(x^{*},K)$ , then the third equation of (28) together with (44) yield the critical equilibrium, as well as the critical value of the active parameter $K$ . In particular, by setting $(r,c,s,d)=(0.28,1,0.5,0.38)$ , we get

$\begin{array}{l}{x}_{cr}{}^{*}\simeq 1.0000000,K_{cr}\simeq \mathrm{1.9381794269...}\\ y_{cr}{}^{*}\simeq 1.31579,z_{cr}{}^{*}\simeq 0.0878794\end{array}$ (45)

Also, the corresponding critical eigenvalues are

${\lambda }_{1,2}=\pm 0.209618i,{\lambda }_{3,s}=-0.0318962$ (46)

Note that there is a second solution of the system of the third equation of (28) together with (44), be it $x=K-1$ . 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

${\frac{d\left(p_1{p}_2-p_0\right)}{dK}|}_{K=K_{cr}}\ne 0$

is verified numerically. Then by computing the first Lyapunov coefficient at $K=K_{cr}$ , $l_1\left(K_{cr}\right)=-9.697250<0$ , 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 $E_h$ , will constitute a structurally stable phenomenon for the system of interest for $p=1$ 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 $E\left(x_0,y_0,z_0\right),\left(x_0,y_0\right)\ne (0,0)$ , that is $E_h$ or $E^{*}$ , by setting

$u=x-x_0,v=y-y_0,w=z-z_0$ (47)

expanding the right-hand sides of (42) in Taylor series around $E$ and keeping

terms up to the fourth order (hence, finally 4^th order BC are extracted) we have

$\begin{array}{l}\stackrel{\dot{}}{u}=H_{1}u+H_{2}v+\frac{1}{2}{H}_{11}{u}^2+H_{12}uv+\frac{1}{6}{H}_{111}{u}^3+\frac{1}{2}{H}_{112}{u}^{2}v+\frac{1}{24}{H}_{1111}{u}^4+\frac{1}{6}{H}_{1112}{u}^{3}v\\ \stackrel{\dot{}}{v}=Q_{1}u+Q_{2}v+Q_{3}w+\frac{1}{2}{Q}_{11}{u}^2+Q_{12}uv+Q_{23}vw+\frac{1}{6}{Q}_{111}{u}^3+\frac{1}{2}{Q}_{112}{u}^{2}v\\ +\frac{1}{24}{Q}_{1111}{u}^4+\frac{1}{6}{Q}_{1112}{u}^{3}v\\ \stackrel{\dot{}}{w}=L_{2}v+L_{3}w+L_{23}vw\end{array}$ (48)

where

$H_1=\frac{K-1}{K}{x}_0-\frac{1}{K}{x}_0^2,H_2=-\frac{x_0}{y_0}\left(1-\frac{x_0}{K}\right),H_{11}=2\frac{K-1}{K}-\frac{K+2}{K}{x}_0+\frac{1}{K}{x}_0^2$

$H_{12}=\frac{1}{y_0}\left(-1+\frac{K+1}{K}{x}_0-\frac{1}{K}{x}_0^2\right),H_{111}=-3+\frac{K+3}{K}{x}_0-\frac{1}{K}{x}_0^2$

$H_{11\text{}\text{}2}=\frac{1}{y_0}\left(2-\frac{K+2}{K}{x}_0+\frac{1}{K}{x}_0^2\right),H_{1111}=4-\frac{K+4}{K}{x}_0+\frac{1}{K}{x}_0^2,H_{1112}=\frac{H_{111}}{y_0}$

$Q_1=-c{y}_0{H}_{12},Q_2=-r-z_0+c\frac{x_0}{y_0}\left(1-\frac{x_0}{K}\right),Q_3=-y_0$

$Q_{11}=-c{y}_0{H}_{112},Q_{12}=-c{H}_{12},Q_{23}=-1,Q_{111}=-c{H}_{111},Q_{112}=-c{H}_{112}$

$Q_{1111}=-c{H}_{1111},Q_{1112}=\frac{-c{H}_{111}}{y_0},L_2=d{z}_0,L_3=-s+d{y}_0,L_{23}=d$

Then, for $\left(x_0,y_0,z_0\right)=\left(x_h,y_h,0\right)$ , where $\left(x_h,y_h\right)$ are given by (43), by substituting (14) in (48) with $k=4$ , $\left({\xi }_1,{\xi }_2,{\xi }_3\right)=\left(u,v,w\right)$ and $\left(u_j,v_j,w_j\right),j=1,2,3,4$ denoting the successive approximations of the state variables, we take the following scale order systems:

1^st order approximation

$\begin{array}{l}{\stackrel{\dot{}}{u}}_1=H_1{u}_1+H_2{v}_1\\ {\stackrel{\dot{}}{v}}_1=Q_1{u}_1+Q_2{v}_1+Q_3{w}_1\\ {\stackrel{\dot{}}{w}}_1=L_2{v}_1+L_3{w}_1\end{array}$ (49)

2^nd order approximation

$\begin{array}{l}{\stackrel{\dot{}}{u}}_2=H_1{u}_2+H_2{v}_2+\frac{1}{2}{H}_{11}{u}_1^2+H_{12}{u}_1{v}_1\\ {\stackrel{\dot{}}{v}}_2=Q_1{u}_2+Q_2{v}_2+Q_3{w}_2+\frac{1}{2}{Q}_{11}{u}_1^2+Q_{12}{u}_1{v}_1+Q_{23}{v}_1{w}_1\\ {\stackrel{\dot{}}{w}}_2=L_2{v}_2+L_3{w}_2+L_{23}{v}_1{w}_1\end{array}$ (50)

3^rd order approximation

$\begin{array}{l}{\dot{u}}_3=H_1{u}_3+H_2{v}_3+H_{11}{u}_1{u}_2+H_{12}\left(u_1{v}_2+u_2{v}_1\right)+\frac{1}{6}{H}_{111}{u}_1^3+\frac{1}{2}{H}_{112}{u}_1^2{v}_1\\ {\dot{v}}_3=Q_1{u}_3+Q_2{v}_3+Q_3{w}_3+Q_{11}{u}_1{u}_2+Q_{12}\left(u_1{v}_2+u_2{v}_1\right)\\ +Q_{23}\left(v_1{w}_2+v_2{w}_1\right)+\frac{1}{6}{Q}_{111}{u}_1^3+\frac{1}{2}{Q}_{112}{u}_1^2{v}_1\\ {\dot{w}}_3=L_2{v}_3+L_3{w}_3+L_{23}\left(v_1{w}_2+v_2{w}_1\right)\end{array}$ (51)

4^th order approximation

For $\left(r,c,s,d\right)=\left(0.28,1,0.5,0.38\right)$ , the homoclinic orbit of interest is numerically computed within the truncated time interval $\left[-831.2371,831.2371\right]$ with the value of the active parameter being equal to $K_h\simeq \mathrm{2.080116..}$ . Then, we have $E_h\left(1.93101,0.494357,0\right)$ and the respective homoclinic orbit is presented in Figure 7.

6. 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, 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 $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 $K_{cr}$ (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.

Acknowledgements

P.S. Douris is pleased to acknowledge his financial support from “Andreas Mentzelopoulos Scholarships for the University of Patras”.

Competing Interests

The authors declare that there is no conflict of interest regarding the publication of this paper.

Cite this paper

Douris, P.S. and Markakis, M.P. (2018) High Order Boundary Conditions Technique for the Computation of Global Homoclinic Bifurcation. Journal of Applied Mathematics and Physics, 6, 554-572. https://doi.org/10.4236/jamp.2018.63049

References