Section 4.4.3: Numerical Solution of an ODE

Finite differences

The key to any numerical method for ODEs is to make a finite difference approximation to the derivative.

m₀ = slope of tangent at t₀
m₁ = slope of tangent at t₁
m = slope of secant from t₀ to t₁

$\begin{eqnarray*}m_{0} = \frac{dx}{dt} \vert _{t=t_{0}} \\ m_{1} = \frac{dx}{dt} \vert _{t=t_{1}} \\ m = \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \end{eqnarray*}$

m is well approximated by m₀ or m₁ if

curvature of x(t) slight over the time step t₀ to t₁;
no discontinuities or sharp kinks in x(t).

Forward difference

$\begin{eqnarray*}\frac{dx}{dt} \vert _{t=t_{0}} \approx \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \\ x_{1} \approx x_{0} + (t_{1} - t_{0}) m_{0} \end{eqnarray*}$

Backward difference

$\begin{eqnarray*}\frac{dx}{dt} \vert _{t=t_{1}} \approx \frac {x_{1} - x_{0}}{t_{1} - t_{0}} \\ x_{1} \approx x_{0} + (t_{1} - t_{0}) m_{1} \end{eqnarray*}$

Explicit Euler method

Forward difference $\rightarrow$ Explicit Euler method

(Backward difference $\rightarrow$ Implicit Euler - see later)

Solve

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

Specify x₀ at initial time t₀.

Calculate f(x) at t₀ $\rightarrow$ f₀

Finite difference formula:

x₁ = x₀ + h f₀

Algorithm: Euler's method

! The ODE is in general $\frac{dx}{dt} = f(x,t)$
! The right hand side function may
! depend on time t explicitly
! as well as implicitly through
! the behaviour of x.
Specify initial conditions:

time t₀
value of differential variable x₀
Specify final time t_f and step size h.
Initialise counter n = 0 to count steps.
do while t_n < t_f
! we still have some way to go

set t_n+1 = t_n + h
if t_n+1 > t_f then
! truncate the final step set t_n+1 = t_f
! for the last step only redefine h = t_f - t_n
end if
evaluate f_n = f(x_n, t_n)
calculate x_n+1 = x_n + h f_n
increment n by 1
end do
return

Numerical solution of a tank flow model

ODE: $\frac{dM}{dt} = b$ with b = 1 kg/s CORRECTION!

IC: M₀ = 500 kg at t₀ = 0 s

Final time: t_f = 100 s
Step size: h = 20 s

Solution

step number	time t (s)	holdup M (kg)
0	0	500
1	20	520
2	40	540
3	60	560
4	80	580
5	100	600

Numerical = Exact ?

The numerical solution equals the exact solution because the exact solution is linear.

This statement must be qualified because computers do not work in exact arithmetic. Real numbers are stored in a finite form.

Perhaps 500.00000 + 20.000000 is calculated as 519.999999.

Rounding error accumulates in computer arithmetic because of the finite representation of numbers.

Termination logic

A further consequence of rounding error.

The algorithm contains tests to decide whether

to truncate the last step,
to do another step.

The latter may have unexpected effects if, say,

80 + 20 = 99.99999

The algorithm will execute the loop again after reaching `100' (or as the computer thinks, 99.99999).
It will chop h to 0.00001.

Probably not serious, but output may look odd.

Another example: linear ODE

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

This does not mean that the solution is linear. Linear means that x appears in the equation to first degree only. (No x², x³.) Many applications in chemical engineering.

Exact solution is known:

x(t) = x₀ e^-kt

k is a rate constant for decay of x.
Units: (1/seconds).

Numerical solution

Parameter values:

k = 1 s^-1
IC: x₀ = 100
at: t₀ = 0 s
final time: t_f = 2 s
step size: h = 0.25 s

Result

step	time (s)	x(Euler)	x(Exact)	abs error
0	0.00	100.0000	100.0000	0.0000
1	0.25	75.0000	77.8801	2.8801
2	0.50	56.2500	60.6531	4.4031
3	0.75	42.1875	47.2367	5.0492
4	1.00	31.6406	36.7879	5.1473
5	1.25	23.7305	28.6505	4.9200
6	1.50	17.7979	22.3130	4.5151
7	1.75	13.3484	17.3774	4.0290
8	2.00	10.0113	13.5335	3.5222

Errors

Absolute error

e_abs,n = | x_n,Euler - x_n,exact |

Relative error

e_rel,n = | (x_n,Euler - x_n,exact) / x_n,exact |

step	time (s)	abs error	rel error
0	0.00	0.0000	0.0000
1	0.25	2.8801	0.0370
2	0.50	4.4031	0.0726
3	0.75	5.0492	0.1069
4	1.00	5.1473	0.1399
5	1.25	4.9200	0.1717
6	1.50	4.5151	0.2024
7	1.75	4.0290	0.2319
8	2.00	3.5222	0.2603

Summary

Finite difference approximation $\rightarrow$ Euler's method for ODEs
Easy to apply explicit Euler: If $\frac{dx}{dt} = f(x,t)$ , calculate f(x,t) at start of each step where x and t already known
If exact solution is linear, no errors
Linear ODE, exponential solution, $\rightarrow$ numerical solution has errors.

Next - Section 4.4.4: Time Constant, Accuracy and Step Size
Return to Section 4.4 Index
Return to Section 4 Index

Course Organiser

Last modified: Tue Aug 4 15:37:33 BST