Page 192 - Numerical Methods for Chemical Engineering

Overview of ODE-IVP solvers in MATLAB 181

ODE solvers in MATLAB

A list of ODE solvers (and of other routines that act on functions) is returned by help
funfun, and a documentation window is opened by doc funfun. The two main rou-
tines of interest are ode45, an explicit single-step integrator, and ode15s, an implicit
multistep integrator that works well for stiff systems. We demonstrate the use of
ode45 and ode15s for a simple batch reactor with the two elementary reactions
A + B → C and C + B → D
dc A dc B dc C dc D
=−r R1 =−r R1 − r R2 = r R1 − r R2 = r R2
dt dt dt dt
r R1 = k 1 c A c B r R2 = k 2 c C c B (4.143)
c A (0) = c A0 c B (0) = c B0 c C (0) = 0 c D (0) = 0
To use ode45 and ode15s, we must provide at least a routine that returns for input values
of t and x, the value of ˙ x = f (t, x),
function f = calc f(t, x, P1, P2, . . .);

P1, P2, . . . are optional fixed parameters. For (4.143), the routine is
function f = batch reactor ex calc f(t, x, k1, k2);
f = zeros(size(x));
% extract concentrations from state vector
cA = x(1); cB = x(2); cC = x(3); cD = x(4);
% compute rates of each reaction
*
*
*
*
r1=k1 cA cB; r2 = k2 cC cB;
f(1) = -r1; % mole balance on A
f(2) = -r1 -r2; % mole balance on B
f(3) = r1 - r2; % mole balance on C
f(4) = r2; % mole balance on D
return;
For k 1 = 1, k 2 = 0.5, c A0 = 1, c B0 = 2, (4.143) is simulated until an end time t end = 10
by calling ode45 with the code

% set initial state
cA 0=1;cB 0=2;x 0 = [cA 0; cB 0; 0; 0];
k 1=1;k 2 = 0.5; % set fixed parameters
t end = 10; % set end time
% call ode45 to perform simulation
[t traj,x traj] = ode45(@batch reactor ex calc f, ...
[0 t end], x 0,[],k 1, k 2);
% plot results
figure; plot(t traj,x traj(:,1));
hold on; plot(t traj,x traj(:,2),‘--’);
plot(t traj,x traj(:,3),‘-.’); plot(t traj,x traj(:,4),‘:’);
xlabel(‘time t’); ylabel(‘concentration c j(t)’);