MATH3888 M-Lab 5:
One-parameter continuation with MatCont
The purpose of this M-Lab is to familiarise ourselves with the numerical bifurcation package
MatCont.1 By the end of the M-lab you should know how to:
(i) Define a system of ODEs in MatCont;
(ii) Numerically ‘integrate’ solutions using standard numerical ODE solvers;
(iii) Visualise solutions graphically, as well as numerically;
(iv) Identify equilibria numerically;
(v) Numerically ‘continue’ equilibria and limit cycles in 1 parameter.
We assume that MATLAB and MatCont are already installed on your personal computer.
If they are not, do this now by following the installation steps outlined in the “Installing
MatCont” document available on the course webpage (Ed Resources).
To open MatCont, first open MATLAB, make sure that the current working folder is your
MatCont folder, and type ‘matcont’ into the MATLAB command line. A window titled
‘MatCont GUI’ should appear.
1 Example: the FitzHugh-Nagumo model
We are now going to define a system of ODEs in MatCont. We will use the FHN model
dv
dt
= (f(v)− w + Iapp)/ε
dw
dt
= v − γw − v0,
(1)
with
f(v) = v(1− v)(v − 0.2)
and the parameter set
ε = 0.2, Iapp = 0, γ = 1, v0 = 0 ,
which you all studied with BDToolbox last week.
1These notes are primarily a compilation of notes by Yu. A. Kuznetsov (Utrecht University) and
N. Neirynck (Ghent University), which can be found at https://webspace.science.uu.nl/~kouzn101/NBA/
1
1.1 Defining a system of ode’s
In the MatCont GUI, click on ‘Select’, then ‘System’, then ‘new’. Hereafter we will write a
sequence of commands like this as Select|System|New. This opens the System window,
which contains several fields and buttons. To identify the system, type for example
FHN
in the Name field (it must be one word).
Input names of the Coordinates: w,v, the Parameters: Iapp,epsilon,gamma,v0 and
the name for the Time: t (default).
If shown, select symbolic generation of the 1st, 2nd and 3rd order derivatives2; see Figure
1.
Figure 1: Specifying a new model.
Finally, in the large input field, enter the FitzHugh-Nagumo system as
w’=v-gamma*w-v0
v’=(v*(1-v)*(v-0.2)-w+Iapp)/epsilon
Avoid typical mistakes:
• Make sure to write multiplication with a *.
• Specify differential equations in the same order that you specify coordinates.
Your system should look as the one in Figure 1. Once it does, press OK. Assuming
you made no typing mistakes, the System window disappears, and you will see in the main
2Symbolic computation of derivatives is usually preferred for ‘simple enough’ systems. For the systems we
consider, choosing symbolic computation of 1st, 2nd and 3rd (but not 4th and 5th) is a pretty good rule of
thumb for good numerical performance.
2
MatCont window that FHN becomes the current System of MatCont. If selected, the
information field Derivatives shows the string SSSNN meaning that the symbolic 1st, 2nd
and 3rd order derivatives of the RHS will be used.
If you want to edit an existing system, click Select|Systems|Edit/Load/Delete, select
this system in the appearing Systems window, and press Edit button. The inputed system
can be completely erased by selecting Action|Delete there.
1.2 Selection of solution type
To tell MatCont that we want to integrate the system, i.e., to numerically solve the ini-
tial value problem for it, we have to specify the type of the initial point and the type of
the curve to compute. To select the initial point type, input Type|Initial Point|Point,
which means that the initial point has no special properties. To select the curve type, click
Type|Curve|Orbit. The information fields in the MatCont main window will reflect the
selections, and two more windows appear. These are the corresponding Starter and Inte-
grator windows. Move them, if necessary, to make both visible.
1.3 Setting initial data for integration
The Starter window is curve-dependent and is used here to specify initial data for the
integration. Let us input the following initial coordinate and parameter values (Figure 2):
w 0.1
v 0.4
Iapp 0.0
epsilon 0.2
gamma 1.0
v0 0.0
Figure 2: Starter window for integration.
We will use the default numerical parameters of the Integrator, except of the Interval.
Input in the Integrator window:
Interval 100
3
Figure 3: Integrator window.
to get Figure 3. Check that the default method of integration is ode45. You can inspect
available methods by clicking the Method menu button.
1.4 2d visualisation
To visualise the orbits, we have to open at least one more window and setup plotting at-
tributes.
1.4.1 2D Graphic Window
Select Window/Output|Graphic|2D plot in the main MatCont window. A new Plot2D
window appears.
Using MatCont|Layout, set w as the Abscissa and coordinate v as the Ordinate and
input new visibility limits along the axes:
w -0.1 1.0
v -1.0 1.0
The layout window is presented at the right of Figure 4.
Figure 4: Coordinates and range limits in a 2D plot.
4
1.5 Integrating orbits
Now we are ready to start numerical integration. Input the Compute|Forward menu com-
mand in the main window to start integration. A Control window is opened automatically.
The computation can be paused or stopped by clicking the corresponding button in that
window. Also, you can interrupt the computation at any time by pressing the Esc key on
the keyboard.
The integration will start and shortly after produces Figure 5. The orbit tends to a stable
equilibrium.
Figure 5: An orbit in the Plot2D window converging to a stable equilibrium.
Figure 6: The Control window when the computation is finished.
When the computation is finished, the Control window looks like Figure 6. You can
now either close the Control window or press the View Result button to inspect the output
in a Data Browser window, see Figure 7. By pressing the View CurveData button in
the Data Browser window one obtains a spreadsheet view of all points computed during
the time integration.
1.6 Continuation of Equilibria
We can continue the equilibrium found by integration with respect to the parameter Iapp.
Open the Data Browser window by pressing Select|Initial Point in the main MatCont
5
Figure 7: Inspect the output of the computation in a Data Browser window.
window. Then click on Last Point and Select Point, as shown in Figure 7.
To tell MatCont to continue an equilibrium curve, we have to specify the initial point
type and the curve type. Selecting Type|Initial Point|Equilibrium we set the point type
EP and the same (default) curve type, so that we prepare to compute the EP EP equilibrium
curve, as indicated in the main MatCont window. The Starter window is now modified and
the Integrator window is replaced by a Continuer window.
The initial values of the parameters and the equilibrium coordinates are visible in the
Starter window. Press the button next to parameter Iapp, indicating that it will be a
bifurcation parameter (i.e. activate parameter Iapp). The resulting Starter window is
shown in Figure 8. There you can also see which singularities will be monitored (all), and
whether the equilibrium eigenvalues will be computed (yes).
Figure 8: Starter window for the equilibrium continuation: Parameter Iapp is activated.
In the Continuer window, several default numerical parameters related to the continu-
ation are listed. Change only (see Figure 9)
MaxStepSize 0.02
6
Figure 9: Continuation parameters in the Continuer window.
Use MatCont|Layout menu in the Plot2D window to select the parameter Iapp as
abscissa and the coordinate v as ordinate with the visibility limits:
Abscissa -0.0 1.0
Ordinate -0.5 1.5
and clear the window by pressing MatCont|Clear in the Plot2D window. Now click
Compute|Forward to get part of a curve as in Figure 10.
Figure 10: Equilibrium manifold in the (Iapp, v)-plane: H denotes an Andronov-Hopf point.
During the computations a Control window is opened and the continuation is stopped
each time when a singularity is found. Then press Resume in the Control window to
advance the continuation. Along the forward branch, MatCont stopped at two Andronov-
Hopf points labelled H.
7
After finishing this continuation press the View Result button in the Control window.
This opens a Data Browser window. Press the View CurveData button at the bottom
of the Data Browser window. This opens a spreadsheet view of all computational values
that will be stored as data of the curve.
To determine stability of the equilibria and read the bifurcation parameter values, open
Window/Output|Numeric and select Layout in the Numeric window.
The Layout menu offers the choice to visualise during the continuation the coordinates
(i.e. state variables), active parameters, testfunctions, eigenvalues, current stepsize, user-
functions (when userfunctions are present), and current number of computed points. Check
only the Coordinates, Parameters and Eigenvalues.
This makes the state variables, active parameters and eigenvalues of the equilibrium
visible in the Numeric window (upon resizing it).
Clear the Plot2D window with MatCont|Clear and repeat Compute|Forward. In
the Numeric window you can read the H parameters:
Iapp = 0.23006 . . . , Iapp = 0.47393 . . .
(see Figure 11; your actual numerical values might slightly differ). Note, the real part of the
complex conjugate eigenvalue is (close to) zero in each case.
Figure 11: Numeric window at the first and second Andronov-Hopf point.
In the MATLAB Command Window, the value of the first Lyapunov coefficient is shown
at each Andronov-Hopf point:
label = H, x = ( 0.236700 0.236700 0.230069 )
First Lyapunov coefficient = -4.999990e+00
label = H, x = ( 0.563299 0.563299 0.473930 )
First Lyapunov coefficient = -5.000002e+00
Since each Lyapunov coefficient is negative, each Andronov-Hopf bifurcation is supercrit-
ical. Finally, clear the Plot2D window once more and recompute the equilibrium curve for-
ward with Options|Suspend Computation mode At Each Point. Resume computations
after each computed point and monitor the eigenvalues in the Numeric window. Restore
the original pause option (At Special Points) via Options|Suspend Computation and
close the Numeric window.
8
1.6.1 Andronov-Hopf bubble
Next, we select one of the Andronov-Hopf bifurcations via Select|Initial Point, by load-
ing a point labeled H within a curve labeled EP EP. Notice the Initializer changes auto-
matically to init H LC mode, indicating that MatCont will continue limit cycles. Use
Compute|Forward to fill out the family of limit cycles emerging at one AH bifurcation
and terminating at the other.
Figure 12: Andronov-Hopf bubble in the (Iapp, v)-plane: H denotes an Andronov-Hopf point.
2 Example: cusp normal form
Your aim is to analyse the following system with Matcont:
dx
dt
= α+ βx− x3
dy
dt
= −y − x2,
(2)
with the starting parameter set α = −0.4, β = 0.5.
(a) Implement system (2) in MatCont.
(b) Use initial conditions (x0, y0) = (0.7,−0.3) and integrate the system forward in time
to detect a stable equilibrium. Make sure that you have converged sufficiently close.
(c) Continue this stable equilibrium with respect to the parameter α ∈ (−0.5, 0.5). What
bifurcation(s) do you detect?
(d) Save the bifurcation diagram as a png-file.
Remark: Please have a look at the online material at https://webspace.science.uu.nl/
~kouzn101/NBA/.
9