Example: an autonomous electronic circuit

We consider the following model of an autonomous electronic circuit [18] where $x,y$ and $z$ are state variables and $\beta,\gamma,r,a_3,b_3,\nu$ are parameters :


\begin{displaymath}
\left\{\begin{array}{rcl}
\dot x&=&(-(\beta+\nu)x+\beta y-...
...ta+\gamma)y-z-b_3(y-x)^3\\
\dot z&=&y
\end{array}\right.
\end{displaymath} (80)

A torus bifurcation in this system is described in the manual of [14]. It is found by starting an equilibrium curve from the trival equilibrium point ($x=y=z=0$) at $\beta=0.5$, $\gamma=-0.6$, $r=0.6$, $a_3=0.32858$, $b_3=0.93358$, $\nu=-0.9$. The free parameter is $\nu$ and the branch is the trivial one ($x=y=z=0$). On this branch a Hopf bifurcation is detected at $\nu=-0.58934$. On the emerging branch of limit cycles a branch point of limit cycles is found; by continuing the newly found branch one detects a torus bifurcation of periodic orbits.

We proceed in a somewhat different way; to avoid the branch point of periodic orbits we start with a slightly perturbed system where the last equation of (80) is replaced by

\begin{displaymath}
\dot z=y+\epsilon
\end{displaymath}

where $\epsilon$ is a new (artificial) parameter with initially $\epsilon=0.001$. The perturbed system is introduced in MATCONT in the odefile torBPC.m in the directory Testruns/TestSystems. It contains a user function $\epsilon$ so that the zeroes of the user function correspond to solutions of the unperturbed system (80).

The trivial solution is replaced by the equilibrium solution $x=0.00125$, $y=-0.001$, $z=0.00052502$ of the perturbed system torBPC.m and we compute a branch of equilibria with free parameter $\nu$ in the now standard way (note that $\nu$ is the sixth parameter, $\epsilon$ is the seventh):

>> p=[0.5;-0.6;0.6;0.32858;0.93358;-0.9;0.001];
>> [x0,v0]=init_EP_EP(@torBPC,[0.00125;-0.001;0.00052502],p,[6]);
>> opt=contset; opt= contset(opt,'Singularities',1);
>> opt=contset(opt,'MaxNumPoints',10);
>> [x,v,s,h,f]=cont(@equilibrium,x0,[],opt);
first point found
tangent vector to first point found
label = H , x = ( 0.005604 -0.001000 0.002702 -0.589286 )
First Lyapunov coefficient = -4.548985e-001

elapsed time  = 0.1 secs
npoints curve = 10
A Hopf point is found for $\nu=-0.58928$ for the state values $x=0.0056037$, $y=-0.001$, $z=0.0027021$. We start a curve of periodic orbits from this Hopf point, using $N=25$ test intervals and $m=4$ collocation points and choosing $\nu$ as the free parameter. The results can be seen in Figure 23.
>> x1=x(1:3,s(2).index);p(6)= x(end,s(2).index);
>> [x0,v0]=init_H_LC(@torBPC,x1,p,[6],0.0001,25,4);
>> opt=contset; opt= contset(opt,'Singularities',1);
>> opt=contset(opt,'MaxNumPoints',50);
>> opt=contset(opt,'Multipliers',1);
>> [x,v,s,h,f]=cont(@limitcycle,x0,v0,opt);
first point found
tangent vector to first point found
Limit point cycle (period = 8.411855e+000, parameter = -5.844928e-001)
Normal form coefficient = 1.788080e-001
Neimark-Sacker (period = 8.861103e+000, parameter = -5.957506e-001)
Normal form coefficient = 2.674115e-003
Period Doubling (period = 9.256846e+000, parameter = -6.146817e-001)
Normal form coefficient = -6.068973e-003

elapsed time  = 28.3 secs
npoints curve = 50
>> plotcycle(x,v,s,[1 2]);
Figure 23: Computed limit cycle curve started from a Hopf bifurcation
\includegraphics[scale=0.7]{ex/Figure23.eps}

The previous computations are done by running the script testtorBPC1, including drawing Figure 23 in which the Hopf point and the three bifurcations of limit cycles are clearly visible. The axis labels were added manually.

In particular, we detect a torus bifurcation point at $\nu=-0.59575$. To recover the torus bifurcation point of (80) we continue the torus bifurcation in two parameters $\nu,\epsilon$. With the aid of the user function $\epsilon$ we will now locate a NS bifurcation of (80). The starting vector x1 is calculated from the NS on this branch using init_NS_NS. Continuation is done using a call to the standard continuer with neimarksacker as curve definition file:

p=[0.5;-0.6;0.6;0.32858;0.93358;-0.9;0.001];
x=[0.00125;-0.001;0.00052502];
[x0,v0]=init_EP_EP(@torBPC,x,p,[6]);
opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxNumPoints',10);
[x,v,s,h,f]=cont(@equilibrium,x0,[],opt);
x1=x(1:3,s(2).index);
p(6)=x(end,s(2).index);
[x0,v0]=init_H_LC(@torBPC,x1,p,[6],0.0001,25,4);
opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'Multipliers',1);
opt=contset(opt,'MaxNumPoints',50);
[x,v,s,h,f]=cont(@limitcycle,x0,v0,opt);
[x1,v1]=init_NS_NS(@torBPC,x,s(3),[6 7],25,4);
opt=contset;
opt=contset(opt,'VarTolerance',1e-4);
opt=contset(opt,'FunTolerance',1e-4);
opt=contset(opt,'Userfunctions',1);
UserInfo.name='epsilon0';
UserInfo.state=1;
UserInfo.label='E0';
opt=contset(opt,'UserfunctionsInfo',UserInfo);
opt=contset(opt,'Backward',1);
opt=contset(opt,'MaxNumPoints',16);
[xns1,vns1,sns1,hns1,fns1]=cont(@neimarksacker,x1,v1,opt);
plotcycle(xns1,vns1,sns1,[1 2]);

This continuation is done by running the script testtorBPC2. The output is the following:

>> testtorBPC2
first point found
tangent vector to first point found
label = H , x = ( 0.005604 -0.001000 0.002702 -0.589286 )
First Lyapunov coefficient = -4.549030e-01

elapsed time  = 0.7 secs
npoints curve = 10
first point found
tangent vector to first point found
Limit point cycle (period = 8.411870e+00, parameter = -5.844928e-01)
Normal form coefficient = 1.788366e-01
Neimark-Sacker (period = 8.861100e+00, parameter = -5.957504e-01)
Normal form coefficient = 2.674115e-03
Period Doubling (period = 9.256846e+00, parameter = -6.146817e-01)
Normal form coefficient = -6.068982e-03

elapsed time  = 11.5 secs
npoints curve = 50
first point found
tangent vector to first point found
label = E0, x = ( 0.046835 0.141698 0.046209 ... 
... 0.141698 0.046209 8.793572 -0.591612 -0.000006 0.822406 )

elapsed time  = 9.1 secs
npoints curve = 16
We note that the point with label E0 is a torus bifurcation point of (80) with $\nu=-0.591612,$ $\epsilon=-0.000006\approx 0,$ and period $0.822406.$ The orbits are shown in Figure 24. The axis labels were added manually.
Figure 24: Computed torus bifurcation curve started from a torus bifurcation of limit cycles
\includegraphics[scale=0.8]{ex/Figure24.eps}

Finally, we continue numerically the NS orbit with two free parameters $\beta,\nu$ and find, interestingly, that the NS orbit shrinks to a single point. The results are plotted using the standard plot function plotcycle where the fourth argument is used to select the coordinates. The plot is shown in Figure 25 (the axis labels were added manually). This computation is executed by running testtorBPC3:

p=[0.5;-0.6;0.6;0.32858;0.93358;-0.9;0.001];
x=[0.00125;-0.001;0.00052502];
[x0,v0]=init_EP_EP(@torBPC,x,p,[6]);
opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'MaxNumPoints',10);
[x,v,s,h,f]=cont(@equilibrium,x0,[],opt);
x1=x(1:3,s(2).index);
p(6)=x(end,s(2).index);
[x0,v0]=init_H_LC(@torBPC,x1,p,[6],0.0001,25,4);
opt=contset;
opt=contset(opt,'Singularities',1);
opt=contset(opt,'Multipliers',1);
opt=contset(opt,'MaxNumPoints',50);
[x,v,s,h,f]=cont(@limitcycle,x0,v0,opt);
[x2,v2]=init_NS_NS(@torBPC,x,s(3),[1 6],50,4);
opt=contset; opt=contset(opt,'VarTolerance',1e-4);
opt=contset(opt,'FunTolerance',1e-4);
opt=contset(opt,'MaxNumPoints',50);
[xns2,vns2,sns2,hns,fns]=cont(@neimarksacker,x2,v2,opt);
plotcycle(xns2,vns2,sns2,[1 2]);
The output is as follows:
>> testtorBPC3
first point found
tangent vector to first point found
label = H , x = ( 0.005604 -0.001000 0.002702 -0.589286 )
First Lyapunov coefficient = -4.549030e-01

elapsed time  = 0.0 secs
npoints curve = 10
first point found
tangent vector to first point found
Limit point cycle (period = 8.411870e+00, parameter = -5.844928e-01)
Normal form coefficient = 1.788366e-01
Neimark-Sacker (period = 8.861100e+00, parameter = -5.957504e-01)
Normal form coefficient = 2.674115e-03
Period Doubling (period = 9.256846e+00, parameter = -6.146817e-01)
Normal form coefficient = -6.068982e-03

elapsed time  = 14.1 secs
npoints curve = 50
first point found
tangent vector to first point found

elapsed time  = 124.5 secs
npoints curve = 50
Here $xns2$ is a $607 \times 50$ matrix. The last four components of each column are, in that order, the period of the orbit, the values of the two free parameters, and the value of the additional variable $\kappa$ that is introduced in (75).

Figure 25: Computed torus bifurcation curve that shrinks to a single point
\includegraphics[scale=0.8]{ex/Figure25.eps}