Software for integrating ODEs

Background

Numerically integrating ODE's is something you should have mastered in your pre-requisite 3E4 course. If you are not familiar with this topic, it is your responsibility to quickly catch up, because we will use this intensively for the rest of the course.

Here is a tutorial a wrote a few years ago when I taught 3E4 in 2010. The tutorial below is similar, but uses a reactor design example.

Example

Consider the pair of differential equations which you have or will cover in your Reactor Design course: <rst> <rst-options: 'toc' = False/> <rst-options: 'reset-figures' = False/> .. math::

\dfrac{dX}{dW} &= \dfrac{-r_A'}{F_{A0}} \\ \dfrac{dy}{dW} &= -\dfrac{\alpha}{2y}\left(1 + \varepsilon X\right)

Some terminology (recap from your pre-requisite math courses)

* The independent variable is \(W\) * The two dependent variables (variables being integrated with respect to the independent variable) are \(X\) and \(y\)

Since there are two dependent variables, we require initial conditions for each variable. In this case:

* :math:`X(0) = 0.0` * :math:`y(0) = 1.0`

We also need to specify initial **and final conditions** for the independent variable:

* :math:`W` initially is 0.0 at the reactor entrance * :math:`W` finally is 28.0 kg at the reactor exit

But there are a few other unknowns we must first specify:

* :math:`r_A' = -k' \dfrac{(1-X)y}{(1+\varepsilon X)}` * :math:`q = \dfrac{(1 + \varepsilon X)}{y}` * :math:`F_{A0} = 0.1362\,\text{mol.s}^{-1}` * :math:`k' = 0.0074\,\text{mol.s}^{-1}\text{.(kg catalyst)}^{-1}` * :math:`\alpha = 0.0367\,\text{kg}^{-1}` * :math:`\varepsilon = -0.15` </rst>

With MATLAB we create two files: a file containing the model equation, and another that drives the ODE integrator. Once you have both files on your computer, call ODEdriver from the MATLAB command line.

Model file: pbr.m

function d_depnt__d_indep = pbr(indep, depnt)
 
% Dynamic balance for the packed bed reactor (PBR); demo problem class 05C
% 
%    indep: the independent ODE variable, such as time or length
%    depnt: a vector of dependent variables
% 
%    X = depnt(1) = the conversion
%    y = depnt(2) = the pressure ratio = P/P_0 = y
% 
%    Returns d(depnt)/d(indep) = a vector of ODEs
 
% Assign some variables for convenience of notation
X = depnt(1);
y = depnt(2);

% Constant(s)
FA0   = 0.1362;  % [mol/s]
kdash = 0.0074;  % [mol/(kg catalyst . s)]
alpha = 0.0367;  % [1/kg]
eps   = -0.15;   % [-]

% Algebraic equation(s)
rAdash = -kdash * (1-X)/(1+eps*X) * y;

% Output from this ODE function must be a COLUMN vector, with n rows
n = numel(depnt);
d_depnt__d_indep = zeros(n,1);
d_depnt__d_indep(1) = -rAdash / FA0;
d_depnt__d_indep(2) = -alpha/(2*y) * (1+eps*X);

ODE integrator: ODEdriver.m

% Integrate the ODE
% -----------------
 
% The independent variable always requires an initial and final value:
indep_start = 0.0;   % kg
indep_final = 28.0;  % kg
 
% Set initial condition(s): for integrating variables (dependent variables)
X_depnt_zero = 0.0;   % i.e. X(W=0) = 0.0
y_depnt_zero = 1.0;   % i.e. y(W=0) = 1.0
 
% Integrate the ODE(s):
[indep, depnt] = ode45(@pbr, [indep_start, indep_final], [X_depnt_zero, y_depnt_zero]);
 
% Plot the results:
clf;
plot(indep, depnt(:,1), 'b')
grid('on')
hold('on')
plot(indep, depnt(:,2), 'g')
xlabel('Catalyst weight, W [kg]')
ylabel('X and y')
legend('X', 'y')