Page 197 - Numerical methods for chemical engineering

Overview of ODE-IVP solvers in MATLAB 183

We supply a routine that computes the Jacobian,
function Jac = batch reactor ex calc Jac(t, x, k1, k2);
% allocate memory for the Jacobian matrix
N = length(x); Jac = zeros(N,N);
% extract state variables
cA = x(1); cB = x(2); cC = x(3); cD = x(4);
% set non zero values of the Jacobian
% row 1 – mole balance on A
*
*
Jac(1,1) = -k1 cB; Jac(1,2) = -k1 cA;
% row 2 – mole balance on B
*
*
*
*
Jac(2,1) = -k1 cB; Jac(2,2) = -k1 cA-k2 cC; Jac(2,3) = -k2 cB;
% row 3 – mole balance on C
*
*
*
*
Jac(3,1) = k1 cB; Jac(3,2) = k1 cA-k2 cC; Jac(3,3) = -k2 cB;
% row 4 - mole balance on D
*
*
Jac(4,2) = k2 cC; Jac(4,3) = k2 cB;
return;
and inform ode15s that this routine has been provided through odeset,
OPTIONS = odeset(‘Jacobian’, ‘batch reactor ex calc Jac’);
[t traj, x traj] = ode15s(@batch reactor ex calc f, ...
[0 t end], x 0, OPTIONS, k 1, k 2);
Other flags in OPTIONS allow us to modify the tolerances or to provide an Events function
that identifies when specified functions approach zero, and optionally stop the simulation
whenever this happens. Type help odeset for more information.

Example. Stiffness and the QSSA in chemical kinetics

We contrast the behavior of the explicit ode45 and implicit ode15s ODE solvers for a
simple chemical system that exhibits stiffness for certain parameter values. Consider a
reaction mechanism for A → B in which A first collides with an energetic third body M
*
to form an activated species A . This activated species then decomposes to the product B.
The mechanism

A + M → A + M r R1 = k 1 c M c A
∗
(4.145)
∗
A → B r R2 = k 2 c A ∗
yields the batch kinetics
dc A
= f 1 =−k 1 c M c A c A0 = 1
dt
dc A ∗
= f 2 = k 1 c M c A − k 2 c A ∗ c A 0 = 0 (4.146)
∗
dt
dc B
= f 3 = k 2 c A ∗ c B0 = 0
dt