- Information
- AI Chat
Matcont tutorial
University of Florida
Preview text
Matcont Tutorial: ODE GUI version
Hil Meijer
Exeter, Feb, 2014
‘ ‘If you want to get credit for solving a complicated mathematical problem, you will have to provide a full proof. But if you’re trying to make something as easy as possible, you want to make it foolproof–so simple even a fool could couldn’t screw it up.” 1
1 Introduction
This tutorial tries to explain the basics of how to use the numerical bifurcation package MATCONT by going through an example. We aim at explaining how you can use the software so we assume a basic knowledge of bifurcation theory. As with many things one learns, you will only master it by performing the procedures yourself. So hands on! We try to cover the following procedures in this tutorial
Installation, including steps on mex-files
System definition
Simulation
Continuation of equilibria in one parameter
Continuation of codim 1 bifurcations of equilibria in two parameters
Starting Limit Cycle from a Hopf point
Starting Limit Cycle from a time simulation
Computing homoclinic orbits
Optional Computing PRCs and their derivatives
This is a condensed overview of the capabilities of MatCont. It is not intended to be a course on Dynamical Systems nor on Numerical Bifurcation Analysis. For the latter, we refer to the lab sessions 1-5 and the manual at sourceforge. Good reference material may be found at
staff.science.uu/~kouzn101/NBA/index
. 1 Every one computer is different from the other. So, although these instructions have been checked step by step, do not be surprised if your computations take a different direction. Just try again but now in the other direction.
2 Getting started/Installation
- Please let Matlab be installed, and download matcont (latest version):
sourceforge/projects/matcont/files/matcont/matcont5p3/
. If you end up at the main page, then go to Browse all files|matcont|matcont5p3 and select the zip file. This is the GUI-version, and should not be confused with the version for maps or command line. Unzip the files to the directory where you want it.
1a. Check that you have a file matcont. If not, you have selected the wrong package.
- Then one should type ”mexext” on the matlab command line. For mexw32, mexw64, mex- maci64 go to sourceforge/projects/matcont/files/Auxiliaries/MEXfiles/ and se- lect the appropriate directory. Download all these files and put them into the directory LimitCycle of your matcont directory.
2a. For other platforms you need a compiler: For Linux you have gcc, or you should look at your package manager. 2b. MAC: If at mex -setup you can select gcc without problems, proceed. Otherwise you have to install the XCode package.
- Switch to the matcont directory and type ”matcont”. Some windows should appear.
Figure 1: A typical matcont startup screen
Now several windows will open with a standard system called adapt2. A typical screenshot is shown in Figure 1. One window is titled matcont and has several menu options. For instance, to end your matcont session, choose ‘Select’ in the matcont window and then ‘Exit’. Hereafter we will indicate this with matcont:Select|Exit. Alternatively one clicks on the close button of the matcont window.
- If your model has some strange behavior like ever increasing variables, check whether you have specified the differential equations in the same order as the order of the coordinates.
3 Time Simulations
Click ‘(matcont):Type|Initial point|Point’ to initialize the computation of an orbit starting from a point. In the matcont window the curve type is now P O, every curve type has a similar meaning. Two extra windows ‘Starter’ and ‘Integrator’ appear, the first is to specify initial condi- tions/points, the second to select the numerical integrator and change its settings. In the starter window set A = 3. 25 , a = 100, Bs = 8. 8 , bs = 20, Bf = 44, bf = 100, e0 = 2. 5 , r = 0. 56 , v0 = 6 , theta = 0. 5 , CC = 220, II = 20. We also want to see how the orbits look like and for this we open a window ‘(matcont):Window|Plot|2D-plot’. In ‘Layout|Variables on axes’ select y 0 as abscissa(standard) and y 1 as ordinate, press ‘OK’ and let the axes range from - to .1 for the abscissa and 0 to 3 for the ordinate using ‘(2Dplot:1):Layout|Plotting region’. In ‘Integrator’ change ‘Interval’ to 1 and ‘MaxStepsize’ to 0.
Figure 3: An orbit converging to a stable fixed point for CC = 220, II = 20. Left: (y 0 , y 1 )-phase plane, Right: simulation showing y 0 (t).
After these preparations we are ready to compute some orbits: Select ‘(matcont): Compute|Forward’. Starting from ~y = (y 0... , y 7 ) = (0,... , 0) converges to (. 0014 , 1. 79 , 4. 24 , 4. 24 , 0 , 0 , 0 , 0). You can monitor the orbit and its final values by opening via numeric window. Select ‘(matcont): Window|Numeric’ and recompute the orbit. You can try also other initial points and check that these also converge to this point. So we have a(n apparent) global attractor. To see how y 0 evolves over time, select ‘(2Dplot:1): Layout|Variables on axes’ and select time and y 0. Re- draw the curve via ‘(2Dplot:1): Plot|Redraw curve’ and adjust the plot-ranges to obtain Figure 3(right). Now change to II = 135 and again compute an orbit via ‘(matcont):Compute|Forward’. The result is quite different. You will have to adjust the axes to see the full range of y 0. You can use the pull-down menu as before, but the zoom-option (looking glass) in the menu bar is also useful, or manually from the command line with ylim([0 0]). The result should like Figure 4(right). We see some periodic behavior that has a frequency of 3Hz and looks like spike-wave discharge patterns seen in EEG.
3 Continuation of Equilibria in one parameter
We are now going to study the effect on the stability of the equilibrium when varying either II. Se- lect the last point via ‘(matcont):Initial point’ and then choose the last point on the previous curve.
Figure 4: An orbit converging to a stable fixed point for CC = 220, II = 135. Left: (y 0 , y 1 )-phase plane, Right: simulation showing y 0 (t).
Change the curve type to EP EP by selecting ‘(matcont):Type|Initial point|Equilibrium’. Ac- tivate the parameter II by clicking on the button next to it. In the 2Dplot:1 window change the abscissa to the parameter II via ‘(2Dplot:1):Layout|Variables on axes’. Also change the plotting region for II to -250 to 300 and for y 0 to 0 to 0. To speed up the continuation we will set ‘Continuer):MaxStepsize=5’ ‘Continuer):MaxNumPoints=600’. The stability changes upon the variation of II thus we would also like to see how the eigenvalues of the Jacobian at the equilibrium develop. For this we select ‘(Numeric):Window|Layout’ and highlight the option ‘eigenvalues’( it changes to capitals). Since we have (too) parameters you may want to unselect the parameters (changes it to small letters). Press ‘OK’. Select ‘(matcont): Compute|Forward’. NOTE: In the lower left corner a small window appears with buttons ‘Pause’,‘Resume’ and ‘Stop’. Whenever a bifurcation is detected, MATCONT pauses and shows some relevant informa- tion. In the sequel a few such points will be detected, but we will not always indicate to resume the computation. At II = 109 the message ‘Limit point’ appears in the status field of the matcont window. Here one eigenvalue has zero real part: Check that in the numeric window!. The curve(branch) of equilibria has a turning point, on one side the equilibria are stable on the other unstable. Press the ’Resume’ button in the lower left of your screen or press space to continue the computation of the now unstable branch. In the meantime MatCont find a special point, i. a test function is zero, but that is due to two real eigenvalues, one positive and one negative, summing up to zero, not a complex pair. This is mentioned in the main MatCont window as a neutral saddle. At II = − 211 .513481 another Limit point is found. Here the branch does not become stable, just one more eigenvalue has now positive real part. When II = − 188 .997488 is reached two complex conjugate eigenvalues have zero real part and a Hopf bifurcation occurs. The 2Dplot should now look like Figure 5. If the computation would be fast to follow, you can select to pause the computation of the EP EP curve after each computed point. For this choose ‘(matcont):Options|Pause|At each point’. (Do not forget to change it to ‘At special points’ after this computation. Now we turn to the command line, i. the Matlab window, where you most probably now have the following information:
first point found tangent vector to first point found label = LP, x = ( 0 0 6 0 5 0 5 0 109.
Select one of the Limit points found in the equilibrium continuation as initial point ‘(matcont):Select|Initial point’. Activate II and CC as active parameters. Close the 3D-plot and re-open the 2D-plot. Change the ordinate of the 2D-plot to the parameter CC and its range from 0 to 250. Now select ‘(matcont):Compute|Forward’ and when finished also ‘Backward’. At (CC, II) ≈ (59. 11 , 168 .7) a Cusp bifurcation (CP) and at (II, CC) = (1. 6840517 , 1 .6829686) a spurious Zero-Hopf bifur- cation is detected. You have now computed the wedge, see Figure 6 within which you have bistability. This is not the full story as we also found already Hopf bifurcations, implying that it is more co-existence and not bistability of equilibria.
Figure 6: Limit Point bifurcation curves in the (II, CC)-plane with codimension 2 points with a wedge of co-existence of equilibria.
To compute a curve of Hopf bifurcations: select a Hopf point analogously as for the Limit point as initial point. Change the curve type to H H via ‘(matcont):Type|Curve|Hopf’. Take here II and CC as active parameters. Press ‘Compute|Forward’ and ‘Compute|Backward’ to obtain the diagram as in Figure 7.
4 Continuation of Limit Cycles
4 Starting from a Hopf point
We are now going to compute a branch of limit cycles starting from the Hopf point. Select the second Hopf point of the curve ‘equi’ using ‘(matcont):Select|Initial point’ as new initial point. Activate the parameter II and the period by clicking on the buttons in the Starter window. Next we select ”yes” for Fold under the heading monitor singularities in the starter window. By default, and for good reasons (many spurious warnings otherwise), detection of bifurcations is switched off. To speed up computations we will set MaxStepsize=200. This is large, but for this tutorial this is acceptable. Close the 2DPlot (for speed reasons: plotting the cycles takes too much time now) and ‘(matcont):Compute|Forward’. At II = 60 we see a Fold bifurcation and the periodic orbit becomes unstable. Resume the computation and stop when you are near II = −188. Here the branch terminates in the subcritical Hopf bifurcation. Open a 3D plot with ‘(matcont):Window|Plot|3D-Plot’. Choose
Figure 7: Limit Point and Hopf bifurcation curves in the (II, CC)-plane with codimension 2 points. The ZH is spurious, the GH indicates where the Hopf bifurcation changes from sub- to supercritical. The Hopf curve emerges from a Bogdanov-Takens point. To detect the lower BT-point, you should set the MaxStepsize to 0 in the Continuer window.
OX=z, OY=F , OZ=y in the draw3 attributes window. Change the plotting the region to II ∈ [− 200 , 100], y 0 ∈ [0, .15], y 1 ∈ [0, 30]. Your 3D-plot will resemble the one in Figure 8.
Figure 8: Equilibria and Limit cycles in (II, y 0 , y 1 )-space. The points and limit cycles which undergo a bifurcation are indicated with red dots.
s A structure that specifies special points, i. detected bifurcations and the first and last point.
5 Plotting the profile
To plot the profile, i. the evolution of a variable during the whole period, is nontrivial as it is stored in a condensed way. Suppose we know how many mesh-points were used. This can be determined from the f variable, simply look up the first entry equal to 1. Now we need
the index of the periodic orbit we want to plot: say ii.
the profile of one variable, say the jj-th: x(jj:ndim:end-2,ii)
the time-mesh tau=f(1:ntst+1,ii), with 0 = τ 0 < τ 1 < ... < τntst = 1.
the period T= x(end-1,ii) to scale the mesh back from [0,1] to [0,T]
So we could plot:
ii=50;jj=1;ntst=40;ndim=8; plot(f(1:ntst+1,ii)*x(end-1,ii),x(jj:ndim:end-2,ii))
You will get an error message because the computed curve x contains also data about the orbit between mesh-points. While you are trying to plot the profile only coarsely at the mesh-points. A first correction would be
ii=50;jj=1;ndim=8;ntst=40;ncol=4; plot(f(1:ntst+1,ii)x(end-1,ii),x(jj:ndimncol:end-2,ii))
This is usually not good enough, so we actually want the detailed grid. For this we need to know that between to mesh-points the orbit is stored at equi-distant time-points ti,j = ti + mj (τi+1 − τ ), j = 0, 1 , ..., m. So the correct commands are:
ii=10;jj=1;ndim=8;ntst=40;ncol=4; mesh=sort([reshape(repmat(f(1:ntst,ii),1,ncol)+... repmat((0:ncol-1)/ncol,ntst,1).repmat(diff(f(1:ntst+1,ii)),1,ncol),ntstncol,1); 1]); plot(mesh*x(end-1,ii),x(jj:ndim:end-2,ii))
5 Identifying Spike Wave Discharges
One particular feature of this model as that it has Spike-Wave Dynamics. That is, on top of the 3Hz wave, there is an additional spike that can be detected as a local extremum. This is also called a ‘false’-bifurcation. Here we will compute these local extrema.
load Systems\JR_slow\diagram\lc1 out=x(3:8:end-2,:)-(x(5:8:end-2,:)+x(7:8:end-2,:))/2; extrema=nan(6,size(out,2)); for i=1:size(out,2) diff_tr=diff(out(1:end-1,i)); diff_tr_shift=diff(out(2:end,i)); temp=[out(find([ (diff_tr>0).(diff_tr_shift<0) ]),i);... out(find([ (diff_tr<0).(diff_tr_shift>0) ]),i)]; extrema(1:size(temp,1),i)=temp; end figure(1);hold on;plot(x(end,:),sort(extrema));
Do this again for the part lc2 and your plot will look like Figure 9.
Figure 9: Extrema in the output signal. The appearance of local extrema indicate spikes in the waveform. Near the fold bifurcations the figure suffers from the large stepsizes.
6 Connecting Orbits
In the following instructions we will see how to start the computation of homoclinic orbits. The main difficulty is to get good starting data for the continuation. We will show two methods that may be useful to accomplish this.
6 Homoclinic orbits starting from a limit cycle
Sometimes when following a limit cycle you will notice that the parameters do not change anymore. Meanwhile the period goes to infinity. Typically, in a graphic window plotting the continuation parameter versus the period you see something as in Figure 10right. In such cases you may suspect that you are close to a connecting orbit. You will have to look at the orbit to be sure. But suppose you notice that the orbit is close to a steady state at some point, wouldn’t it be nice to take this as starting data for the computation of connecting orbit? This method is implemented in Matcont for homoclinic orbits only. Along a limit cycle we can make a guess for an equilibrium by checking at which point of the cycle the vectorfield has smallest norm. This goes well for a single equilibrium, but not for two. We will use the last point on the curve lc1, which you already computed. Take the last point on that curve via ‘matcont:Select|Initial Point’ and browse in the menu. Change the curve type ‘matcont:Type|Curve|’ and select “homoclinic to saddle”. The Starter window should show 108. 4 < II < 108 .6 and the period is at least 0. If not, select Compute|Forward or backward to make sure you are closer to the homoclinic bifurcation. Select the last point of the LC LC curve just computed as initial point. Next select Curve|Homoclinic to saddle and select II, CC and eps 0 , eps1 as free parameters. It turns out that eps0 and eps1 are the correct auxilary parameter to free, but in general this needs experience or trial-and-error, by trying and singlet or doublet of T,eps0,eps1 sequentially. This is the hardest part for the continu- ation. Set Maxstepsize to 200 again in the continuer window and Compute|Forward. Compute a few points and stop the continuation when you want to see results. For this open a 2Dplot and choose Layout|Variables on axes to II and CC for the abscissa and ordinate, respectively, and change Layout|Plotting region the range to [0, 200] × [0, 250] and zoom in as appropriate
If you are interested in just one PRC, just set MaxNumPoints in the Continuer window to
The derivative is also computed as it is related to phase-locking and synchronization, see [3].
References
[1] New features of the software MatCont for bifurcation analysis of dynamical systems. A. Dhooge, and W. Govaerts, Yu. Kuznetsov, H.G. Meijer and B. Sautois, Mathematical and Computer Modelling of Dynamical Systems, Vol. 14, No. 2, pp 147-175, 2008.
[2] Goodfellow M., Schindler K. and Baier G., Intermittent spikewave dynamics in a heteroge- neous, spatially extended neural mass model, NeuroImage 2011, Vol 55(3), pp. 920-932.
[3] W. Govaerts and B. Sautois, Phase Response Curves, Delays and Synchronization in Matlab, Lecture Notes in Computer Science 3992 (2006) 391-398.
[4] AJ Homburg and B Sandstede: Homoclinic and heteroclinic bifurca- tions in vector fields. Handbook of Dynamical Systems III (Edited by H Broer, F Takens and B Hasselblatt). Elsevier (2010) 379-524. available at dam.brown/people/sandsted/publications/survey-homoclinic- bifurcations
Matcont tutorial
Course: Nonlinear Control (EML 6350)
University: University of Florida
This is a preview
Access to all documents
Get Unlimited Downloads
Improve your grades
This is a preview
Access to all documents
Get Unlimited Downloads
Improve your grades
Why is this page out of focus?
This is a preview
Access to all documents
Get Unlimited Downloads
Improve your grades
Why is this page out of focus?
This is a preview
Access to all documents
Get Unlimited Downloads
Improve your grades
Why is this page out of focus?
This is a preview
Access to all documents
Get Unlimited Downloads
Improve your grades