Introductory Programming
Fall 2005

Evaluation 8 on Thursday


Using ode45 for systems of equations
------------------------------------

We started out using ode45 to solve a single, first-order equation:

ode45(@f, [0, 1], 0)

where the arguments are the name of the function file, the time
range, and the initial value of y.

In f.m, we wrote a function that computes the time derivative of y.

function dydt = f(t, y)
dydt = blah, blah
end


Then we figured out how to rewrite a second order equation as
a system of first order equations:

ode45(@freefall, [0, 30], [4000, 0])

The initial condition is a vector that (we understand) contains
the initial position and velocity.

function dxdt = f(t, X)
r = X(1);      % the first component is position
v = X(2);      % the second component is velocity
 
drdt = v;
dvdt = acceleration(t, r, v);
 
dxdt = [drdt; dvdt];    % pack the results in a column vector
end
 
function a = acceleration(t, r, v)
g = 9.8;      % acceleration of gravity in m/s^2
a = -g;
end



Then we solved a system of three first-order equations:

[T, M] = ode45(@lorenz, [0,30], [3, 2, 23]);

The initial condition is a vector with the values x, y and z.

function dVdt = f(t, V)
x = V(1);
y = V(2);
z = V(3);
 
sigma = 10;
b = 8/3;
r = 28;
 
dxdt = sigma * (y-x);
dydt = x * (r-z) - y;
dzdt = x*y - b*z;
 
dVdt = [dxdt; dydt; dzdt];    % pack the results in a column vector
end


The result from ode45 is a vector T that contains the time
values, and a matrix M that contains one column for each
variable x, y and z.


By generalizing from these examples, you should be able to
use ode45 to solve an arbitrary system of second-order equations:

x'' = fx(x, y, x', y')
y'' = fy(x, y, x', y')

You should write your program in a way that can handle any functions
fx and fy, but for debugging, you might want to try a simple
projectile:

x'' = 0
y'' = -g



How to proceed:

0) write a sketch of your program using the examples
   above to cut, paste, and modify chunks of code

1) start with a working program that is similar to what you want.
   freefall.m is probably the closest thing we have.

   Make a copy of the file with a different name, then run it
   and make sure that it works (and that you are running the
   program you think you are running).

2) whenever possible, make small, testable changes, and then test
   them.  Construct a scenario (functions fx and fy, and initial
   conditions) where you know what the right answer looks like.

3) run your function from the command line and make sure that
   it works before you try to run it from ode45

4) think carefully about the input and output parameters of each
   function, and make sure that when you call a function you
   provide the right kind of data (right size matrix)