Initialization by time integration

The initialization of a continuation of limit cycles is a nontrivial issue. Probably the most often used method is to start from a Hopf point that is detected on a curve of equilibria, see §7.6. Another powerful method is to compute orbits and see whether some of them converge to (stable) periodic orbits. If such an orbit is found (the easiest way is by graphical inspection) then one uses an (approximation to) a point of the orbit and integrates again over a time interval that is somewhat larger than the period of the orbit but smaller that twice the period. The routine initOrbLC.m in the directory LimitCycle then can initialize the continuation of periodic orbits with a user-chosen free parameter, with the period of orbit as the second free parameter. An example is provided in the testrun testselectcycle.m in the directory Testruns:

OPTIONS = [];
hls = adaptx;
OPTIONS=odeset('RelTol',1e-8);
[t,y] = ode45(hls{2},[0 300],[0.3 0.5 -0.1],OPTIONS,1,0.8);

x1 = y(end,:);

[t,y] = ode45(hls{2},[0 10],x1,OPTIONS,1,0.8);

figure
plot(y(:,1),y(:,2))

p=[1;0.8];
ap=[2];

tolerance=1e-2;
[x0,v0]=initOrbLC(@adaptx,t,y,p,ap,20,4,tolerance);
opt=contset;
opt=contset(opt,'MaxNumPoints',50);
%opt=contset(opt,'TSearchOrder',0);
%opt=contset(opt,'Backward',1);
[xlcc,vlcc,slcc,hlcc,flcc]=cont(@limitcycle,x0,v0,opt);

figure
axes
plotcycle(xlcc,vlcc,slcc,[size(xlcc,1) 1 2]);

In this script adaptx is the dynamical system (54) in §7.7 whose odefile adaptx.m is stored in the directory Testruns/TestSystems. It has three state variables and two parameters. Starting from the point $(0.3;0.5;-0.1)$ and with parameter values $(1,0.8)$ it is integrated over the time span $300.$ The endpoint is called $x1$ and from this point another time integration over the shorter span $10$ is performed. It is then checked graphically that this new integration contains a closed circle, see Figure 10 for a plot in a two-dimensional space.

Figure 10: A closed curve in the adaptx - system.
\begin{figure}
\centering
\centerline{\psfig{figure=ex/SelCyc1.eps,height=6cm}}
\end{figure}

The call

[x0,v0]=initOrbLC(@adaptx,t,y,p,ap,20,4,tolerance);
now initializes the continuation of limit cycles. Here $t$ is a column vector whose entries are the time points (between $0$ and $10$) and $y$ is a matrix whose rows contain the coordinates of the points computed along the orbit. Also, $p=(1;0.8)$ is the parameter vector and $ap=[2]$ contains the indices of the free parameters, so in this case the second parameter is free. Also, the number of test functions is $20$ and the number of collocation points in the discretization of the limit cycle is $4.$ Finally, the tolerance is a threshold for accepting an initial part of the computed orbit as an approximation to a limit cycle.

The last command plots the continuation of a branch of limit cycles with $50$ computed limit cycles in the state space, see Figure 11 (the axis labels were added manually).

Figure 11: 50 limit cycles in the adaptx - system.
\begin{figure}
\centering
\centerline{\psfig{figure=ex/SelCyc2.eps,height=6cm}}
\end{figure}

Note: when the lines

%opt=contset(opt,'TSearchOrder',0);
%opt=contset(opt,'Backward',1);
in testselectcycle are made active, then the initialization of the continuation of limit cycles involves a decreasing order of the unit vectors in continuation space instead of an increasing order. This accidentally leads to a continuation in the opposite direction; therefore cds.options.Backward is set to reverse the continuation.