Section 4.4.6: Systems of ODEs

Flow through two tanks $\Longrightarrow$ two ODEs

$\begin{displaymath}\frac{dM_1}{dt} = F_0 - F_1 \end{displaymath}$

$\begin{displaymath}\frac{dM_2}{dt} = F_1 - F_2 \end{displaymath}$

Assume F₀ = a, constant;
$F_1 = k_1 \sqrt{M_1}$ , $F_2 = k_2 \sqrt{M_2}$
Two ODEs in two unknowns M₁, M₂:

$\begin{displaymath}\frac{dM_1}{dt} = a - k_1 \sqrt{M_1} \end{displaymath}$

$\begin{displaymath}\frac{dM_2}{dt} = k_1 \sqrt{M_1} - k_2 \sqrt{M_2} \end{displaymath}$

Explicit Euler for M ODEs

General problem:

$\begin{displaymath}\frac{dx_1}{dt} = f_1(x_1,x_2,..,x_M,t) \end{displaymath}$

$\begin{displaymath}\frac{dx_2}{dt} = f_2(x_1,x_2,..,x_M,t) \end{displaymath}$

$\begin{displaymath}\frac{dx_M}{dt} = f_M(x_1,x_2,..,x_M,t) \end{displaymath}$

Approximate the derivative at the start of the step (t_n):

$\begin{displaymath}\frac{dx_i}{dt}\vert _{t_n} = \frac{x_i(n+1) - x_i(n)}{h} \end{displaymath}$

NEW NOTATION!

x_i(n) is the ith differential variable at time step n.

Explicit Euler:

x₁(n+1) = x₁(n) + h f₁(x₁(n),x₂(n),...,x_M(n),t_n)

x₂(n+1) = x₂(n) + h f₂(x₁(n),x₂(n),...,x_M(n),t_n)

x_M(n+1) = x_M(n) + h f_M(x₁(n),x₂(n),...,x_M(n),t_n)

Algorithm: explicit Euler

Specify initial conditions: time t₀, value of each differential variable x_i(0).
Specify final time t_f and step size h.
Initialise step counter n = 0.
do while t_n < t_f
set t_n+1 = t_n + h
if t_n+1 > t_f then
set t_n+1 = t_f
redefine h = t_f - t_n
end if
! DIFFERENT
do i=1,M
calc $f_i(n) = f_i({\bf x}(n), t_{n})$
calc x_i(n+1) = x_i(n) + h f_i(n)
end do
increment n by 1
end do
return

Vector notation: ${\bf x} = (x_1,x_2,...,x_M)$

Example: two tanks

a = 4 kg/s,
k₁ = 0.2 kg^0.5/s,
k₂ = 0.4 kg^0.5 /s.
Initial conditions:
M₁(0) = 100 kg,
M₂(0) = 225 kg.

Choose step size for initial local error of 1 kg on M₁ and M₂.
Calculation for M₁ $\Longrightarrow$ h = 10s.
Calculation for M₂ $\Longrightarrow$ h = 7.746s.

Details in section 4.4.6.1

So choose h < 7.746 s, say h = 7 s.

Results
Note minimum in M₁.
Dynamic simulation reveals interesting, possibly unexpected behaviour.

Stiff ODE systems

Large tanks and vessels

$\rightarrow$ large time constant
$\rightarrow$ want large steps.

Small volumes or fast rate processes

$\rightarrow$ small time constant
$\rightarrow$ need small steps.

Step size is limited by the shortest time constant in the process:

(a) for accuracy of the fast dynamics,
(b) for numerical stability,even if fast dynamics no longer needs to be accurate.

If a system contains widely different time constants it is said to be stiff. Usually this means a range of at least 3 orders of magnitude in the time constants.

Two tanks with recycle

Tank 1 : a small mixing device
Tank 2 : a large storage tank.

Assume total holdups in steady state. Recycle flow (F₃) is fraction r of feed F₀.
Model mole fraction of a dilute solute. Assume tanks are well stirred with mole fractions x₁ and x₂. Let T₁ and T₂ be tank residence times.

The solute mole fractions are modelled by

$\begin{displaymath}T_1 \frac{dx_1}{dt} = x_0 + r x_2 - (1 + r) x_1 \end{displaymath}$

$\begin{displaymath}T_2 \frac{dx_2}{dt} = (1 + r) (x_1 - x_2) \end{displaymath}$

Model derivation in section 4.4.6.1

Simulate with T₁ = 10s, T₂ = 2000s, r = 0.6, x₀ = 0.1 = constant.

ODEs:

$\begin{displaymath}\frac{dx_1}{dt} = 0.01 + 0.06 x_2 - 0.16 x_1 \end{displaymath}$

$\begin{displaymath}\frac{dx_2}{dt} = 0.0008 (x_1 - x_2) \end{displaymath}$

Results with explicit Euler

1. h = 1 s ( = 10% of smaller time constant)

step	time	x(1)	x(2)
0	0.0	0.00000	0.00000
1	1.0	0.01000	0.00000
2	2.0	0.01840	0.00001
3	3.0	0.02546	0.00002
4	4.0	0.03138	0.00004
5	5.0	0.03637	0.00007

10	10.0	0.05160	0.00024
20	20.0	0.06075	0.00069
40	40.0	0.06295	0.00167

100	100.0	0.06410	0.00457
200	200.0	0.06585	0.00922
300	300.0	0.06751	0.01364

800	800.0	0.07469	0.03271
900	900.0	0.07592	0.03599

2. h = 200 s ( = 10% of larger time constant)

step	time	x(1)	x(2)
0	0.0	0.000	0.000
1	200.0	2.000	0.000
2	400.0	-60.000	0.320
3	600.0	1865.840	-9.331
4	800.0	-57951.014	290.696

Numerical instability destroys x(1) and then x(2).

Implicit Euler: M simultaneous ODEs

Approximate the derivative at the end of the step (t_n+1):

$\begin{displaymath}\frac{dx_i}{dt}\vert _{t_{n+1}} = \frac{x_i(n+1) - x_i(n)}{h} \end{displaymath}$

Implicit Euler: solve a system of M NLAEs in M unknowns, x₁(n+1),..,x_M(n+1):

$\begin{displaymath}{\bf g}({\bf x}(n+1)) = {\bf0} \end{displaymath}$

where $g_i({\bf x}(n+1))$ $\; = x_i(n+1) - x_i(n) - h f_i({\bf x}(n+1), t_{n+1})$

Algorithm: implicit Euler

Specify initial conditions: time t₀, and each x_i(0).
Specify t_f and h.
Initialise step counter n = 0.
do while t_n < t_f
set t_n+1 = t_n + h
if t_n+1 > t_f then
set t_n+1 = t_f
redefine h = t_f - t_n
end if
! DIFFERENT FROM
! EXPLICIT EULER
solve simultaneous system of M NLAEs
$g_1({\bf x}(n+1),t_{n+1}) = 0$
...
$g_M({\bf x}(n+1),t_{n+1}) = 0$
increment n by 1
end do
return

Algorithm for NLAE solution insection 4.4.6.1

Example: two tanks with recycle.Linear ODEs and initial conditions as for explicit Euler.

h = 200 s

step	time	x(1)	x(2)
0	0.0	0.00000	0.00000
1	200.0	0.06381	0.00880
2	400.0	0.06875	0.01707
3	600.0	0.07163	0.02460
4	800.0	0.07421	0.03144

5	1000.0	0.07655	0.03766
10	2000.0	0.08543	0.06126
15	3000.0	0.09094	0.07592
20	4000.0	0.09437	0.08504

Implicit Euler solves the problem very effectively with a large time step.
Values are slightly behind (lower than) explicit Euler with h = 1 s.

Summary

Systems of ODEs

explicit Euler algorithm is like single ODE version.
implicit Euler needs solution of simultaneous NLAE system

Stiff ODEs

characterised by wide range of time constants
explicit methods must take steps $\approx$ shortest time constant; impractical
implicit methods highly effective even with step $\gg$ shortest time constant

Next - Section 4.4.7: Other Methods - Euler is not the Only Method
Return to Section 4.4 Index
Return to Section 4 Index

Course Organiser

Last modified: Wed Aug 5 12:33:01 BST