代写program程序、代写MATLAB语言编程

日期：2022-08-20 10:07

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