Difference between revisions of "Software for integrating ODEs"

From Process Control: 3P4
Jump to navigation Jump to search
(Created page with "== 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 res...")
 
 
(7 intermediate revisions by the same user not shown)
Line 3: Line 3:
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.
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 [http://modelling3e4.connectmv.com/wiki/Software_tutorial/Integration_of_ODEs taught 3E4 in 2010]. The tutorial below is similar, but uses a reactor design example.
Here is a tutorial a wrote a few years ago when I [https://learnche.org/3E4/Software_tutorial/Integration_of_ODEs taught 3E4 in 2010]. The tutorial below is similar, but uses a reactor design example.
 


== Example ==
== Example ==
The example we will work with is a common modelling reaction: a liquid-based stirred tank reactor, with (initially) constant physical properties, a second order chemical reaction where species A is converted to B according to  \({\sf A} \rightarrow {\sf B} \), with reaction rate \(r = k C_{\sf A}^2\).  One can find the time-dependent mass balance for this system to be:
\[ \frac{dC_{\sf A}(t)}{dt} = \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - k C_{\sf A}^2 \]
where \(C_{{\sf A},\text{in}} = 5.5\) mol/L, we will initially assume constant volume of \(V = 100\) L and constant inlet flow of \(F(t) = 20.1 \) L/min.  The reaction rate constant is 0.15 \(\frac{\text{L}}{\text{mol}.\text{min}}\).  We must specify an initial condition for every differential equation: we will let \( C_{\sf A}(t=0) = 0.5\) mol/L.
In the code below we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 10.0\) minutes and plot the time-varying trajectory of concentration in the tank.
== MATLAB code ==
In a file called '''<tt>tank.m</tt>''':
<syntaxhighlight lang="matlab">
function d_depnt__d_indep = tank(indep, depnt)
% Dynamic balance for a CSTR
%
% C_A = y(1) = the concentration of A in the tank, mol/L
%
% Returns dy/dt = F/V*(C_{A,in} - C_A) - k*C_A^2


Consider the pair of differential equations which you have or will cover in your Reactor Design course:
%    indep: the independent ODE variable, such as time or length or the reactor
<rst>
%    depnt: a VECTOR of dependent variables
<rst-options: 'toc' = False/>
%
<rst-options: 'reset-figures' = False/>
%    Returns:
.. math::
%
%    d(depnt)
%  ---------- = a vector of ODEs
%    d(indep)
% Assign some variables for convenience of notation: one row per DEPENDENT variable
CA = depnt(1);


\dfrac{dX}{dW} &= \dfrac{-r_A'}{F_{A0}} \\ \dfrac{dy}{dW} &= -\dfrac{\alpha}{2y}\left(1 + \varepsilon X\right)
% Constant and other equations
F = 20.1;      % L/min
Some terminology (recap from your pre-requisite math courses)
CA_in = 2.5;  % mol/L
V = 100;      % L
k = 0.150;    % L/(min.mol)


* The independent variable is \(W\)
% The output from the ODE function must be a COLUMN vector, with n rows
* The two dependent variables (variables being integrated with respect to the independent variable) are \(X\) and \(y\)
% n = how many ODE's in this system?
n = numel(depnt); d_depnt__d_indep = zeros(n,1);
Since there are two dependent variables, we require initial conditions for each variable. In this case:


* :math:`X(0) = 0.0`
% Specify every element in the vector below: 1, 2, ... n
* :math:`y(0) = 1.0`
d_depnt__d_indep(1) = F/V*(CA_in - CA) - k*CA^2;
</syntaxhighlight>
We also need to specify initial **and final conditions** for the independent variable:


* :math:`W` initially is 0.0 at the reactor entrance
In a separate file (any name), for example: '''<tt>ode_driver.m</tt>''', which will "drive" the ODE solver:
* :math:`W` finally is 28.0 kg at the reactor exit
<syntaxhighlight lang="matlab">
% Integrate the ODE
But there are a few other unknowns we must first specify:
% -----------------


* :math:`r_A' = -k' \dfrac{(1-X)y}{(1+\varepsilon X)}`
% Set the time range:
* :math:`q = \dfrac{(1 + \varepsilon X)}{y}`
t_start = 0;
* :math:`F_{A0} = 0.1362\,\text{mol.s}^{-1}`
t_final = 10.0;
* :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 <tt>ODEdriver</tt> from the MATLAB command line.
% Set the initial condition(s):
CA_t_zero = 0.5;


=== Model file: <tt>pbr.m</tt>===
% Integrate the ODE(s):
[t, y] = ode45(@tank, [t_start, t_final], [CA_t_zero]);


<syntaxhighlight lang="MATLAB">
% Plot the results:
function d_depnt__d_indep = pbr(indep, depnt)
clf;
plot(t, y)
% Dynamic balance for the packed bed reactor (PBR); demo problem class 05C
grid('on')
xlabel('Time [minutes]')
ylabel('Concentration of A [mol/L]')
</syntaxhighlight>
 
A plot from this code shows the system stabilizing after about 9 minutes.
[[Image:Single-ode-MATLAB.png | 550px]]
 
 
Let's take a look at the MATLAB output:
<syntaxhighlight lang="matlab">
>> t(1:10)
ans =
        0
    0.0689
    0.1378
    0.2067
    0.2757
    0.5257
    0.7757
    1.0257
    1.2757
    1.5257
 
>> y(1:10)
ans =
    0.5000
    0.5248
    0.5490
    0.5726
    0.5956
    0.6742
    0.7452
    0.8090
    0.8662
    0.9171
</syntaxhighlight>
MATLAB places the points at unequal step sizes of time.  This is what their documentation has to say:
 
:The MATLAB ODE solvers utilize these methods by taking a step, estimating the error at this step, checking to see if the value is greater than or less than the tolerance, and altering the step size accordingly. These integration methods do not lend themselves to a fixed step size. Using an algorithm that uses a fixed step size is dangerous since you can miss points where your signal frequency is greater than the solver frequency. Using a variable step ensures that a large step size is used for low frequencies and a small step size is used for high frequencies. The ODE solvers within MATLAB are optimized for a variable step, run faster with a variable step size, and clearly the results are more accurate.
 
== Example with coupled ODEs ==
 
We will now expand our example to a system of two coupled ODEs. Still in the same reactor we are now considering the rate constant to be a function of temperature: \(k = 0.15 e^{- E_a/(RT)}\) L/(mol.min), with \(E_a = 5000\) J/(mol) and R = \(8.314 \) J/mol.K, and \(T(t)\) is the time-varying tank temperature, measured in Kelvin.  Furthermore, the reaction is known to release heat, with \(\Delta H_r = -590\) J/mol.  The time-dependent mass and energy balance for this system is now:
 
\[ \begin{align*}\frac{dC_{\sf A}(t)}{dt} &= \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - 0.15 e^{- E_a/(RT)} C_{\sf A}^2 \\ \frac{dT(t)}{dt} &= \frac{F(t)\rho C_p(T)}{V \rho C_p(T)}\bigg(T_\text{in} - T(t) \bigg) - \frac{0.15 e^{- E_a/(RT)} C_{\sf A}^2 V (\Delta H_r)}{V \rho C_p}\end{align*}\]
 
where \(C_{{\sf A},\text{in}} = 2.5\) mol/L, \(T_\text{in} = 288\) K; we will initially assume constant volume of \(V = 100\) L and constant inlet flow of \(F(t) = 20.1 \) L/min.  Also, \(C_p(T) = 4.184 - 0.002(T-273) \) J/(kg.K), the molar heat capacity of the liquid phase system, a weak function of the system temperature. The density of the liquid phase is \(\rho=1.05\) kg/L. We need two initial conditions as well:
* \( C_{\sf A}(t=0) = 0.5\) mol/L
* \(T(t=0) = 295 K \)
 
In the code below we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 45.0\) minutes and plot the time-varying trajectory of concentration in the tank.
 
The following code appears in an m-file called: '''<tt>tank_coupled.m</tt>'''
<syntaxhighlight lang="matlab">
function d_depnt__d_indep = tank_coupled(indep, depnt)
 
% Dynamic balance for a CSTR
%    C_A = y(1) = the concentration of A in the tank, [mol/L]
%    T  = y(2) = the tank temperature, [K]
%
%    Returns dy/dt = [F/V*(C_{A,in} - C_A) - k*C_A^2      ]
%                    [F/V*(T_in - T) - k*C_A^2*HR/(rho*Cp) ]
%
 
%    indep: the independent ODE variable, such as time or length or the reactor
%    depnt: a VECTOR of dependent variables
%  
%  
%    indep: the independent ODE variable, such as time or length
%    Returns:
%    depnt: a vector of dependent variables
%
%  
%    d(depnt)
%    X = depnt(1) = the conversion
%   ---------- = a vector of ODEs
%   y = depnt(2) = the pressure ratio = P/P_0 = y
%    d(indep)
%
%    Returns d(depnt)/d(indep) = a vector of ODEs
   
   
% Assign some variables for convenience of notation
% Assign some variables for convenience of notation: one row per DEPENDENT variable
X = depnt(1);
CA = depnt(1);
y = depnt(2);
T = depnt(2);


% Constant(s)
% Constants
FA0  = 0.1362;  % [mol/s]
F = 20.1;    % L/min
kdash = 0.0074; % [mol/(kg catalyst . s)]
CA_in = 2.5;  % mol/L
alpha = 0.0367;  % [1/kg]
V = 100.0;    % L
eps  = -0.15;   % [-]
k0 = 0.15;   % L/(mol.min)
Ea = 5000;    % J/mol
R = 8.314;    % J/(mol.K)
Hr = -590;    % J/mol
T_in = 288;  % K
rho = 1.050;  % kg/L
% Algebraic equations
k = k0 * exp(-Ea/(R*T));  % L/(mol.min)
Cp = 4.184 - 0.002*(T-273); % J/(kg.K)


% Algebraic equation(s)
% Output from this ODE function must be a COLUMN vector, with n rows
rAdash = -kdash * (1-X)/(1+eps*X) * y;
% n = how many ODEs in this system?


% Output from this ODE function must be a COLUMN vector, with n rows
n = numel(depnt); d_depnt__d_indep = zeros(n,1);
n = numel(depnt);
 
d_depnt__d_indep = zeros(n,1);
d_depnt__d_indep(1) = F/V*(CA_in - CA) - k*CA^2;
d_depnt__d_indep(1) = -rAdash / FA0;
d_depnt__d_indep(2) = F/V*(T_in - T) - (Hr*k*CA^2)/(rho*Cp);
d_depnt__d_indep(2) = -alpha/(2*y) * (1+eps*X);
</syntaxhighlight>
</syntaxhighlight>
 
In a separate file (any name), for example: '''<tt>ode_driver.m</tt>''', which will "drive" the ODE solver:
=== ODE integrator: <tt>ODEdriver.m</tt> ===
<syntaxhighlight lang="matlab">
<syntaxhighlight lang="MATLAB">
% Integrate the ODE
% Integrate the ODE
% -----------------
% -----------------
 
% The independent variable always requires an initial and final value:
% Set the time range:
indep_start = 0.0;   % kg
t_start = 0;
indep_final = 28.0; % kg
t_final = 45.0;
 
% Set initial condition(s): for integrating variables (dependent variables)
% Set the initial condition(s):
X_depnt_zero = 0.0;   % i.e. X(W=0) = 0.0
CA_t_zero = 0.5;
y_depnt_zero = 1.0;   % i.e. y(W=0) = 1.0
T_t_zero = 295;
 
% Integrate the ODE(s):
% Integrate the ODE(s):
[indep, depnt] = ode45(@pbr, [indep_start, indep_final], [X_depnt_zero, y_depnt_zero]);
[t, y] = ode45(@tank_coupled, [t_start, t_final], ...
                              [CA_t_zero, T_t_zero]);
 
% Plot the results:
% Plot the results:
clf;
clf;
plot(indep, depnt(:,1), 'b')
subplot(2, 1, 1)
plot(t, y(:,1))
xlabel('Time [minutes]')
ylabel('Concentration of A [mol/L]')
grid('on')
grid('on')
hold('on')
 
plot(indep, depnt(:,2), 'g')
subplot(2, 1, 2)
xlabel('Catalyst weight, W [kg]')
plot(t, y(:,2), 'r')
ylabel('X and y')
xlabel('Time [minutes]')
legend('X', 'y')
ylabel('Temperature in tank [K]')
grid('on')
 
print('-dpng', 'coupled-ode-MATLAB.png')
</syntaxhighlight>
The trajectories produced by the above code are shown below.  Do they make engineering sense?
[[Image:Coupled-ode-MATLAB.png | 550px]]
 
 
;'''Important notes'''
 
:* Each trajectory (ODE) must have an initial condition.
:* The function that specifies the ODE must return a ''column'' vector with \(n\) entries, one for each ODE.
:* Algebraic equations may be included in the ODE function.  For example:
:** \(C_p = 4.184 - 0.002(T-273)\) is an algebraic function
:** If the inlet flow varied over time according to \(F(t) = 20.1 + 2\sin(t)\) then this algebraic equation could be added.
:* Numerical integration allows for discontinuities. For example, a step input at \(t=50\) minutes in the inlet flow can be coded as
::
<syntaxhighlight lang="matlab">
if t < 50
    F = 20.1
else
    F = 25.1
end
</syntaxhighlight>
</syntaxhighlight>

Latest revision as of 15:04, 16 September 2018

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

The example we will work with is a common modelling reaction: a liquid-based stirred tank reactor, with (initially) constant physical properties, a second order chemical reaction where species A is converted to B according to \({\sf A} \rightarrow {\sf B} \), with reaction rate \(r = k C_{\sf A}^2\). One can find the time-dependent mass balance for this system to be:

\[ \frac{dC_{\sf A}(t)}{dt} = \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - k C_{\sf A}^2 \]

where \(C_{{\sf A},\text{in}} = 5.5\) mol/L, we will initially assume constant volume of \(V = 100\) L and constant inlet flow of \(F(t) = 20.1 \) L/min. The reaction rate constant is 0.15 \(\frac{\text{L}}{\text{mol}.\text{min}}\). We must specify an initial condition for every differential equation: we will let \( C_{\sf A}(t=0) = 0.5\) mol/L.

In the code below we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 10.0\) minutes and plot the time-varying trajectory of concentration in the tank.

MATLAB code

In a file called tank.m:

function d_depnt__d_indep = tank(indep, depnt)

% Dynamic balance for a CSTR
%
% C_A = y(1) = the concentration of A in the tank, mol/L
%
% Returns dy/dt = F/V*(C_{A,in} - C_A) - k*C_A^2

%    indep: the independent ODE variable, such as time or length or the reactor
%    depnt: a VECTOR of dependent variables
% 
%    Returns:
%
%    d(depnt)
%   ---------- = a vector of ODEs
%    d(indep)
 
% Assign some variables for convenience of notation: one row per DEPENDENT variable
CA = depnt(1);

% Constant and other equations
F = 20.1;      % L/min
CA_in = 2.5;   % mol/L
V = 100;       % L
k = 0.150;     % L/(min.mol)

% The output from the ODE function must be a COLUMN vector, with n rows
% n = how many ODE's in this system?
n = numel(depnt); d_depnt__d_indep = zeros(n,1);

% Specify every element in the vector below: 1, 2, ... n
d_depnt__d_indep(1) = F/V*(CA_in - CA) - k*CA^2;

In a separate file (any name), for example: ode_driver.m, which will "drive" the ODE solver:

% Integrate the ODE
% -----------------

% Set the time range:
t_start = 0;
t_final = 10.0;

% Set the initial condition(s):
CA_t_zero = 0.5;

% Integrate the ODE(s):
[t, y] = ode45(@tank, [t_start, t_final], [CA_t_zero]);

% Plot the results:
clf;
plot(t, y)
grid('on')
xlabel('Time [minutes]')
ylabel('Concentration of A [mol/L]')

A plot from this code shows the system stabilizing after about 9 minutes. Single-ode-MATLAB.png


Let's take a look at the MATLAB output:

>> t(1:10)
ans =
         0
    0.0689
    0.1378
    0.2067
    0.2757
    0.5257
    0.7757
    1.0257
    1.2757
    1.5257

>> y(1:10)
ans =
    0.5000
    0.5248
    0.5490
    0.5726
    0.5956
    0.6742
    0.7452
    0.8090
    0.8662
    0.9171

MATLAB places the points at unequal step sizes of time. This is what their documentation has to say:

The MATLAB ODE solvers utilize these methods by taking a step, estimating the error at this step, checking to see if the value is greater than or less than the tolerance, and altering the step size accordingly. These integration methods do not lend themselves to a fixed step size. Using an algorithm that uses a fixed step size is dangerous since you can miss points where your signal frequency is greater than the solver frequency. Using a variable step ensures that a large step size is used for low frequencies and a small step size is used for high frequencies. The ODE solvers within MATLAB are optimized for a variable step, run faster with a variable step size, and clearly the results are more accurate.

Example with coupled ODEs

We will now expand our example to a system of two coupled ODEs. Still in the same reactor we are now considering the rate constant to be a function of temperature: \(k = 0.15 e^{- E_a/(RT)}\) L/(mol.min), with \(E_a = 5000\) J/(mol) and R = \(8.314 \) J/mol.K, and \(T(t)\) is the time-varying tank temperature, measured in Kelvin. Furthermore, the reaction is known to release heat, with \(\Delta H_r = -590\) J/mol. The time-dependent mass and energy balance for this system is now:

\[ \begin{align*}\frac{dC_{\sf A}(t)}{dt} &= \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - 0.15 e^{- E_a/(RT)} C_{\sf A}^2 \\ \frac{dT(t)}{dt} &= \frac{F(t)\rho C_p(T)}{V \rho C_p(T)}\bigg(T_\text{in} - T(t) \bigg) - \frac{0.15 e^{- E_a/(RT)} C_{\sf A}^2 V (\Delta H_r)}{V \rho C_p}\end{align*}\]

where \(C_{{\sf A},\text{in}} = 2.5\) mol/L, \(T_\text{in} = 288\) K; we will initially assume constant volume of \(V = 100\) L and constant inlet flow of \(F(t) = 20.1 \) L/min. Also, \(C_p(T) = 4.184 - 0.002(T-273) \) J/(kg.K), the molar heat capacity of the liquid phase system, a weak function of the system temperature. The density of the liquid phase is \(\rho=1.05\) kg/L. We need two initial conditions as well:

  • \( C_{\sf A}(t=0) = 0.5\) mol/L
  • \(T(t=0) = 295 K \)

In the code below we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 45.0\) minutes and plot the time-varying trajectory of concentration in the tank.

The following code appears in an m-file called: tank_coupled.m

function d_depnt__d_indep = tank_coupled(indep, depnt)

% Dynamic balance for a CSTR
%    C_A = y(1) = the concentration of A in the tank, [mol/L]
%    T   = y(2) = the tank temperature, [K]
%
%    Returns dy/dt = [F/V*(C_{A,in} - C_A) - k*C_A^2       ]
%                    [F/V*(T_in - T) - k*C_A^2*HR/(rho*Cp) ]
%

%    indep: the independent ODE variable, such as time or length or the reactor
%    depnt: a VECTOR of dependent variables
% 
%    Returns:
%
%    d(depnt)
%   ---------- = a vector of ODEs
%    d(indep)
 
% Assign some variables for convenience of notation: one row per DEPENDENT variable
CA = depnt(1);
T = depnt(2);

% Constants
F = 20.1;     % L/min
CA_in = 2.5;  % mol/L
V = 100.0;    % L
k0 = 0.15;    % L/(mol.min)
Ea = 5000;    % J/mol
R = 8.314;    % J/(mol.K)
Hr = -590;    % J/mol
T_in = 288;   % K
rho = 1.050;  % kg/L
% Algebraic equations
k = k0 * exp(-Ea/(R*T));  % L/(mol.min)
Cp = 4.184 - 0.002*(T-273);  % J/(kg.K)

% Output from this ODE function must be a COLUMN vector, with n rows
% n = how many ODEs in this system?

n = numel(depnt); d_depnt__d_indep = zeros(n,1);

d_depnt__d_indep(1) = F/V*(CA_in - CA) - k*CA^2;
d_depnt__d_indep(2) = F/V*(T_in - T) - (Hr*k*CA^2)/(rho*Cp);

In a separate file (any name), for example: ode_driver.m, which will "drive" the ODE solver:

% Integrate the ODE
% -----------------

% Set the time range:
t_start = 0;
t_final = 45.0;

% Set the initial condition(s):
CA_t_zero = 0.5;
T_t_zero = 295;

% Integrate the ODE(s):
[t, y] = ode45(@tank_coupled, [t_start, t_final], ...
                              [CA_t_zero, T_t_zero]);

% Plot the results:
clf;
subplot(2, 1, 1)
plot(t, y(:,1))
xlabel('Time [minutes]')
ylabel('Concentration of A [mol/L]')
grid('on')

subplot(2, 1, 2)
plot(t, y(:,2), 'r')
xlabel('Time [minutes]')
ylabel('Temperature in tank [K]')
grid('on')

print('-dpng', 'coupled-ode-MATLAB.png')

The trajectories produced by the above code are shown below. Do they make engineering sense? Coupled-ode-MATLAB.png


Important notes
  • Each trajectory (ODE) must have an initial condition.
  • The function that specifies the ODE must return a column vector with \(n\) entries, one for each ODE.
  • Algebraic equations may be included in the ODE function. For example:
    • \(C_p = 4.184 - 0.002(T-273)\) is an algebraic function
    • If the inlet flow varied over time according to \(F(t) = 20.1 + 2\sin(t)\) then this algebraic equation could be added.
  • Numerical integration allows for discontinuities. For example, a step input at \(t=50\) minutes in the inlet flow can be coded as
if t < 50
    F = 20.1
else
    F = 25.1
end