From Sutherland_wiki
Revision as of 16:53, 24 February 2009 by 00033394 (talk | contribs) (Solving ODEs in Matlab)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Solving ODEs in Matlab

To solve an ODE (or a system of ODEs) in Matlab, you need to build a function that evaluates the right-hand-side of the ODEs. For example, let's consider a simple ODE that describes first-order decay

\frac{\mathrm{d}c}{\mathrm{d}t} = -kc.

This can easily be separated and solved analytically to find

c = c_0 \exp\left(-kt\right),

where c_0 is the concentration at t=0.

To solve this in Matlab, we create a function like the following:

function rhs = decay( k, c )
% calculates the rate of decay given the
% rate constant (k) and concentration (c)
rhs = -k*c;

The ODE solvers in Matlab have the following requirements:

  • Provide a function that calculates the RHS of the ODE. This function must take two arguments: the independent variable (e.g. time), and the dependent variable(s).
  • Provide the time interval over which you desire the solution.
  • Provide the value of the dependent variables at the beginning of the time interval (the initial condition).

We can now solve our problem above using the function contained in decay.m to integrate this ODE. The Matlab code to do this is outlined below.

c0 = 10;
k = 0.1;
Define the initial concentration and rate constant.
rhsfun = @(t,c)( decay(k,c) );
Define a dummy anonymous function that can be used with Matlab's time integrator. This function simply calls through to the decay function defined above, supplying the appropriate arguments. The arguments of rhsfun are mandated by Matlab's ODE solver.
[t,c] = ode45( rhsfun, [0,50], c0 );
Solve the ODE on the interval [0,50]. The solution is stored in the vectors t and c.
plot(t,c); xlabel('t'); ylabel('c');
Plot the concentration as a function of time.

Of course, for this ODE we could solve this without a separate function for the RHS as follows:

c0 = 10;
k = 0.1;
rhs = @(t,c)( -k*c );
[time,conc] = ode45( rhs, [0,50], c0 );
plot(time,conc); xlabel('t'); ylabel('c');

Systems of ODEs in Matlab

Consider a system of ODEs,

  \frac{\mathrm{d} \phi_i}{\mathrm{d}t} = \mathcal{F}( \phi_1 \ldots \phi_n ).

To solve a system of ODEs in Matlab, you need to create a function that calculates the RHS of the system as a vector. The function that Matlab's integrator will call to evaluate the RHS must take two arguments:

  1. The independent variable. In the system depicted above, this is time, t
  2. The vector of values of the dependent variable, \phi_i

In the event that you have a function that requires additional parameters, you can often use an anonymous function to supply these parameters.

Example: Kinetics

Consider the following reactions,

A + B &\stackrel{k_1}{\rightarrow} C \\
A + C &\stackrel{k_2}{\rightarrow} D

Given k_1=1, k_2=0.1, and initial conditions C_A=1, C_B=0.6, determine the composition of the species as a function of time.

To solve this problem, we need to provide a function that calculates the RHS. We have the following system of ODEs implied by the mechanism:

\frac{\mathrm{d}C_A}{\mathrm{d}t} &= -k_1 C_A C_B -k_2 C_A C_C \\
\frac{\mathrm{d}C_B}{\mathrm{d}t} &= -k_1 C_A C_B \\
\frac{\mathrm{d}C_C}{\mathrm{d}t} &= k_1 C_A C_B - k_2 C_A C_C \\
\frac{\mathrm{d}C_D}{\mathrm{d}t} &= k_2 C_A C_C

Reflecting this in Matlab, we have

function rhs = kinetics( c )

% set rate constants
k1 = 1;
k2 = 0.1;

% calculate rates
r1 = k1*c(1)*c(2);
r2 = k2*c(1)*c(3);

% set the rhs.
rhs = zeros(4,1);
rhs(1) = -r1-r2;
rhs(2) = -r1;
rhs(3) = r1-r2;
rhs(4) = r2;

Note that our function has violated the rules for the interface required for Matlab's time integrators mentioned above. We either revise our function, or wrap it with an anonymous function as follows:

% set initial conditions
%      A    B    C  D
c0 = [ 1.0, 0.6, 0, 0 ];

tend = 100;  % ending time for the ODE solution

% build an anonymous function that can be used with Matlab's ODE solver
rhsfun = @(t,c)( kinetics(c) );

% solve the system of ODEs
[t,c] = ode45( rhsfun, [0,tend], c0 );

plot(t,c); xlabel('t'); ylabel('concentration');