next up previous contents
Next: Branch switching for a Up: HomCont Demo : Homoclinic Previous: HomCont Demo : Homoclinic   Contents


Branch switching at an inclination flip in Sandstede's model.

Consider the system [#!Sa:95b!#]

\begin{displaymath}\begin{array}{rcl} \dot{x} & = & a x + b y - a x^2 - \alpha z...
...= & c z + \mu x + 3 x z + \alpha (x^2 (1-x) - y^2). \end{array}\end{displaymath} (27.1)

as given in the file sib.c, where for simplicity we have set $ \tilde\mu=0$, $ \beta=1$ and $ \gamma=3$.

We study an inclination flip that exists for $ a=0.375$, $ b=0.625$ and $ c=-0.75$. This corresponds to the situation where the eigenvalues of the equilibrium at the origin are $ a+b=1$, $ a-b=-0.25$ and $ c=-0.75$. Hence, the corresponding bifurcation diagram consists of a complicated structure involving a fan of infinitely many $ n$-periodic and $ n$-homoclinic orbits for arbitrary $ n$ and a region with horseshoe dynamics; see also HoKr:00 and the references therein.

This computation starts from an equilibrium at $ (2/3,0,0)$, which exists for $ a=\mu=\alpha=0$. Also, $ b$ is set to $ 0.625$ (the value we would like it to be) and $ c$ is set to $ -2.5$ in stpnt. Choosing $ c=-2$ at this stage leads to convergence problems. This equilibrium is not the one corresponding to the homoclinic orbit, but it is an equilibrium with complex eigenvalues, that we can follow until it reaches a Hopf bifurcation. A periodic orbit emanates from this Hopf bifurcation and can be followed to the homoclinic orbit. However, first we need to change $ a$ from 0 to $ 0.375$.

All the following commands, except for demo('sib') are contained within the file 'sib.auto' which you can either execute in a batch mode by entering
> auto sib.auto
or step by step using
AUTO> demofile('sib.auto').

We start by copying the demo to the current work directory and running the first step

demo('sib')
ld('sib')
rn()
sv('1')
The equilibrium is followed in $ a$ until $ a$ (or PAR(1)) is at our desired value, $ 0.375$.
  BR    PT  TY  LAB    PAR(1)     ...    U(1)          U(2)          U(3)     
   1     1  EP    1   0.00000E+00 ... 6.66667E-01   0.00000E+00   0.00000E+00
   1     5  EP    2   3.75000E-01 ... 6.66667E-01  -1.33333E-01   0.00000E+00
The output is saved in the files b.1, s.1 and d.1. Next we continue in $ \alpha$ (PAR(4)) until a Hopf bifurcation is found:
rn(c='sib.2',s='1')
sv('2')
or, alternatively,
cc("IRS",2)
cc("ICP",[4])
rn(s='1')
sv('2')
  BR    PT  TY  LAB    PAR(4)     ...    U(1)          U(2)          U(3)     
   1     6  HB    3   3.18429E-01 ... 6.54375E-01  -1.34754E-01   7.70102E-02
The output is saved in the files b.2, s.2 and d.2. This Hopf bifurcation can then be continued into a periodic orbit. The periodic orbit eventually reaches a homoclinic bifurcation. We continue in $ \mu$=PAR(5) and PAR(11), which corresponds to the period, and stop when the period is equal to $ 35$.
rn(c='sib.3',s='2')
sv('3')
  BR    PT  TY  LAB    PAR(5)        L2-NORM    ...    PERIOD    
   3    10        5  -2.41881E-03   6.70569E-01 ...  1.08975E+01
                                                ...
   3    40        8  -1.29495E-02   6.14547E-01 ...  1.41297E+01
                                                ...
   3    81  EP   13  -1.04657E-04   4.01829E-01 ...  3.50000E+01
The output is saved in the files b.3, s.3 and d.3. Note that $ \mu$ first decreases and then increases towards 0, which is precisely what we expect in this model, as homoclinic orbits occur on the line $ \mu=0$ in the $ (\alpha,\mu)$-plane. It is now instructive to look at a phase space diagram to see what is going on.
plot('3')
Selecting 'solution' for Type, [5,6,7,8,9,10,11,12,13] for Label, [0] for X and [1] for Y, we obtain the diagram depicted in Figure 27.1(a), where the periodic orbit grows from the Hopf equilibrium to a homoclinic orbit.
Figure 27.1: Periodic orbit growing from a Hopf bifurcation to a homoclinic orbit (a). The unshifted homoclinic orbit (b).

\begin{picture}(400,180)
\put(0,0){
\put(157,148){(a)}
\includegraphics[scale=0....
...157,148){(b)}
\includegraphics[scale=0.5]{include/notshifted.eps}}
\end{picture}

Note however, that the homoclinic orbit has the wrong left-hand and right-hand end points. This can be seen by plotting the solution corresponding to Label [13] using 't' vs. 'x' (coordinate [0]), as depicted in Figure 27.1(b).

Hence, in order to continue this as a real homoclinic we have to give HomCont special instructions, to do a phase-shift in time. This can be done by setting ISTART=4. Moreover, since we have not specified the value of the equilibrium at the origin in sib.c, we need to set IEQUIB=1 to let HomCont detect the equilibrium. Note that in this case this is not strictly necessary; however, we do this for instructional purposes.

Now we use HomCont to continue the homoclinic orbit in $ c$ and $ \mu$ (PAR(3), PAR(5)), to get the desired value $ c=-2.0$.

rn(c='sib.4',h='sib.shift',s='3')
sv('4')
  BR    PT  TY  LAB    PAR(3)        L2-NORM      ...    PAR(5)     
   3    51  EP   14  -2.00000E+00   4.01890E-01   ...  2.66146E-09
The output is saved in the files b.4, s.4 and d.4. Note that PAR(5)=$ \mu$ remains zero, which is exactly what we expect.

Next we want to add a solution to the adjoint equation to this solution. This is achieved by making the change ITWIST = 1 saved in h.sib.twist. Also, we set ISTART to 1 to tell HomCont that it is should not try to shift the orbit anymore.

rn(c='sib.5',h='sib.twist',s='4')
sv('5')
or, alternatively,
cc("IRS",14)
cc("ICP",[5,8])
cc("NMX",2)
hch("ITWIST",1)
hch("ISTART",1)
rn(s='4')
sv('5')
where hch means ``change HomCont constant''. The output is stored in b.5, s.5 and d.5.
  BR    PT  TY  LAB    PAR(5)        L2-NORM    ...    PAR(8)     
   3     2  EP   15   2.66146E-09   4.01890E-01 ...   1.00000E-02
Here PAR(8) is a dummy (unused) parameter and $ \mu$ just stays where it is. Now that we have obtained the solution of the adjoint equation, we are able to detect inclination flips. This can be achieved by setting NPSI to 1, IPSI(1) to 13, and monitoring PAR(33).
rn(c='sib.6',h='sib.if',s='5')
sv('6')
  BR    PT  TY  LAB    PAR(4)        L2-NORM    ...  PAR(5)        PAR(33)
   3    19  UZ   16   7.11774E-02   4.01890E-01 ... 1.24376E-11  -2.36702E-07
The output is stored in b.6, s.6 and d.6. Hence an inclination flip was found at $ \alpha=0.711774$.

Now we are ready to perform homoclinic branch switching, using the techniques described in [#!OlChKr:03!#]. Our first aim is to find a 2-homoclinic orbit. The ingredients we need are: a homoclinic orbit where $ n$-homoclinic orbits are close by, and the solution to the adjoint equation to obtain the Lin vector. Since both ingredients are there, we can now continue in $ \mu$, $ \varepsilon _1$ and $ T_1$, to obtain the initial Lin gap. Recall from Chapter 20 that the Lin gaps $ \varepsilon_i$ correspond to PAR(20+i*2) and the time intervals $ T_i$ correspond to PAR(21+i*2). We stop when $ \varepsilon _1=0.2$. We need to specify ITWIST=2, to tell HomCont we aim to find a 2-homoclinic orbit, so that it will split it up in three parts with two potential Lin gaps. We effectively have a 9-dimensional system at this point.

rn(c='sib.7',h='sib.hbs2',s='6')
sv('7')
  BR    PT  TY  LAB    PAR(21)       L2-NORM    ...  PAR(22)       PAR(5)     
   3    10       18   3.45897E+01   4.46818E-01 ... 7.87712E-07  -1.55885E-11
   3    20       19   2.73699E+01   4.46818E-01 ... 2.91119E-05  -1.63974E-09
   3    30       20   1.73720E+01   4.46817E-01 ... 4.42273E-03  -3.10167E-05
   3    38  EP   21   1.01451E+01   4.46796E-01 ... 2.00000E-01  -1.48615E-02
The output is stored in b.7, s.7 and d.7. Here we see that $ T_1$, the time it takes to make the first loop with respect to the Poincaré section, decreases. This is illustrated in Figure 27.2. Next we are ready to close this gap, by continuing in $ \alpha$, $ \mu$, and $ \varepsilon _1$, while keeping $ T_1$ at a constant value.
Figure 27.2: Behaviour of the second piece of the `broken homoclinic orbit' when creating a Lin gap (a). Projection of the ``broken homoclinic orbit'' onto the $ (x,y)$-plane, where $ \varepsilon _1=0.2$. To include all the pieces necessary to obtain this figure, the ``X'' box must contain [0,3,6] and the ``Y'' box must contain [1,4,7] (b).

\begin{picture}(400,180)
\put(0,0){
\put(157,148){(a)}
\includegraphics[scale=0....
...put(157,148){(b)}
\includegraphics[scale=0.5]{include/broken.eps}}
\end{picture}
rn(c='sib.8',h='sib.hbs2',s='7')
ap('6')
  BR    PT  TY  LAB    PAR(4)        L2-NORM    ...   PAR(5)        PAR(22)    
   3     3  UZ   22   7.40000E-02   4.46781E-01 ... -1.43162E-02   1.93746E-01
   3    32  EP   23   1.98414E-01   4.46590E-01 ... -6.05495E-03   2.29300E-06
The output is appended to b.6, s.6 and d.6. Now we have obtained a 2-homoclinic orbit at label 24. However, the homoclinic orbit is still split in three parts. We can switch back to a normal orbit by setting ITWIST back to 0 and continuing in the usual way. Here we continue back to the inclination flip point in $ \alpha$ and $ \mu$.
rn(c='sib.9',h='sib.hom',s='6')
ap('6')
  BR    PT  TY  LAB    PAR(4)        L2-NORM    ...   PAR(5)     
   3     7  UZ   24   1.50000E-01   4.94490E-01 ... -3.60248E-03
   3    30  EP   25   7.61403E-02   4.98746E-01 ... -2.64847E-06
So the 2-homoclinic orbit converges back to the 1-homoclinic orbit at the inclination flip bifurcation. The output is appended to b.6, s.6 and d.6. The resulting 2-homoclinic orbits can be seen using
plot('6')
and is depicted in Figure 27.3(a).
Figure 27.3: The 2-homoclinic orbit as $ a$ is changed (a). The two different 3-homoclinic orbits (b).

\begin{picture}(400,180)
\put(0,0){
\put(157,148){(a)}
\includegraphics[scale=0....
...
\put(157,148){(b)}
\includegraphics[scale=0.5]{include/hom3.eps}}
\end{picture}

Next, we aim to find a 3-homoclinic orbit. To do so, we restart at the inclination flip point at label 16 and set ITWIST=3. Moreover, we need to continue in one more gap, $ \varepsilon_2$=PAR(24) and, once again, stop when $ \varepsilon _1$=PAR(22)=0. Note that the dimension of the boundary value problem we continue is now equal to 12. This is not to be confused with the setting of NDIM=3 in the parameter file, because HomCont handles this internally.

rn(c='sib.10',h='sib.hbs3',s='6')
sv('10')
  BR    PT  TY  LAB    PAR(21)    ...   PAR(22)       PAR(24)       PAR(5)     
   3    10       26   3.45896E+01 ...  7.87894E-07   6.42157E-07  -1.06346E-11
   3    20       27   2.73699E+01 ...  2.91126E-05   6.51591E-07  -1.63655E-09
   3    30       28   1.73719E+01 ...  4.42289E-03   1.44090E-04  -3.10188E-05
   3    38  EP   29   1.01451E+01 ...  2.00000E-01   6.97445E-02  -1.48615E-02
The output is stored in b.10, s.10 and d.10. Now we need to subsequently close the Lin gaps. Our strategy is to keep $ T_1$ fixed. We first continue in $ \alpha$, $ \mu$, $ \varepsilon _1$ and $ \varepsilon_2$ until $ \varepsilon_1=0$.
rn(c='sib.11',h='sib.hbs3',s='10')
ap('6')
  BR    PT  TY  LAB    PAR(4)     ...   PAR(5)        PAR(22)       PAR(24)
   3     6  UZ   30   8.20000E-02 ... -1.29790E-02   1.76995E-01   6.37184E-02
   3    32  EP   31   1.98414E-01 ... -6.05495E-03   2.30717E-06   3.62449E-02
The output is appended to b.6, s.6 and d.6. Note that this continuation is very similar to the one where we found a 2-homoclinic orbit. In fact we have now found a 2-homoclinic orbit (numerically) followed by a `broken' 1-homoclinic orbit; only the mesh is not aligned.

The next step is to close the gap corresponding to $ \varepsilon_2$ to obtain a 3-homoclinic orbit. We replace the continuation parameter $ \varepsilon _1$ by $ T_2$, because $ T_2$ (PAR(23)) still has to be decreased from its high value (35) and $ \varepsilon _1$ needs to stay at 0.

rn(c='sib.12',h='sib.hbs3',s='6')
ap('6')
  BR    PT  TY  LAB    PAR(4)     ...   PAR(5)        PAR(23)       PAR(24)
   3    16  UZ   32   1.98395E-01 ... -6.05536E-03   2.01311E+01   1.82491E-08
   3    24  UZ   33   1.80000E-01 ... -6.50293E-03   1.27554E+01  -3.14294E-02
   3    30  UZ   34   1.66990E-01 ... -6.89269E-03   9.41745E+00  -1.03179E-06
   3    32  EP   35   1.78172E-01 ... -6.55364E-03   9.50300E+00  -7.20367E-02
The output is appended to b.6, s.6 and d.6. Note that we have found two zeros of PAR(24), at labels 32 and 34, respectively. The two zeros correspond to two different 3-homoclinic orbits, which, when followed from periodic orbits, both emanate from from the same saddle-node bifurcation. These two 3-homoclinic orbits are depicted in Figure 27.3(b). We can follow both of these back to the inclination flip point, by setting ITWIST back to 0:
rn(c='sib.13',h='sib.hom',s='6')
ap('6')
  BR    PT  TY  LAB    PAR(4)        L2-NORM    ...   PAR(5)     
   3    13  UZ   36   1.29999E-01   5.04807E-01 ... -2.33902E-03
   3    30  EP   37   9.27258E-02   5.06560E-01 ... -2.76788E-04
rn(c='sib.14',h='sib.hom',s='6')
ap('6')
  BR    PT  TY  LAB    PAR(4)        L2-NORM    ...   PAR(5)     
   3     7  UZ   38   1.45000E-01   5.47347E-01 ... -4.79400E-03
   3    30  EP   39   8.39401E-02   5.52605E-01 ... -7.36741E-05
All the output is appended to b.6, s.6 and d.6. The bifurcation diagram and the paths we followed when closing the Lin gaps are depicted in Figure 27.4. It is possible and straightforward to obtain $ 4, 5, 6, \dots$-homoclinic orbits by extending the above strategy.
Figure 27.4: Parameter space diagram near an inclination flip. The curve through label 17 corresponds to a 1-homoclinic orbit. The opening of the Lin gaps occurs along the vertical line from label 16 to label 23. The curves through labels 23 and 30 denote the path that is followed when closing the Lin gaps. The (approximately overlaid) curves though labels 25 and 35 correspond to the 2- and one of the 3-homoclinic orbits. Finally, the curve through label 37 corresponds to the other 3-homoclinic orbit, which was obtained for PAR(22)= $ T_2=12.03201$.
\includegraphics[scale=0.5]{include/parspace.eps}


next up previous contents
Next: Branch switching for a Up: HomCont Demo : Homoclinic Previous: HomCont Demo : Homoclinic   Contents
Gabriel Lord 2007-11-19