Section 4.4.5: Numerical Instability

High accuracy needs small steps: $hk \ll 1$ But if accuracy is not important, can we take large steps?

Example:

$\begin{displaymath}\frac{dc}{dt} = 0.004 (0.1 - c) \end{displaymath}$

with c₀ = 0.099 kmol/m³ at t = 0 (nearly at steady state).
Try h = 900s = 15 min.

step	time (s)	c
0	0.0	0.099000
1	900.0	0.102600
2	1800.0	0.093240
3	2700.0	0.117576
4	3600.0	0.054302
5	4500.0	0.218814
6	5400.0	-0.208916
7	6300.0	0.903181
8	7200.0	-1.988271

Numerical instability if time step is large.

Reason for the instability:

1 - hk negative if $hk > 1 \Longrightarrow$ oscillations;
1 - hk < -1 if $hk > 2 \Longrightarrow$ oscillations grow.

Implicit Euler method

Avoid extrapolation of derivative from start of step.

Approximate dx/dt by backward difference.

$\begin{displaymath}\frac{dx}{dt} \vert _{t=t_{1}} \approx \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \end{displaymath}$

$\begin{displaymath}\frac{x_1 - x_0}{h} \approx f(x_1) \end{displaymath}$

Implicit Euler (for first step):

x₁ = x₀ + h f(x₁)

Must be treated as an equation for x₁. If f(x) is linear, solution is easy. If f(x) is non-linear, we have a NLAE; solution needs iteration.

$\begin{displaymath}\frac{dx}{dt} = - k x \end{displaymath}$

Implicit Euler:

x₁ = x₀ - h k x₁

$\begin{displaymath}x_1 = \frac{x_0}{1 + hk} \end{displaymath}$

Does this overcome the problem of instability?

Implicit Euler:

$\begin{displaymath}\frac{x_1}{x_0} = \frac{1}{1 + hk} \end{displaymath}$

Always between one and zero for positive k.
Tends to zero at large k.
`Strongly A-stable' method.

Example: large steps

Rerun

$\begin{displaymath}\frac{dc}{dt} = 0.004 (0.1 - c) \end{displaymath}$

with c₀ = 0.099 kmol/m³ at t = 0 and h = 900s.

step	time (s)	Implicit Euler	exact	error
0	0.0	0.099000	0.099000	0.000000
1	900.0	0.099783	0.099973	0.000190
2	1800.0	0.099953	0.099999	0.000047
3	2700.0	0.099990	0.100000	0.000010
4	3600.0	0.099998	0.100000	0.000002
5	4500.0	0.100000	0.100000	0.000000
6	5400.0	0.100000	0.100000	0.000000
7	6300.0	0.100000	0.100000	0.000000
8	7200.0	0.100000	0.100000	0.000000

Non-linear equation: solve a NLAE at each time step.

For

$\begin{displaymath}\frac{dx}{dt} = f(x,t) \end{displaymath}$

implicit Euler at time step n+1 is

x_n+1 - x_n = h f(x_n+1,t_n+1)

We need to solve

g(x_n+1) = x_n+1 - x_n - h f(x_n+1,t_n+1) = 0

for x_n+1.

Choice of NLAE method:

Newton
bisection
secant
direct substitution

For short time steps direct substitution may work well enough:

x^(k) = x_n + h f(x^(k-1),t_n+1)

where x^(k) is the kth iterate for solving x_n+1.
NB: x_n and t_n+1 are known and stay fixed while iterating.
For a linear ODE the limit at which direct substitution fails to converge is hk = 1.

Newton's method is recommended. Use the value at x_n as a starting guess for x_n+1, i.e.

x⁽⁰⁾ = x_n

Algorithm is similar to explicit Euler except for the point within the time step loop where we calculate x_n+1.

Algorithm: Implicit Euler method

Specify initial conditions: t₀ and x₀.
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
! THIS IS DIFFERENT:
solve
$g(x_{n+1}) = \\ x_{n+1} - x_n - h f(x_{n+1},t_{n+1}) = 0$ for x_n+1.
increment n by 1
end do
return

Newton algorithm for NLAE

Initialise x⁽⁰⁾ = x_n.
Do k = 1,max_iterations
! max_iterations is maximum number of Newton iterations at each time step

call deriv_and_fn
(x^(k-1),t_n+1, f,dfdx)
! you must supply a subroutine
! deriv_and_fn (x,t,f,dfdx) which
! when given x and t will calculate
! and return
! f = f(x,t) and dfdx = $\frac{\partial f}{\partial x}(x,t)$
calculate g = x^(k-1) - x_n - h f
if |g| <nlae_tolerance exit
calculate $x^{(k)} = x^{(k-1)} - \frac{g}{1 - h * \frac{\partial f}{\partial x}}$

end do
set x_n+1 = x^(k-1)

Remarks

nlae_tolerance is the equation-solving tolerance for converging Newton's method. This number should be very small compared with the desired integration accuracy for the dynamic simulation (2 or 3 orders of magnitude smaller).
Usually Newton's method converges quickly in 2 or 3 iterations, because the starting guess is good.
Set max_iterations to be about 10 or 20.
For each specific f(x,t) a different version of the subroutine deriv_and_fn (x,t,f,dfdx) will be needed. Given x and t, return f and dfdx.

Example

$\begin{displaymath}\frac{dM}{dt} = 4 - 0.2 \sqrt{M} \end{displaymath}$

with M₀ = 100 kg at t = 0.
Integrate to t = 500 s with h = 50 s.
NLAE tolerance = 0.01 kg.

step	time	M	iterations
0	0.0	100.000	0
1	50.0	169.719	3
2	100.0	221.043	3
3	150.0	259.846	3
4	200.0	289.654	3
5	250.0	312.794	3
6	300.0	330.890	3
7	350.0	345.117	3
8	400.0	356.346	3
9	450.0	365.235	3
10	500.0	372.280	2

Total number of iterations 29

Summary

With large steps (h k > 2) explicit Euler is unstable.

Implicit Euler uses backward difference and is always stable for problems which tend to a steady state, no matter how large the time step.

If the ODE is non-linear, then with implicit Euler we must solve a NLAE at each time step.

Newton's method will converge quickly in most cases, using x at the start of the step as a guess for the end of the step.

Next - Section 4.4.6: Systems of ODE's
Return to Section 4.4 Index
Return to Section 4 Index

Course Organiser

Last modified: Wed Aug 5 10:35:12 BST