• Second 3-week block: Example 1

Systems of non-linear differential equations

In this 3 week block we are going to address a major issue in using differential equations to model physical or biological phenomena: the nature of the solution curves to a system of coupled non-linear differential equations.

These types of systems of differential equations arise in most serious modeling of physical and biological phenomena. There are very few known exact methods of solution for these systems of differential equations. As a result, people working with these types of systems of differential equations have to resort to numerical methods for obtaining approximations to solution curves.

You have already used one such method in the first week: Euler’s method. We will see, in more detail in this 3-week block, that Euler’s method does not produce an accurate enough solution for long enough values of the variables, to be a serious method for obtaining approximate numerical solutions.

More powerful methods – with smaller error -  are needed and we will focus on Runge-Kutta methods.

Edward Lorenz’s equations and the Lorenz attractor

Edward Lorenz (born in New England – West Hartford, Connecticut in 1917, and died in April 2008 in Cambridge, Massachusetts, aged 90) set up a simplified model of convection rolls arising in the equations of the atmosphere, in 1963.

Edward Lorenz

His simplified model involves 3 linked differential equations of 3 variables. Specifically, the equations are:

In these equations, all parameters are positive.

The parameters and   are related to physical properties of a fluid in motion: kinematic viscosity and heat flow (the Prandtl number and Rayleigh number, respectively). In these equations it is common to take .

Different values of give qualitatively different solution curves for this system of equations. For example, for the solution is chaotic, meaning that nearby points diverge approximately exponentially rapidly (so prediction with time is difficult to impossible). The solution curves for these values form an object in 3-space that is known as the Lorentz attractor:

The Lorenz attractor

For other values of the system has solutions that are periodic, but knotted in space.  A classic example of a knotted solution occurs when

Euler approximation to solution curves

Excel

The method for setting up an Euler approximation to a system of first-order differential equations is scarcely more difficult than for a single first-order eqation. We just have to remember that in the Euler approximation formula the variables are linked.

Just as you did for a single differential equation, we can write Excel formulas to carry out Euler’s method for the 3 linked Lorenz equations. Below is such a spreadsheet, with the formulas shown (to enlarge the image, click on it):

Excel does not have a built-in feature that allows XYZ scatter plots, so we cannot easily plot the points in three-dimensional space, as varies.

However, there are some things we can do:

• We can plot each of as functions of :

• We can plot each of
  • y(t) versus x(t):

  

  • z(t) versus x(t):

• z(t) versus y(t):

MATLAB

MATLAB comes equipped with a 3D plot capability, so the first thing for us to do if we are to understand how a 3D plot arises from an Euler approximation to the Lorenz system of differential equations, is to implement Euler’s method for a system of differential equations in MATLAB .

To do this, we can either do it “inline” following the MATLAB prompt, or we can define an M-file to run. We will do the latter, as on pp. 94-95 of the text Differential Equations With MATLAB, modified for a system of differential equations rather than a single differential equation.

A very well laid out explanation of Euler’s method in MATLAB for a system of differential equations appears on the Website Simple ODE Methods. The code below is based on the code at that Web site.

Before we write code to numerically approximate solutions to a system of differential equations, we will first write the code for a single differential equation. Because MATLAB handles vectors so easily, we will then simply modify the code for a single differential equation to a system of equations by replacing a single solution function by a vector of solution functions.

We wiil first write a MATLAB M file called euler.m (saved as “euler”),  the code for which appears below:

% Euler’s method for a single differential equation.

% The differential equation is dy/dt = f(t,y).

function [ t, y ] = euler ( f, t_range, y_initial, nstep )

% We set a range of time values, from t(1) to t(2).

t(1) = t_range(1);

% We define dt by dividing the time range into an equal number of specified steps.

dt = ( t_range(2) – t_range(1) ) / nstep;

% We set the initial value of y at the beginning time t(1).

y(1) = y_initial;

% We use Euler’s method to update the value of y at new time steps.

%  “feval” is used instead of “eval” because we are passing the name of f to the program.

for i = 1 : nstep
t(i+1) = t(i) + dt;
y(i+1) = y(i) + dt * feval ( f, t(i), y(i) );
end

% The following command plots y as a function of t.

plot(t,y)

The M File euler.m returns two data vectors (= lists):

1. t, the value of the independent variable at nstep+1 equally spaced points;
2. y, the value of the dependent variable at each t.

and plots y as a function of t over the specified time interval.

Let’s implement this for the first order differential equation for and .

For this example, .

We set up a simple M file, let’s call it test_example.m (saved as “test_example) that describes this function:

function yprime = test_example ( t, y )
yprime = y*(2./t-1);

Now to run this in the command line we set up the initial condition, the time range, and the number of steps to take, and call “test_example” as an argument to euler.m:

>> y_init = 0.1;
>> [ t, y ] = euler ( ‘test_example’, [ 0.1, 9.0 ], y_init, 200 );

Note the following things:

• we specify the initial value of y at the beginning time;
• the name of the function test_example appears in single quotes;
• the time interval is specified;
• the number of steps is specified.

After running this program we get the following graphic output:

We now want to modify euler.m so that it can calculate a numerical approximation to the Lorenz equations. We will do this in a general way by writing the code for a system of first-order equations, in an M file called euler_system.m (saved as “euler_system”):

% Euler’s method for a system of differential equations.

% The differential equations are dy/dt = f(t,y) where y is a vector of unknown functions.

function [ t, y ] = euler_system( f, t_range, y_initial, nstep )

% We set a range of time values, from t(1) to t(2).

t(1) = t_range(1);

% We define dt by dividing the time range into an equal number of specified steps.

dt = ( t_range(2) – t_range(1) ) / nstep;

% We set the initial value of the vector y at the beginning time t(1).

y(:,1) = y_initial;

% We use Euler’s method to update the value of y at new time steps.

%  “feval” is used instead of “eval” because we are passing the name of f to the program.

for i = 1 : nstep
t(i+1) = t(i) + dt;
y(:,i+1) = y(:,i) + dt * feval ( f, t(i), y(:,i) );
end

% The following command plots the components of y as functions of t.

plot(t,y)

% We also want a 3D plot of the vector components as they vary with t.

plot3(y(1,:),y(2,:),y(3,:))

There is no need to specify how many component differential equations there are -  MATLAB figures that from the form of the variables in the prescription for the function f, which is in a separate M file.

The colon  in y(:,1) = y_initial is used because y is a column vector (that’s what the “1″ tells us) with as many components (= rows) as there are unknown functions.

Now we write an M file, lorenz_system.m, which returns the derivative function as a column vector:

function yprime = lorenz_system ( t, y )
yprime = [ 10.0* (y(2)-y(1)); y(1)*(28.0-y(3))-y(2);y(1)*y(2)-8*y(3)/3 ];

To run this in the command line we set up the vector of initial condition, the time range, and the number of steps to take, and call “lorenz_system” as an argument to euler_system.m:

>> y_init = [ rand(); rand(); rand() ];

>> [ t, y ] = euler_system ( ‘lorenz_system’, [ 0.0, 20.0 ], y_init, 1000 );

After running this program we get the following graphic output:

MATLAB’s Numerical ODE solvers

MatLab has a number of built-in packages for obtaining approximate numerical solutions to differential equations and systems of differential equations. The one we use to get a plot of the Lorenz attractor is ode45.

To get background information about ode45 (and other ode solvers):

1. click on the MATLAB help link at the top of the command window;
2. enter “ode45″ in the search window and press “GO”;
3. in the Demo Results window click on the first line containing the phrase “ode45″;
4. a new window opens with the background information.

To use this for the Lorenz system we set up two M-files (ref. this Web site).

The first one, g.m (save it as “g”), is as follows:

function xdot = g(t,x)

xdot = zeros(3,1);

sig = 10.0;

rho = 28.0;

bet = 8.0/3.0;

xdot(1) = sig*(x(2)-x(1));

xdot(2) = rho*x(1)-x(2)-x(1)*x(3);

xdot(3) = x(1)*x(2)-bet*x(3);

The second, lorenz_demo.m (save it as “lorenz_demo”), is as follows:

function lorenz_demo(time)

% Usage: lorenz_demo(time)

% time=end point of time interval

% This function integrates the lorenz attractor

% from t=0 to t=time

[t,x] = ode45(‘g’,[0 time],[1;2;3]);

disp(‘press any key to continue …’)

pause

plot3(x(:,1),x(:,2),x(:,3))

print -deps lorenz.eps

Save these as M-files and then enter “lorenz_demo(200)” at the MATLAB prompt.

MATLAB’s built-in demonstration of the Lorenz attractor

MATLAB has a built-in demonstration of the dynamics of the Lorenz attractor in time. Just type “lorenz” after the MATLAB prompt.

To see the code, do the following:

1. click on the MATLAB help link at the top of the command window;
2. enter “lorenz” in the search window and press “GO”;
3. in the Demo Results window click on “Lorenz Attractor Animation”;
4. when the new window opens click on “Open lorenz.m in the Editor”.

MATLAB help

In addition to the book Differential Equations with MATLAB, you might find the following pages helpful:

• A Practical Introduction to MATLAB
• MATLAB help

Advertisement