[AMath 383] Using Matlab to solve systems of differential equations

[David Brian Walton <dbwalton@u.washington.edu>]

Dear class members,

I realize that a number of you are trying to analyze differential 
equations that involve a system of differential equations.

For those of you who are working with systems of two differential 
equations, you may find that the Java applets that we have links for on 
the class webpage will do an adequate job. (They are a bit difficult to 
extract the images sometimes, but do quite a lot.)

But for those of you who have more than two equations, you may want to 
consider using a computer program to help out. Matlab is one excellent 
option. If this is your choice, then you will probably want to use the 
matlab command 'ode45'. Matlab help explains how this will work. Here 
is a simple example:

For our chemostat problem, we had the following nondimensionalized 
system:

dS/dt = Sin - S - S B
dB/dt = S B - B - B P
dP/dt = B P - P

We are interested in finding numerical solutions for S(t), B(t) and 
P(t) for a specific number of times, starting at t=0 to t=10, examined 
every 0.1 time units. Also, we have an initial condition S0, B0 and P0.

Matlab requires that we have a function defined for the derivative 
returned as a vector. We will define the vector C=[S B P] so that 
C(1)=S, C(2)=B and C(3)=P. Here is the definition of the function for 
the growth rate:

--> Begin Matlab code <---
function dC = chemostat(t, C)
   % dC is a temporary variable that is where we assign the values of the
   % different derivatives. It will be a vector.

   % the extra parameter t allows the growth rate to depend on time as 
well as C.
   % it is required by the ode45 solver, even if you don't use it.

   global Sinput; % We declare this as a global parameter so that we can 
define it
	% somewhere else.

   dC = zeros(3,1);	% a column vector that matches the dimension of C

   % Now we put the 3 growth rate values into the vector dC.
   % Instead of calling them S, B and P, we must use elements of the 
vector C.
   dC(1) = Sinput - C(1) - C(1)*C(2);
   dC(2) = C(1)*C(2) - C(2) - C(2)*C(3);
   dC(3) = C(2)*C(3) - C(3);
---> End Matlab code <---

Assume this is saved in the 'm' file 'chemostat.m'

Now, in the matlab interactive mode, we do the following:

---> Begin Matlab interactive typing <---- (Ignore comments.)
global Sinput
Sinput = 2;	% We try solving with this particular value.
C0 = [1 1 1];	% This is our initial condition as a row
[t,C] = ode45('chemostat', [0 10], C0);	% Solve the system from time 0 
to 10.
% We use 'chemostat' to let the ode solver know where to look for the
% growth rate function.
plot(t,C(:,1));	% Plot C(1) as a function of time, that is, plot S(t)
plot(t,C(:,2));	% Plot C(2) as a function of time, that is, plot B(t)
plot(t,C(:,3));	% Plot C(3) as a function of time, that is, plot P(t)
---> End Matlab interactive typing <----

If you want to see how S and B relate to each other, you might try 
plotting
plot(C(:,1), C(:,2))

Or you might want to play with 3-d plot using
plot3(C(:,1),C(:,2),C(:,3))  % look at 3-d phase-space
plot3(C(:,1),C(:,2),t)  % look at how points move in 2-d portion as 
function of time

If any of you have problems working with this, let me know.

Sincerely,


D. Brian Walton
AMath 383 Instructor
Guggenheim 408C
(206) 685-9298