MATH 437/537 – XPP/AUTO Tutorial
by Sam Jamaleddine
Winter 2023
Department of Physiology, McGill University, Montreal, Canada
Contents
1 Introduction 2
2 Creating an ODE file 2
3 Opening XPP/AUTO 4
3.1 Running XPP/AUTO on Windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 Running XPP/AUTO on Mac or Linux . . . . . . . . . . . . . . . . . . . . . . . . . . 5
4 Simulations and analysis in XPP 6
4.1 Time series solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
4.2 Phase portrait . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.3 Determining equilibria and plotting manifolds . . . . . . . . . . . . . . . . . . . . . 12
5 Parameter bifurcations in AUTO 14
5.1 Single-parameter steady-state bifurcations . . . . . . . . . . . . . . . . . . . . . . . . 14
5.2 Continuation of periodic solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
5.3 Two-parameter bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
XPP/AUTO Tutorial
1 Introduction
To the experienced user, XPP/AUTO can be a quick and efficient way to analyze the basics of
a nonlinear system of ordinary differential equations (ODEs). However, the learning curve can
be quite steep at first and the user interface may sometimes crash and/or behave unpredictably,
all of which can easily frustrate a new user. The goal of this document and tutorial is to walk
through some of the essential functions that can be carried out using XPP/AUTO, and should
hopefully offer some clarity into its usage.
I will outline some of the key points that you need to get you started with XPP and to use it
for the course, although there are many more things that you can do with it that I won’t cover.
Don’t worry too much if it’s understandably overwhelming at first; with enough practice, the
keyboard shortcuts, optimal parameters, etc. become engraved in your brain and eventually
they become automatic without having to even think about them anymore.
Make sure to first download XPP/AUTO at http://www.math.pitt.edu/~bard/xpp/xpp.html
(follow the download instructions specific to your operating system).
If you have any further questions about the program that are not covered here, or if you
would like more in-depth explanations, or if you want to see more examples, please refer
to the online documentation at http://www.math.pitt.edu/∼bard/xpp/help/xpphelp.html
2 Creating an ODE file
In order to analyze the system of ODEs in question, you need to write a special file called an
ODE file (file extension: .ode). You can write one using a standard text file editor, making sure
to save it as using a .ode file extension (I recommend using Sublime Text but use whichever text
editor you like).
Say we want to analyze the following system of ODEs describing a simple SIR model (standard
epidemiological model of susceptible, infected, and recovered individuals), with the R variable
2
XPP/AUTO Tutorial
set to quasi-steady state:
dS
dt
= pi − aSI − µS (1)
dI
dt
= aSI − bI − µI (2)
dR
dt
= bI − µR = 0 ⇒ R(t) = b
µ
I(t) (3)
The system described by Eqs. (1)-(3) contains two differential equations (S and I with respect to
time), one so-called “auxilliary” variable (R with respect to I), and 4 constant parameters (pi, a,
b, and µ). A typical ODE file written to analyze this system in XPP/AUTO could look like this
(see SIR.ode):
1 # Equations
2 S' = p - a*S*I - mu*S
3 I' = a*S*I - b*I - mu*I
4
5 # Functions
6 aux R(I)=b*I/mu
7
8 # Parameters
9 par p=3
10 p a=1, mu=1
11 p b=1
12
13 # Numerics
14 @ TOTAL=100,DT=.01,xlo=0,xhi=100,ylo=-0.5,yhi=5
15 @ NPLOT=1,XP1=t,YP1=S
16 @ MAXSTOR=10000000
17 @ BOUNDS=100000
18 @ dsmin=1e-5,dsmax=.1,parmin=-.5,parmax=.5,autoxmin=-.5,autoxmax=.5
19 @ autoymax=.4,autoymin=-.5
20
21 # Initial conditions
22 S(0)=0.5
23 I(0)=0.5
24
25 done
3
XPP/AUTO Tutorial
• Equations: this is where you put your differential equations. You can write these in dif-
ferent ways, either as S’ = ... (as shown above), or more explicitly as dS/dt = ...; both
work exactly the same.
• Functions: these are the auxilliary functions, i.e. the quantities which are functions of one
or more of your time-dependent variables (in this case R is a function of I(t) as in Eq. (3)).
Including the aux prefix allows you to still plot R in XPP (but not in AUTO). WARNING:
Do not forget to write R(I) (instead of just R); XPP will still compile your file, and it will
seem like it works, but sometimes it won’t! (learned the hard way)
• Parameters: these are the constant parameters of your model (4 of them in this case). You
can start the line with either p or par or even parameter, and can put multiple parameters
on one line separated by commas as I did on line 10.
• Numerics: These let you modify the internal settings within XPP/AUTO that have to do
with various things such as how you want XPP to display time series/bifurcation plots
(which variables to plot, the range, etc.), and also numerical parameters like min/max
step size (dsmin, dsmax), upper bound, etc. We’ll go over these in more detail later. You
technically don’t need to set any of them beforehand since they can all be changed within
the program, but often it is helpful to do so since you save yourself the trouble of resetting
them manually each time you open a file. IMPORTANT: exceptionally in the numerics
settings, you cannot add spaces before or after the commas and equal symbols.
• Initial Conditions: Set the value of your time-dependent variables at the first time point,
i.e. at t = 0. If you leave these out, they default to 0. You could also write init S=0.5, e.g..
Note that variable and parameter names are NOT case sensitive, so for example you cannot
name one variable as R and another as r because XPP will register these as the same variable/-
parameter. Consequently, you should not name a variable T, because T (i.e., t) is reserved for
the independent variable (time). Other keywords (like ’pi’) are also reserved for other purposes.
There are a lot more things you can write into .ode files that XPP can handle (e.g. step functions,
impulse functions, delay differential equations, noise, and much more). For full documentation,
see http://www.math.pitt.edu/~bard/xpp/help/xppodes.html.
3 Opening XPP/AUTO
How you go about opening a .ode file may vary depending on your operating system, i.e. if
you are using a Windows, Mac, or Linux machine. First, make sure you’ve followed specific
4
XPP/AUTO Tutorial
instructions for your OS outlined at http://www.math.pitt.edu/~bard/xpp/xpp.html.
Note that sometimes, the file will fail to open because of syntax or other errors in your .ode file.
If this is the case, make sure you read the output text on the XPP Shortcut/Terminal window,
there should be an explanation of the error at the end of this text.
3.1 Running XPP/AUTO on Windows
If you are running a Windows machine, then you first need to run Xming before opening your
.ode file (attempting to run the ODE file before doing so will result in an error). Then, in order
to start XPP/AUTO, you can click and drag the .ode file, e.g. from your File Explorer, and
drop it into the XPP Desktop shortcut (instructions for setting up the Desktop shortcut found at
http://www.math.pitt.edu/~bard/xpp/installonwindows.html).
So, for example, to open the ODE file called SIR.ode, first ensure that Xming is running, then
drag the file SIR.ode and drop it into the XPP Shortcut on the Desktop.
3.2 Running XPP/AUTO on Mac or Linux
To run the program if using a Linux machine or a more recent version of Mac OS, you should
go through the system terminal. Notably, for Mac users, the ability to drag and drop ODE
files into the Desktop shortcut is no longer available (because Mac stopped supporting 32-bit
programs for some reason). If you can get it to work, great, but if not, follow the instructions
below.
First, open a terminal window (for Mac users that may not be familiar with terminal, do a Spot-
light search and type in ’Terminal’ and it will come up). Then, you’ll want to navigate to the
directory containing the .ode file, using the cd command (‘change directory’). For example, on
my Mac I have the two ODE files, SIR.ode and dengue.ode in a subfolder called ’XPPtutorial’
within my ’Documents’ folder. To open the file named SIR.ode using XPP/AUTO, I would then
type the following commands
1 cd Documents/XPPtutorial
2 /Applications/xppaut SIR.ode
5
XPP/AUTO Tutorial
The first line of code brings me into the directory containing the file of interest; this line of code
will vary from person to person depending where you store the file called SIR.ode. The sec-
ond line varies for Mac and Linux users; for Linux omit /Applications/ and just write xppaut
SIR.ode.
Of course, later on when you want to open different .ode files, make sure you change the name
of the file from SIR.ode to whatever you named the file in question.
A few notes for those of you who might be completely unfamiliar with Shell commands:
• You can always use the ls command in the terminal window; typing “ls” will display all
files and folders within the current directory, which can be useful if for example you forget
the file or folder names. You can use the Tab key to auto-complete in some cases.
• If you want to list only files of a given type (e.g. only the .ode files in the directory) you
can type ls *.ode.
• If you mistakenly went into the wrong directory and want to go back to the previous one,
type “cd ..” (two dots). Alternatively, “cd ~” (tilde) brings you back to your home direc-
tory.
That’s probably as much Shell as you’ll need to know to run XPP/AUTO, but feel free to look
up basic Shell commands for any additional information.
4 Simulations and analysis in XPP
Now that you’ve hopefully got the XPP window open, it should look something like in Fig. 1.
Take a look at 6 buttons in the top row:
• ICs: View and set the initial conditions, restore initial conditions to default values (from
the .ode file)
• BCs: View and modify boundary conditions (won’t be needed in this course)
• Delay: If using delay differential equations, set the initial conditions of delayed variables
(i.e. prior to t=0)
6
XPP/AUTO Tutorial
Figure 1: The XPP window that should appear when the program successfully reads SIR.ode
• Param: View and modify, or reset, model parameters. Tip: if the parameter list is long and
cut off, press the full screen button to show all parameters at the same time.
• Eqns: View model equations (ODEs). Do not attempt to modify the model equations
within XPP, this has never worked for me. If you need to change it, close the program,
modify the source file (.ode), save, and recompile in XPP.
• Data: View the data of each variable at each time point of the integration (note: this will be
blank if you haven’t run the simulation yet).
On the left hand side, you have all the options for running simulations, plotting nullclines, mod-
ifying numerics (e.g. step size, integration method, etc.), modifying view options (e.g. which
variables to plot, time series vs. phase portrait), adding curves, solving for fixed points, plotting
manifolds, and much more. KEYBOARD SHORTCUTS: Note that each button in this column
has one uppercase letter, e.g. (I)nitialconds or ph(A)sespace. You can navigate these different
options and tools using your cursor, but for the most part and as you get used to the program,
you learn to use these uppercase letters as keyboard shortcuts to perform tasks very quickly and
efficiently. We’ll now go through some of the most important of these functions.
4.1 Time series solutions
Notice at the top of the program above the empty plot, the title reads ”S vs T” meaning variable
S with respect to time. This is because under the ’Numerics’ section in the ODE file, on line 15,
7
XPP/AUTO Tutorial
we set XP1=t (i.e. plot time on the x-axis) and YP1=S (plot S on the y-axis). If these are left out,
then XPP will simply plot the first variable in the equations list w.r.t. time.
1. To run the simulation/integrate the solution using the default initial conditions (ICs) and
parameters set in the ODE file, press (I)nitialconds → (G)o, or simply type I then G on
your keyboard. You should get the following trace to appear (Fig. 2).
Figure 2: By default, XPP is plotting YP1 vs. XP1, which we set to S and t, respectively, in the ODE file.
2. To view the trace for I(t) on the same graph, press (G)raphic stuff → (A)dd curve, then
set the Y-axis variable to I, and change the colour by typing in an integer from 1 to 10
(0 is black, you can play around with the colours to see what you get for each number).
Choosing 1 yields red as in Fig. 3.
Figure 3: Notice the new red trace representing I(t)
8
XPP/AUTO Tutorial
3. As you see in Fig. 3, S(t) (black) and I(t) (red) converge to a stable steady-state. A very
useful feature, and one way to see the numerical values of S and I at steady state, is to press
(I)nitialconds→ (L)ast (do this a few times). This makes the ICs of the new simulation the
last data point of the previous simulation, in this case near the steady state, so each time
you perform (I)→ (L) you get closer and closer to the steady state (Fig.4).
Figure 4: Performing the (I) → (L) operation 3 times, solutions converge closer and closer to the true
steady state value at S = 2 and I = 0.5 for the default parameter values. Window at top-left is the ’ICs’
window.
4. The sliders that you see at the bottom of the XPP window allow you to vary up to 3 pa-
rameters at a time on a sliding scale, between a lower and upper range that you set. For
example, try setting a slider for the parameter a between 0 and 1 (Fig. 5), and dragging the
small black square across that range. You should see S(t) and I(t) vary as you vary a.
Figure 5: Slider window to quickly scan over parameter value.
9
XPP/AUTO Tutorial
A few tips:
• To delete any/all the extra curves and prevent them from being plotted in the future (e.g.
the red trace in previous Figs), press/type (G)raphic stuff → (D)elete last, or (G)raphic
stuff→ (R)emove all (to clear everything but the original trace).
• Running new simulations simply plots new traces on top of the old ones. To tidy up,
simply press/type (E)rase and they should be cleared off the plot.
4.2 Phase portrait
Up until now, we’ve only been plotting solutions with respect to time. Often, however (and
especially with 2D systems like the one we effectively have), it is useful to plot the variables with
respect to each other. Before proceeding, make sure you set parameters back to their defaults
(click ’Param’ at the top, and ’Default’ in the pop-up window).
1. To change the axes, press (V)iewaxes→ (2)D. A window will appear allowing you to set
the axes, ranges, and labels; set them as shown below in Fig. 6.
Figure 6: Changing what’s being plotted on the x- and y-axes.
2. In this view, we can do a few cool things. For example, we can get a good idea of the vector
field by pressing (I)nitialconds→ m(I)ce. This allows you to click anywhere on the graph
to plot many traces from different initial conditions, and eventually a vector field reveals
itself as in Fig. 7 and we can see all solutions converging to a single point in the graph.
When you’re satisfied, press the ‘Escape’ key on your keyboard to exit.
10
XPP/AUTO Tutorial
Figure 7: Solution flows in phase space reveals behaviour of system (evolution to stable focus).
3. Let’s plot nullclines! First, we need to tweak some numerics to get smoother, more accurate
estimates from XPP as to where the nullclines actually are. Go into the n(U)merics menu,
and you will see a bunch of options on the left. In this menu, you can change things like
integration method (M), step size ∆t in the integration (D), total integration time (T), upper
allowable bound for variables to assume (B), and more.
4. Right now, we are interested in nullcline numerics; press (N)cline ctrl and increase “ncline
mesh’ at the top from 40 to at least 400 (the larger this value, the more mesh points used
in computing the nullclines and the smoother and more accurate they will be). Then press
’Escape’ to go back to the main menu
5. To plot the nullclines, press (N)ullclines → (N)ew. They should appear as in Fig. 8, with
x-axis nullclines shown in orange and y-axis nullclines in green.
Figure 8: One S (horizontal)-nullcline plotted in orange, two I (vertical)-nullclines in green.
11
XPP/AUTO Tutorial
4.3 Determining equilibria and plotting manifolds
There are a few ways to find steady states. In anticipation of performing continuation analysis
with AUTO (i.e. plotting bifurcation diagrams), it is necessary that the data stored in the ICs
column contains the values of the variables at, or near enough to, a fixed point; failure to do
so will result in convergence error once you get to the AUTO stage. Below are a few ways to
accomplish this:
1. We’ve already seen the first way to reach a stable fixed point, which is to run (I)nitialconds
→ (L)ast a few times until the ICs converge to the fixed point (refer to Fig. 4)
2. A second way is to use the built-in tool that computes fixed points. To do so, click/type
(S)ing pts→ (G)o. It will search for a numerical fixed point in the neighbourhood of the
values stored in the ICs data. Then it will ask if you want to print eigenvalues; if you click
‘Yes’ they will appear in the Terminal window. Important: make sure to click “Import”
in the ‘Equilibria’ window (I like to click it twice); this write the equilibrium point to the
ICs data, effectively replacing the need to perform (I)→ (L) as many times as you would
otherwise.
Figure 9: Using (S) → (G) to find the stable equilibrium point. Eigenvalues are shown in the terminal
window; both are negative and complex indicating that the fixed point at S = 2, I = 0.5 is a stable focus.
Clicking “Import” twice in the ‘Equilibria’ window (top left) sets the ICs to these fixed-point values.
3. Another way to converge to a fixed point involves using the phase portrait configuration
with the nullclines plotted, as in Fig. 8. The two intersections of green and orange nullclines
represent the two fixed points of the system. By clicking on (S)ing pts→ (M)ouse, you use
12
XPP/AUTO Tutorial
your mouse to click on/near the left-most intersection (the stable fixed point). The same
‘Equilibria’ window will pop up as before, as well as the option to display eigenvalues
(Fig. 10). Again, make sure to press ‘Import’ before continuing toward AUTO.
Figure 10: Using (S)→ (M) to click on and converge to the stable equilibrium point. A circle appears over
the stable fixed point when in phase portrait mode (hard to see here but it’s there).
Plotting manifolds of a saddle point can be done by following the same instructions as in step
3 above, only this time click on the saddle point, and make sure to press ‘Yes’ (or type Y) when
prompted to “Draw Invariant Sets?”; the result is shown in Fig. 11.
Figure 11: Using (S) → (M) to click on and converge to the saddle point (eigenvalues: −1 and 1). XPP
marks the saddle point with a triangle. Stable manifolds shown in cyan, unstable manifolds in greenish-
yellow.
13
XPP/AUTO Tutorial
5 Parameter bifurcations in AUTO
Now comes the time to tackle bifurcation diagram generation, using AUTO. Before starting,
make sure you:
• have all parameters set to their default values (Param window: click “Default,” then click
“Close” to close the window (note: pressing your system’s built-in close button won’t work
in XPP/AUTO, you have to use the buttons marked “Close” within the program).
• have imported the stable fixed-point values of S = 2 and I = 0.5 (when a = 1); follow
instructions in the previous section.
Side note: AUTO is actually its own program, but the creator of XPP/AUTO, Bard Ermentrout,
embedded AUTO into XPP in order to use time series simulations and bifurcation analyses in
conjunction with each other.
5.1 Single-parameter steady-state bifurcations
Let’s create a bifurcation diagram with a ∈ [0, 1] as our bifurcation parameter/range. To start:
1. Launch the AUTO window: press (F)ile → (A)uto; the AUTO window should pop up
(Fig 12).
Figure 12: Empty AUTO window.
14
XPP/AUTO Tutorial
2. Set which parameter you want as your bifurcation parameter in the (P)arameter window.
Make sure to set Par1 to the parameter you want to vary, in our case to a (see Fig. 13).
Note that, if you intend to do a two-parameter bifurcation, this is the moment for you
to set Par2 to be your second parameter of interest; if you run the continuation before
checking for this and realize after, you’ll have to restart and no one wants that (more on
two-parameter bifurcations in the last subsection).
Figure 13: Setting Par1 as your primary bifurcation parameter, and Par2 as your second parameter (for
two-parameter bifurcations).
3. Open the (N)umerics window, and set the values such that they match those shown in
Fig. 14, click ‘Ok’ or press Tab key when you’ve finished.
Figure 14: Numerics window. Do not change the values for those parameters boxed in red.
There are many AUTO numerics to set, and each affects the output of your bifurcation diagram
in different ways. Below (Table 1) is a list of what these actually are, as well as soft suggestions
15
XPP/AUTO Tutorial
for what values/ranges to set and when to target them if troubleshooting (visit http://www.
math.pitt.edu/~bard/xpp/help/xppauto.html for complete documentation).
4. Next, set the diagram axes for the bifurcation diagram. Press/type (A)xes → h(I)-lo and
set the values in the ‘AutoPlot’ window to those shown in Fig. 15, then press ‘Ok’/Tab key.
Figure 15: Setting the axes and their limits.
5. We’re ready! To run the bifurcation, press (R)un → (S)teady state. You should get the
bifurcation diagram shown in Fig. 16.
Figure 16: AUTO-generated bifurcation diagram, which reveals a transcritical bifurcation when a ≈ 0.67.
Stable stead-state branches are plotted in red, unstable steady state branches in black.
Pay attention to the output text in the Terminal window (Fig. 17). AUTO produces a label for
points along the computation every Npr points, whenever a special point (e.g. bifurcation point)
is found, and when the endpoints are attained (e.g. Par Min, Par Max, or Norm Max).
16
XPP/AUTO Tutorial
Param. Description Value?
Ntst Number of mesh intervals for discretization of periodic or-
bits; increase if periodic solution fails to converge (MX error
on a periodic branch)
≥ 50
Nmax The maximum number of points that AUTO will compute
before terminating the computation; increase if you find that
the computation terminates prematurely
Depends, ≥ 500
NPr Every Npr steps, AUTO spits out information about where
it’s at; it doesn’t really matter what value you choose, but too
small a number will make the diagram, and Terminal output,
appear cluttered with tags.
50-500 or more
Ds The initial step size and direction in the bifurcation parame-
ter that AUTO will take; this must be made negative if in-
tend to move from right to left; try decreasing its absolute
value if you get an MX error in first step
Depends, (+) or
(–)
Dsmin Smallest allowable step that AUTO takes (step size is adap-
tive); try decreasing if AUTO throws an MX error, or if it ap-
pears to be missing bifurcation points
Depends
EPSL Tolerance criterion (convergence of equation parameters); no
need to worry about what it is but keep it small
10−7 to 10−6
Dsmax Largest allowable step that AUTO takes; increase if going too
slow, but be careful (too large and it can miss bfn points)
Depends
Par Min Minimum allowable value for bfn parameter set as Par1 Depends
Par Max Maximum allowable value for bfn parameter set as Par1 Depends
Norm Min Minimum allowable value for the solution’s L2-norm 0
Norm Max Maximum allowable value for the solution’s L2-norm; in-
crease this if the variables in the ODE are large, though it
might take longer if too big (maybe 10 to 100 times the ex-
pected L2-norm, if known)
Depends
EPSU Tolerance criterion (convergence of solution components),
keep it small and same as EPSL
10−7 to 10−6
EPSS Tolerance criterion (arclength convergence); keep it small,
and about 100-1000 times the value set for EPSU and EPSS
10−5 to 10−4
Table 1: Parameters within AUTO’s numerics window. IMPORTANT: Values/ranges listed in the right-
most column are only my suggestions, please take with a grain of salt and use your judgement depending
on the specific problem you are attempting to solve.
17
XPP/AUTO Tutorial
Figure 17: Terminal output from bifurcation diagram in Fig. 16. BR tells you which branch the points are
on (there are two branches in our diagram); PT tells you the point number; TY tells you the type of point
it is (i.e. why it was labelled); LAB is the label that AUTO assigned that point, PAR(0) is the value of the
bifurcation parameter at this point; L2-NORM is the value of the L2-norm at this point (if this exceeds Norm
Max, the computation stops on that branch); and U(1) and U(2) are the steady-state values of S and I at
that point, respectively.
In Fig. 17, under the TY column, the points are labelled by their type (e.g. the kind of bifurcation).
Different type labels you could see:
• EP: Endpoint of a branch (reached Par Min, Par Max, Norm Min, or Norm Max), or the start
point of the computation.
• LP: Limit point or turning point of a branch (what you’ll see with a saddle-node bifurcation,
for example, or a saddle-node of periodics)
• TR: Torus bifurcation occurring on a periodic branch
• PD: Period doubling bifurcation occurring on a periodic branch
• UZ: User defined function (you can force a label to appear at specific values of parameters
or at fixed periods of a limit cycle, if you wanted to, with (U)sr period)
• MX: Failure to converge (see descriptions in Table 1 for tips on how to troubleshoot)
• BP: Bifurcation or branch point (e.g. the transcritical bifurcation in Fig. 16)
• HB: Hopf bifurcation point (from which you can perform continuation of periodic branches,
see next subsection)
• If a type is missing from a label, it’s just because it also assigns a label every Npr points.
For example, the point labelled 2 in Figs. 16 and 17 is assigned the type BP, meaning bifurcation
point. AUTO stopped computing along the first branch, which went down into negative I,
because the L2-norm exceeded Norm Max = 10 in our numerics setup.
18
XPP/AUTO Tutorial
5.2 Continuation of periodic solutions
We now turn our attention toward a different system, namely the 2-strain dengue virus model
that will be studied in more detail in class. The ODE file (dengue.ode) for this system is shown
below:
1 # Equations.
2 SS'= r*(1-N/K)*N-beta1*SS*IFF-beta2*SS*FI-mu*SS
3 SI'= beta2*SS*FI-beta1*sigma*SI*IFF-(mu+alpha2+gamma2)*SI
4 SR'= gamma2*SI-beta1*sigma*SR*IFF-mu*SR
5
6 IS'= beta1*SS*IFF-beta2*sigma*IS*FI-(mu+alpha1+gamma1)*IS
7 II'= beta1*sigma*SI*IFF+beta2*sigma*IS*FI-(mu+alpha1+alpha2+gamma1+gamma2)*II
8 IR'= gamma2*II+beta1*sigma*SR*IFF-(mu+alpha1+gamma1+delta)*IR
9
10 RS'= gamma1*IS-beta2*sigma*RS*FI-mu*RS
11 RI'= gamma1*II+beta2*sigma*RS*FI-(mu+alpha2+gamma2+delta)*RI
12 RR'= gamma1*IR+gamma2*RI-mu*RR
13
14 # Functions
15 N=SS+SI+SR+IS+II+IR+RS+RI+RR
16 IFF=IS+II+IR
17 FI=SI+II+RI
18
19 # Parameters (sigma and delta are the bifurcation parameters)
20 p sigma=0, delta=7.5
21 p alpha1=0.01, alpha2=0.02, beta1=0.025, beta2=0.05, gamma1=0.9
22 p gamma2=0.7, mu=0.01, r=1, K=1000
23
24 # Numerics
25 @ TOTAL=1000,DT=.01,xlo=0,xhi=1000,ylo=0,yhi=150
26 @ NPLOT=2,XP1=t,YP1=SI,XP2=t,YP2=IS
27 @ MAXSTOR=10000000
28 @ BOUNDS=100000
29 @ dsmin=1e-5,dsmax=.1,parmin=-.5,parmax=.5,autoxmin=-.5,autoxmax=.5
30 @ autoymax=.4,autoymin=-.5
31
32 # Initial conditions
33 init SS=0.5, SI=0.5, SR=0.5, IS=0.5, II=0.5
34 init IR=0.5, RS=0, RI=0, RR=0
19
XPP/AUTO Tutorial
Note: In order to quit XPP/AUTO to switch ODE files, you cannot just click your operating
system’s close button on the XPP/AUTO windows, you either:
• close the Terminal window altogether
• press Ctrl+C (Windows and Linux), or Command+. (period) on Mac, within the Terminal
window
• within the XPP window, press (F)ile→ (Q)uit, and press Yes or type (Y) when prompted if
you’re sure.
In order to see how we can produce periodic branches, first close SIR.ode. Then:
1. Use XPP Shortcut (Windows) or Terminal (Mac and Linux) to open dengue.ode with XP-
P/AUTO. Make sure that parameters sigma=0 and delta=7.5 in the Param window.
2. Run the simulation from the default ICs, using (I)nitialconds→ (G)o, then converge to the
steady state using (I)→ (L) or (S)→ (G) techniques as described previously.
3. Open AUTO by pressing (F)ile→ (A)uto. Make sure that Par1 is set to sigma and Par2 is
set to delta in the (P)arameter window.
4. Set your (N)umerics as shown in Fig. 18.
Figure 18: Numerics window for periodic branch continuation. Of note, Ntst should be larger than
before, or else the periodic continuation will fail to converge at one point along the branch. Norm Max
should also be bigger because the steady-state values go up to an order of magnitude of about 102.
5. Set your axes and bounds as shown in Fig. 19, by pressing (A)xes→ h(I)-lo.
20
XPP/AUTO Tutorial
Figure 19: Setting the axes and their limits.
6. Run the single-parameter bifurcation with respect to the first parameter sigma by press-
ing/typing (R)un→ (S)teady state. You should obtain a plot as in Fig. 20.
Figure 20: The single-parameter bifurcation with respect to sigma. The middle branch changes stability
from stable to unstable at a Hopf bifurcation point (HB, label 25), and goes from unstable back to stable at
another Hopf bifurcation point (HB, label 28). The goal is to plot the periodic branches that arise from the
Hopfs.
7. To plot the periodic branches that begin at the Hopf bifurcations, press/type (G)rab; you’ll
now be able to use your cursor to select points on the bifurcation branches. Select a point
near, but to the left of, the first Hopf bifurcation point (label 25), on the middle branch.
Press Tab on your keyboard until you reach the HB point labelled 25 (you can also use
left/right arrow keys but this is much slower), then press Enter (see Fig. 21).
21
XPP/AUTO Tutorial
Figure 21: Grab point labelled 25 indicated by red arrow.
8. Press/type (R)un; if you’ve done the previous steps properly, you should see a different
menu than before with the title “Hopf pt”, offering different options. Press/type (P)eriodic,
and the periodic branch should slowly start to fill in until you obtain an image like that in
Fig. 22.
Figure 22: Periodic branches are shown together with steady-state branches. Stable periodic branches are
plotted in green, unstable periodic branches are plotted in blue.
On my machine it took about 2 minutes for this branch to fill in (due to the large but necessary
value of Ntst). Note: because the two Hopf bifurcations meet via the periodic branch, the con-
tinuation will just loop around from one side to the other until Nmax is reached. You can force
the continuation to end after one loop by clicking the all-caps “ABORT” button at the bottom of
the left-hand column.
22
XPP/AUTO Tutorial
5.3 Two-parameter bifurcations
The last thing we will do is to produce a two-parameter bifurcation, which can be very useful
in painting a more detailed picture of your system. In a manner analogous to how a single-
parameter bifurcation tells you how the steady-state value of a system variable changes as a
function of a first parameter, a two-parameter bifurcation tells you how a bifurcation point varies
as a function of a second parameter. To pick up from where the previous subsection leaves off:
1. First, you must clear the grabbed data; remember that in the previous subsection, we had
to use our cursor and grab the Hopf bifurcation point (step 7), so we must undo this. Press
the Escape key, then execute (F)ile→ (C)lear grab in the AUTO window.
2. Let’s start by continuing the transcritical bifurcation in two-parameter space. Select (G)rab,
and press Tab until you reach the transcritical bifurcation (label 2, type: BP), then press
Enter.
3. Change to two-parameter view: select (A)xes→ (T)wo par. Set the axis bounds as shown
in Fig. 23. Now, the plot reads delta vs. sigma, two parameters of the model.
Figure 23: Setting the axes where Main Parm should be your Par1 and Secnd Parm should be your Par2.
4. Select (R)un → (T)wo par. You should see a portion of the transcritical branch in two
parameter space appear as in Fig. 24.
Note that, since Ds=0.002 in the (N)umerics menu is positive, the continuation will only progress
in the direction of increasing sigma, i.e. toward the right; hence why only a portion of the curve
is plotted in Fig. 24.
23
XPP/AUTO Tutorial
Figure 24: Portion of the 2-par bifurcation to the right of the starting point (Ds> 0.
5. Clear the grab again by first pressing Escape, then (F)ile→ (C)lear grab. Then, select (G)rab
and press Tab until you reach the same transcritical bifurcation point BP labelled 2 (you will
see the cursor meet with the end of the teal trace on the screen).
6. In the (N)umerics menu, change Ds to -0.002; this will ensure that AUTO will vary sigma
by decreasing its value, i.e. progress toward the left.
7. Select (R)un→ (T)wo par. The trace should now fill in as shown in Fig. 25.
Figure 25: Full branch of the 2-par bifurcation
8. Now, repeat the same steps except with the the Hopf bifurcation HB at label 25. When using
(G)rab in this view, you’ll have to press Tab 24 times until you reach it. After grabbing, run
the two-parameter continuation (R)un→ (T)wo par; do this twice, once with Ds negative,
24
XPP/AUTO Tutorial
and once again with Ds positive. Note: in both cases, it’s a good idea to manually “ABORT”
the computation once the trace makes it off the top of the screen or stalls at the same spot
for a while, since it will keep going for a long time if you leave it alone. You should get a
trace that looks like Fig. 26.
Figure 26: Two-parameter bifurcation with respect to delta and sigma. The Hopf bifurcation points are
shown in blue, and the transcritical bifurcation points in teal.
One final, important note: If you make a mistake and have to restart a bifurcation, you first have
to go back to the XPP window, and reset the parameter(s) to the values(s) you want to start the
bifurcation at, and re-converge to the steady state using (I)→ (L) or (S)→ (G), before going back
to the AUTO window.
This concludes the content of the tutorial. Feel free to post any and all questions you may have
to the Discussion board and I will be happy to answer them!
25

欢迎咨询51作业君