Section 4.4.4: Time constant, Accuracy and Step Size

Error in numerical solution related to curvature in true solution:

second derivative $\frac{d^2x}{dt^2}$

Solve $\frac{dx}{dt} = - k x $
One step, size h.

Euler: x1 = x0 (1 - hk)
Exact: $x_1 = x_0 \exp (-hk) $

Choice of time step

k : rate constant (s-1)
1/k: time constant
or characteristic time (s)
h : time step (s)
$\epsilon$ : desired local truncation error (units of x)


\begin{displaymath}x_0 (hk)^2 /2 = \epsilon \end{displaymath}

or

\begin{displaymath}h = \sqrt{(2 \epsilon / x_{0})} / k \end{displaymath}

General non linear equation: solution x(t)

\begin{displaymath}(h^{2} / 2) d^{2}x/dt^{2} = \epsilon \end{displaymath}

or

\begin{displaymath}h = \sqrt{2 \epsilon/(d^{2}x/dt^{2})} \end{displaymath}

Accumulation of error

k=1, x0 = 100

One step: h = 0.04s
Exact: $x(0.04) = 100 \exp(-0.04) \\ = 96.079$
Euler: $x(0.04) = 100 - 0.04 \times 100 \\ = 96.000$
Error = 0.079


Two steps: h = 0.02s
First step: $x(0.02) = 100 - 0.02 \times 100 \\ = 98.000$ (exact: 98.020)
Second step: $x(0.04) \\ = 98.000 - 0.02 \times 98.000 = 96.040 $
Error = 0.039


Global error

Local truncation error proportional to h2.Gives error over one step. Subsequent steps start with an error. Global error proportional to h.
Decaying exponential (linear ODE):

Local error declines as steady state is approached.
Global error also declines.
Maximum global error at $t \approx 1/k$.


Material balance for dilute species flowing through mixing tank:

V dc/dt = q c1 - q c

or

T dc/dt = c1 - c

where

T = V/q

Example: V = 0.5 m3, q = 0.002 m3/s
so T = 250 s
c1 = cf = 0.1 kmol/m3
c0 = 0 at t0 = 0
Exact solution: $ c = 0.1 ( 1 - \exp (-0.004 t)) $


Choose step:

for $\epsilon = 0.001$ kmol/m3 (1% of st.st)
$ 0.1 \times (0.004h)^{2} / 2 = 0.001 $
h = 35.355


Results

step time c Euler c exact abs error
0 0.000 0.000000 0.000000 0.000000
1 35.355 0.014142 0.013188 0.000954
2 70.710 0.026284 0.024636 0.001648
3 106.065 0.036709 0.034575 0.002134
4 141.420 0.045660 0.043203 0.002457
5 176.775 0.053344 0.050693 0.002652
6 212.130 0.059942 0.057195 0.002747
7 247.485 0.065607 0.062840 0.002767
8 282.840 0.070471 0.067741 0.002731
9 318.195 0.074647 0.071995 0.002652

Summary

For desired local accuracy we can choose a step size which depends on time constant 1/k for linear system, or alternatively some estimate of second derivative.

Errors accumulate

Global error > local error after one step.
Euler's method:
local error $\propto h^2$
global error $\propto h$
linear ODE - decaying exponential
max global error at $t \approx 1/k$

Variable step size

Approximate d2c/dt2 from previous three points.


\begin{displaymath}\frac{d^2c}{dt^2} \approx \frac{2}{h_n + h_{n-1}} [\frac{c_n - c_{n-1}}{h_n} - \frac{c_{n-1} - c_{n-2}}{h_{n-1}}] \end{displaymath}

where

hn = tn - tn-1


hn-1 = tn-1 - tn-2

We then estimate a new step size using

\begin{displaymath}h = \sqrt{2 \epsilon/(d^{2}c/dt^{2})} \end{displaymath}

step time c step size abs error
0 0.000 0.000000    
1 35.355 0.014142 35.355000 0.000954
2 70.710 0.026284 35.355000 0.001648
3 106.065 0.036709 35.355339 0.002134
4 144.222 0.046369 38.156248 0.002533
5 186.208 0.055376 41.986497 0.002857
6 231.751 0.063505 45.542669 0.003079
7 281.040 0.070700 49.289334 0.003193
8 335.044 0.077030 54.003740 0.003210
9 394.952 0.082534 59.907804 0.003135
10 462.030 0.087220 67.078139 0.002974
11 537.973 0.091102 75.943619 0.002729
12 625.322 0.094211 87.348124 0.002409
13 727.867 0.096586 102.545520 0.002025
14 851.443 0.098273 123.575870 0.001592
15 1000.000 0.099299 148.557098 0.001131


Next - Section 4.4.5: Numerical Instability
Return to Section 4.4 Index
Return to Section 4 Index

Course Organiser
Last modified: Tue Aug 4 16:47:20 BST