M-LAB 2:
USING MATLAB TO NUMERICALLY INTEGRATE ODES
The purpose of this M-Lab is to learn the basics of numerical integration of systems of
first-order ODEs using MATLAB. Along the way we provide a tour of the very basic
functionalities of MATLAB. Our goals are to:
(i) learn the basics of coding a system of first-order ODEs and numerically inte-
grating solutions using built-in MATLAB functions,
(ii) learn how to make basic plots of the solutions using built-in MATLAB func-
tions, and
(iii) gain practice in graphically analysing solutions through complementary lenses:
- by plotting the time series of a given solution, and
- by plotting the corresponding solution trajectory in phase space.
Step 0 is to install MATLAB on your personal computer. If you do not have MATLAB
installed, please do this now by following the appropriate installation steps outlined
in the “Installing MATLAB ” document available under ED/Resources.
When you open MATLAB, the main screen that appears is a Command Window
together with a directory called the Current Folder window. You may also see a
Workspace window which will be empty. (If these windows are not visible, go to Lay-
out on your MATLAB screen and ensure that they are checked under SHOW.)
You can already begin making computations in your Command Window. For ex-
ample, typing
a = 3
and pressing the Return key defines the real variable a, which now shows up in
your Workspace. The variable a can be used for further computations.
We want to write more complicated programs involving several steps, so we will in-
stead use m-files.
1
2 M-LAB 2: USING MATLAB TO NUMERICALLY INTEGRATE ODES
Step 0.5 is to download the brusselator.m file from Ed/Resources. Place this file in
your main MATLAB directory. When you start MATLAB, this m-file should appear
in the list in your Current Folder window.
Figure 1
Double-click the m-file: an editor window should appear.
Figure 2
Read through the m-file slowly and carefully to get an idea of what it does. Even if
you do not have much experience with this kind of coding, it will soon become clear
that this simple m-file can be easily modified and adapted to a wide variety of problems!
M-LAB 2: USING MATLAB TO NUMERICALLY INTEGRATE ODES 3
Exercise 0. With brusselator.m open in the Editor window, press the Run button to
verify that the installation went well. The m-file should execute successfully and two
figures should be produced.
A note about MATLAB syntax: When writing/modifying m-files, it is standard
practice to end every line with a semicolon ;. The semicolon suppresses output in the
Command Window when lines of code are executed.
1. ODE Solvers in MATLAB
The brusselator.m file numerically integrates the Brusselator model, a canonical low-
dimensional model for chemical oscillations arising in autocatalytic reactions:
dx
dt
= a + x2y ? bx? x
dy
dt
= bx? x2y,
(1)
where x, y are nonnegative variables representing chemical concentrations and a, b de-
note real parameters. For the time being we will fix a = 1.
The brusselator.m script is divided into the following sections:
? % PREAMBLE
? % INITIALISATION
? % NUMERICAL INTEGRATION
? % FIGURE PLOTTING
? % DEFINITION OF THE ODEs
For the purposes of converting this mathematical definition into lines of code, the core
parts of the m-file are the % NUMERICAL INTEGRATION section and the %
DEFINITION OF THE ODEs section.
In % NUMERICAL INTEGRATION, the line
[t,z] = ode45(@(t,z) odefcn(t,z,a,b), tspan, z0, opts);
executes the ode45 function, which is the standard MATLAB ODE solver.1 The ode45
function makes a call to odefcn, another function which stores the definition of the
equations being integrated.
We also need to specify tspan, the time interval of integration, and the initial condition
z0 which stores (x0, y0). All of this initial data is provided in % INITIALISATION.
1The solver is based on a kind of Runge-Kutta method, one of the most commonly used families
of iterative methods to approximate solutions to ODEs.
4 M-LAB 2: USING MATLAB TO NUMERICALLY INTEGRATE ODES
The optional line
opts = odeset(’RelTol’,1e-3,’AbsTol’,1e-6);
specifies two error tolerances, RelTol and AbsTol, which are measured in slightly dif-
ferent ways. A general rule of thumb is that making these error tolerances smaller will
reduce the size of the discrete timesteps that ode45 takes, resulting in a more accurate
result. On the other hand, the integration will generally take longer to run, and the
data output can be much larger.
In % DEFINITION OF THE ODEs, the Brusselator equations (1) are actually
written in the odefcn function, which is placed at the end of the m-file:
function dzdt = odefcn(t,z,a,b)
dzdt = zeros(2,1);
dzdt(1) = a+(z(1)^2)*z(2)-b*z(1)-z(1);
dzdt(2) = b*z(1)-(z(1)^2)*z(2);
end
Notice that it is easy to modify this function to suit your purposes, for example by
writing entirely new equations; adding more equations for high-dimensional problems;
adding more parameters; including the time variable t explicitly in the equations to
handle nonautonomous problems; and so on.
1.1. Output. At the end of a run, all of the variables that are defined will now ap-
pear in the Workspace, where they can continue to be manipulated and viewed for the
purposes of further m-file/Command Window computations, plotting, etc.
The output of the integration step is a set of two objects:2
? an (nR × 1)-dimensional array t denoting discrete integration times, and
? an (nR × n)-dimensional matrix z denoting the corresponding discrete points
along the solution curve.
Here nR is a number depending on the ODEs being integrated, ODE solver chosen,
the time interval used, the error tolerances, and the initial condition chosen. The
number n is the dimension of the problem. In other words, each column of z lists the
discretisation of each of the coordinates along the solution trajectory, e.g., for a given
solution (x(t), y(t)) of the Brusselator system (1),
? t stores the data for discrete values within the range t ∈ [t0, tend],
? z(:,1) stores the data for the x(t) component at corresponding t values,
? z(:,2) stores the data for the y(t) component at corresponding t values, and
? z(1,1:2) corresponds to the initial condition (x(t0), y(t0)) and z(end,1:2)
corresponds to the (numerical approximation of the) final point (x(tend), y(tend)).
2Arrays are the basic objects that we manipulate in MATLAB, and they are the building blocks
for matrices, tensors, etc. Arrays are analogous to lists in Mathematica. The basic MATLAB
syntax for arrays can be found here: https://au.mathworks.com/help/matlab/learn_matlab/
matrices-and-arrays.html
M-LAB 2: USING MATLAB TO NUMERICALLY INTEGRATE ODES 5
1.2. Choice of ODE solver. MATLAB has several built-in solvers adapted to dif-
ferent purposes. The ode45 solver is quite efficient but will often fail to produce
meaningful results for stiff problems, which typically include problems having strong
scale separation (which are very common in biology!). The ode15s solver is more
computationally expensive to use but can often overcome difficulties associated with
stiff problems.
2. Plotting in MATLAB
The basic MATLAB command for plotting curves and points in the plane is
plot(x,y)
where x and y denote two arrays of equal size. Read the % FIGURE PLOTTING
section thoroughly to observe how the plot function is used to generate figures in our
Brusselator example.
MATLAB also allows you to modify all sorts of visualisation options; browse the ex-
amples in https://au.mathworks.com/help/matlab/ref/plot.html to get an idea
of what you can do.
2.1. Figure windows. After plots are generated within the figure windows, you can
continue formatting various features of the plot by selecting View | Property Edi-
tor in the toolbar (or simply by double-clicking within the figure). Experiment with
adjusting the fonts and ranges of the axes labels; changing the line and marker (point)
styles of the curves; etc.
You can see a list of all the objects that are generated within the figure by select-
ing View | Plot Browser. This window is useful for selecting particular objects, and
also for hiding some objects selectively for the purposes of visualisation.
Finally, you can save figures in both .fig formats (in case you would like to save
your work in MATLAB to continue adjusting the figure later) and various image for-
mats, like .png, .eps, and .pdf.
3. Exercises
Exercise 1. (Analysing time series and phase space plots.)
(a) Fixing the parameters a = 1, b = 1.7 in the brusselator.m file, numerically in-
tegrate four different initial conditions (x0, y0) (with x0, y0 > 0). With the aid
of time series and phase space plots, report your findings.
(b) Now fix a = 1 but vary b for five equally spaced values within the range
[1.7, 2.3]. Comment on any changes in the dynamics as b is varied.
6 M-LAB 2: USING MATLAB TO NUMERICALLY INTEGRATE ODES
Exercise 2. (Plotting practice.)
(a) Plot a graph of the function f(x) = x3 ? x within the range ?1 ≤ x ≤ 1.
Use at least 1000 points in your discrete approximation. (Hint: look up the
help page for the plot command in the mathworks website for examples. The
command linspace can be helpful.)
(b) If
dx/dt = f(x, y)
dy/dt = g(x, y)
is a dynamical system in the plane, the x-nullcline is the set of points such that
f(x, y) = 0 (i.e so that the rate of change dx/dt = 0), and the y-nullcline is the
set of points such that g(x, y) = 0 (i.e so that the rate of change dy/dt = 0).
Generically each of these nullcline sets can be described as a union of curves.
Fix the parameters a = 1, b = 1.7 in the brusselator.m file and numerically
integrate forward the trajectory with initial condition (0.1, 0.1).
Then produce a plot of the phase trajectory with the nullclines of (1) overlaid
in the figure. Plot the nullclines as dashed curves, with the x-nullcline(s) in
blue and the y-nullcline(s) in red. Report your observations concerning the
behavior of the trajectory relative to the locations of the nullclines.
Exercise 3. (Time intervals and accuracy)
(a) Fix the parameters a = 1, b = 1.7 and initial condition (x0, y0) = (0.1, 0.1) in
the brusselator.m file. For this parameter set, there is a unique stable fixed
point (1, 1.7) which attracts all initial conditions (x0, y0) with x0, y0 > 0.
Approximate the location of this fixed point to within six decimal places using
numerical integration. You may need to increase the time interval of integra-
tion tspan and reduce the error tolerances in RelTol and AbsTol to do this;
report the values you used.
(b) Fix the parameters a = 1, b = 1.95 and initial condition (x0, y0) = (0.1, 0.1) in
the brusselator.m file. Integrate forward.
Based on your numerical results, does the trajectory converge to a fixed point or
to a limit cycle (periodic orbit)? Experiment with increasing the time interval
and reducing the error tolerance in RelTol and AbsTol.
(c) Fix the parameter set a = 1, b = 2.2 in the brusselator.m file. Identify a
fixed point numerically. (Hints: integrate forward to get an idea of where a
fixed point might be. MATLAB supports integrating backwards in time– for