Figure 4.9:

Batch reactor with nth-order kinetics, r=kc_A^n, k_0=kc_{A0}^{n-1}=1, n<1.

Code for Figure 4.9

Text of the GNU GPL.

main.m

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

% Revised 7/24/2018

npts = 100;
tfin = 10;
K = 1;
orders = [-2;-1;-0.5;0;0.5;1;2];
norders = length(orders);
store = [];

for i = 1: norders
    t = linspace(0,tfin,npts)';
    n = orders(i);

    if (n == 1)
        c = exp(-K*t);
    else
        % test here for nonnegativity of c

        if (n < 1)
            tstop = 1/((-n+1)*K);
            t = linspace(0,tstop,npts)';
        end

        arg = (1 - (-n+1)*K*t);
        c = arg.^(1/(-n+1));
    end

    store = [store t c];
end

save -ascii order.dat store;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
    plot (store(:,1), store(:,2), store(:,3), store(:,4), ...
          store(:,5), store(:,6), store(:,7), store(:,8), ...
          store(:,9), store(:,10), store(:,11), store(:,12), ...
          store(:,13), store(:,14));
    axis([0 2 0 1])
    % TITLE
end % PLOTTING