Difference between revisions of "ODEs"
(New page: Category:Matlab == 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 examp...) |
|||
Line 77: | Line 77: | ||
=== Systems of ODEs in Matlab === | === Systems of ODEs in Matlab === | ||
− | {{ | + | Consider a system of ODEs, |
+ | :<math> | ||
+ | \frac{\mathrm{d} \phi_i}{\mathrm{d}t} = \mathcal{F}( \phi_1 \ldots \phi_n ). | ||
+ | </math> | ||
+ | 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: | ||
+ | # The independent variable. In the system depicted above, this is time, <math>t</math> | ||
+ | # The vector of values of the dependent variable, <math>\phi_i</math> | ||
+ | |||
+ | In the event that you have a function that requires additional | ||
+ | parameters, you can often use an | ||
+ | [[Matlab_Functions#Anonymous_Functions|anonymous function]] to supply | ||
+ | these parameters. | ||
+ | |||
+ | ==== Example: Kinetics ==== | ||
+ | Consider the following reactions, | ||
+ | :<math> | ||
+ | \begin{align} | ||
+ | A + B &\stackrel{k_1}{\rightarrow} C \\ | ||
+ | A + C &\stackrel{k_2}{\rightarrow} D | ||
+ | \end{align} | ||
+ | </math> | ||
+ | Given <math>k_1=1</math>, <math>k_2=0.1</math>, and initial conditions | ||
+ | <math>C_A=1</math>, <math>C_B=0.6</math>, 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: | ||
+ | :<math> | ||
+ | \begin{align} | ||
+ | \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 | ||
+ | \end{align} | ||
+ | </math> | ||
+ | |||
+ | Reflecting this in Matlab, we have | ||
+ | {| border="2" cellpadding="5" cellspacing="1" | ||
+ | |<source lang="matlab"> | ||
+ | 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; | ||
+ | </source> | ||
+ | |} | ||
+ | |||
+ | 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 | ||
+ | [[Matlab_Functions#Anonymous_Functions|anonymous function]] as | ||
+ | follows: | ||
+ | |||
+ | {| border="2" cellpadding="5" cellspacing="1" | ||
+ | |<source lang="matlab"> | ||
+ | % set initial conditions | ||
+ | % A B C D | ||
+ | c0 = [ 1.0, 0.6, 0, 0 ]; | ||
+ | |||
+ | % 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'); | ||
+ | legend('A','B','C','D'); | ||
+ | </source> | ||
+ | |} |
Revision as of 15:53, 23 February 2009
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
This can easily be separated and solved analytically to find
where is the concentration at .
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*t;
|
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.
|
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,
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:
- The independent variable. In the system depicted above, this is time,
- The vector of values of the dependent variable,
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,
Given , , and initial conditions , , 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:
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 ];
% 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');
legend('A','B','C','D');
|