Section 4.4.6.1: Choosing a Step Size for More than One ODE

For each ODE we have a desired local error, say at t0 = 0.

Calculate step size h to meet the desired error for each equation.
Choose the smallest of the step sizes calculated.

Example
ODEs:

\begin{displaymath}\frac{dM_1}{dt} = 4 - 0.2 \sqrt{M_1} \end{displaymath}


\begin{displaymath}\frac{dM_2}{dt} = 0.2 \sqrt{M_1} - 0.4 \sqrt{M_2} \end{displaymath}

with M1(0) = 100 kg, M2(0) = 225 kg.

Let desired local error be 1 kg for both M1 and M2.

Local errors are:

\begin{displaymath}e_1 \approx (h^2/2) \frac{d^2M_1}{dt^2} \end{displaymath}


\begin{displaymath}e_2 \approx (h^2/2) \frac{d^2M_2}{dt^2} \end{displaymath}

M1:
$d^2M_1/dt^2 = d( 4 - 0.2 \sqrt{M_1} )/dt = - 0.2 d(\sqrt{M_1})/dt = -(0.1/\sqrt{M_1}) dM_1/dt\\ = -(0.1/\sqrt{M_1}) ( 4 - 0.2 \sqrt{M_1} )$.
At t = 0 this gives
$ \vert d^2M_1/dt^2\vert = \vert(0.1/\sqrt 100) ( 4 - 0.2 \sqrt 100 ) \vert = 0.02.$
We require
$ e_1 \approx 1 $, or $ (h^2/2) \times 0.02 = 1 $.
So h = 10 s.

M2:
$d^2M_2/dt^2 = d( 0.2 \sqrt{M_1} - 0.4 \sqrt{M_2} )/dt = (0.1/\sqrt{M_1}) dM_1... ...( 4 - 0.2 \sqrt{M_1} ) - (0.2/\sqrt{M_2}) ( 0.2 \sqrt{M_1} - 0.4 \sqrt{M_2} )$.
At t = 0,
$ \vert d^2M_2/dt^2\vert = \vert (0.1/\sqrt 100) ( 4 - 0.2 \sqrt 100 ) - (0.2/... ... \sqrt 100 - 0.4 \sqrt 225 ) \vert = \vert 0.02 - 0.053333 \vert\\ = 0.033333 $.
We require
$ e_2 \approx 1 $, or $ (h^2/2) \times 0.033333 = 1 $.
So $ h = \sqrt 60 = 7.746 $ s.

Overall, step is restricted to be less than 7.746 s.

 


Two tanks with recycle: derivation of model

Steady state flow balances, no overall accumulation in tanks:

F0 + F3 = F1


F1 = F2 + F3

Adding equations gives F2 = F0 and if we say F3 = r F0 to define r, then

F1 = (1+r) F0

Component balance on solute, Mi being constant:

\begin{displaymath}M_1 \frac{dx_1}{dt} = F_0 x_0 + F_3 x_2 - F_1 x_1 \end{displaymath}


\begin{displaymath}M_2 \frac{dx_2}{dt} = F_1 x_1 - F_2 x_2 - F_3 x_2 \end{displaymath}

Express all flows in terms of F0 and r:

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


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

Introduce time constants: T1 = M1 / F0 and T2 = M2 / F0.
Rewrite equations (after division by F0 and a little rearrangement) as

\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}

 


Algorithm for NLAE solution

do i = 1,M

initialise xi(0) = xi(n).
end do

do k = 1,max_iterations

call multi_der_fn $(M,{\bf x}^{(k-1)},t_{n+1},$ f,DFDX)
do i = 1,M
calculate gi = xi(k-1) - xi(n) - h fi
do j = 1,M
calculate DGDXij = $\delta_{ij} - h $DFDXij
! $\delta_{ij} = 1 $ if i = j, or 0 otherwise
end do
end do
if |gi| <nlae_tolerance for all i = 1 to M, exit
solve the linear equation system

\begin{displaymath}DGDX \; {\bf\Delta x} = -{\bf g} \end{displaymath}

! use solve_linear
do i=1,M
calculate $x_i^{(k)} = x_i^{(k-1)} + \Delta x_i $
end do
end do

do i = 1,M

set xi(n+1) = xi(k-1)
end do




Remark:
You must supply a subroutine

multi_der_fn (x,t,f,DFDX)
which when given x and t will calculate and return f and DFDX where DFDXij = $\frac{\partial f_i}{\partial x_j}$.


Return to Section 4.4.6 Index
Return to Section 4.4 Index

Course Organiser
Last modified: Wed Aug 5 12:28:12 BST