Ordinary Differential Equations - Initial Value Problems (last updated 9/9/99)
The information in this tutorial is located in Section 16.6 of the MATLAB manual, The Student Edition of MATLAB, Version 5, The MATH WORKS Inc. Prentice-Hall, Englewood Cliffs, NJ 07632, 1997.
=============================================================================================================
This tutorial contains the following sections;
Solving a Set of Differential Equations Using MATLAB Routines 'ode23' and 'ode45'
Example Problem - A Non-Isothermal Batch Reactor
Example Problem Solution Using the Built-In MATLAB Routines
============================================================================================================
Single and sets of ordinary differential equations arise in a number of situations in chemical engineering. Remember that higher-order differential equations must first be reduced to a series of first-order equations before being solved by any of the numerical techniques available.
MATLAB has several built-in routines to solve these systems or you could write your own. This tutorial will describe the use of the built-in routines that should suffice for most all problems you encounter. Check the on-line help for routines to design stiff sets of equations.
============================================================================================================
Solving a Set of Differential Equations Using MATLAB Routines 'ode23' and 'ode45'
MATLAB contains several routines to solve single or sets of differential equations. The two that should work for most problems you will encounter are 'ode23 and ode45. 'ode23' uses second and third order Runge-Kutta formulas, and ode45 uses fourth and fifth-order formulas. The syntax for using both routines is identical. Anywhere you see ode23 in this tutorial you could substitute ode45.
The syntax for these routines is;
[ ind_var, dep_var] = ode23('ode_func',t_init, t_final, init_cond)
where ind_var is an array returning values
of the independent variable of the differential equations
dep_var is an array containing values of the dependent variables at the
corresponding values of ind_var
ode_func is the name of a function file containing the differential equations
t_init and t_final are the beginning and ending values of the ind_var for
the integration
init_cond is an array containing the initial values of the dependent variables.
In the function file the equations must be in the same array. Also the dependent variables must be elements of the same array. For example
de(1) = y(1)^2 / y(2)
de(2) = 25*exp(y(1)*y(2) )
============================================================================================================
Example Problem - A Non-Isothermal Batch Reactor
For a non-isothermal batch reactor, the rate at which
the reactant disappears and the temperature in the reactor are interdependent.
For the first-order reaction
A ---> P the following two differential equations describe
the system.
dCA / dt = k0 * CA * exp( - E / (R*T))
dT / dt = h*k0 * CA * exp( - E / (R*T))
where CA = the concentration of A
t = time
T = the reactor temperature
k0 = the first-order rate constant
h = the term expressing the temperature change due to reaction
E = the activation energy
R = the gas constant
These equations are solved below. Values for the constants, iinitial conditions and integration range and units are given in the MATLAB files below.
============================================================================================================
Example Problem Solution Using the Built-In MATLAB Routines
============================================================================================================
Create the following m-file, ode_ivp.
============================================================================================================
disp(' Darin Ridgway')
disp(' Chemical Engineering')
disp(' Last Updated August 25, 1998')
% This program solves a set of two differential
equations
% that describe the concentration and
temperature response
% in a non-isothermal batch reactor
for the first 150 seconds.
% It uses the built-in MATLAB routine
ode23.
% The reaction A ---> P is first-order,
and the first-order
% rate constant is described by the
Arrhenius expression.
% The set of differential equations
is located in function file
% ode_func.m
% Define the constants in the equations and make them
global
global E R k0 h
CA0 = 1.0;
% initial conc of A in mol/L
T0 = 300;
% initial temp in degrees Kelvin
E = 2500;
% the activation energy in J/mol
R = 8.314;
% the gas constant in J/(mol K)
k0 = 0.1;
% the pre-exponential constant in s^(-1)
h = -10;
% a term expressing the T rise due to reaction in (K L)/mol
% Set the initial values in the array initial_cond
initial_cond(1) = CA0;
initial_cond(2) = T0;
% Call the routine ODE23 to solve the system.
% The first column of dep_var contains the concentrations
of A.
% The second column of dep_var contains the temperatures.
% The routine ode45 would have the same syntax.
[time,dep_var] = ode23('ode_func',0,150,initial_cond);
% Determine the length of the arrays time and dep_var
npts = length(time);
% Create the output
fprintf('\n\n\n time conc A temp\n')
fprintf(' (sec) (mol/L) (K)\n\n')
for n = 1:npts
fprintf(' %5.1f %6.2e % 4.1f\n',time(n),dep_var(n,1),dep_var(n,2))
end
===========================================================================================================
Here is the function file, ode_func.m, that contains the differential equations.
===========================================================================================================
% Darin Ridgway
% Chemical Engineering
% Last Updated August 25, 1998
% This is the function file that contains
the two first-order
% differential equations for the example
problem.
% Define the functions as 'de'
% y(1) is the concentration of A
% y(2) is the temperature
function de = ode_func(time,y)
% The global variables are defined in
the main file, ode_ivp.m
global E R k0 h
% The first equation expresses the change
in concentration w.r.t time
% The second equation expresses the
change in temperature w.r.t time
de(1) = - k0*y(1)*exp(-E/(R*y(2)));
de(2) = - h*k0*y(1)*exp(-E/(R*y(2)));
===========================================================================================================
Return to the MATLAB prompt and run the m-file ode_ivp. The response in the Command Window is;
===========================================================================================================
Darin Ridgway
Chemical Engineering
Last Updated August 25, 1998
time
conc A temp
(sec) (mol/L)
(K)
0.0 1.00e+000
300.0
1.2 9.58e-001
300.4
10.5 6.77e-001
303.2
19.9 4.77e-001
305.2
29.3 3.36e-001
306.6
38.7 2.36e-001
307.6
48.0 1.66e-001
308.3
57.4 1.16e-001
308.8
66.8 8.14e-002
309.2
76.2 5.70e-002
309.4
85.5 4.00e-002
309.6
94.9 2.80e-002
309.7
104.3 1.96e-002
309.8
113.7 1.37e-002
309.9
123.0 9.62e-003
309.9
132.4 6.74e-003
309.9
141.8 4.72e-003
310.0
150.0 3.46e-003
310.0