Page 191 - Numerical methods for chemical engineering

Overview of ODE-IVP solvers in MATLAB 177

yields the new state x [k+1] directly without solving an algebraic system. For example, in the
explicit (forward) Euler method
[k]
[k]
[k]
F FE x , f (x; Θ) ≡ f x ; Θ x [k+1] = x [k] + ( t) f x ; Θ (4.119)

The explicit Euler method is not very accurate and is not often used. A more popular
single-step explicit rule is the fourth-order Runge–Kutta (RK4) method
(1)

k (1) = ( t) f x [k] k (2) = ( t) f x [k] + k /2
(2)

k (3) = ( t) f x [k] + k /2 k (4) = ( t) f x [k] + k (3)
(4.120)

x [k+1] − x [k] = 1 k (1) + 2k (2) + 2k (3) + k (4)
6
5 4
local error ∼ O( t ) global error ∼ O( t )
The RK4 method is said to be fourth-order as the global error, the net difference between
the numerical and true solutions over the course of a simulation, is proportional to the fourth
power of the time step t. For a method of order p, after some period of simulation time
t N = t 0 + N( t),
'
t N
p
x(t N ) − x(t 0 ) = ˙ x(t)dt = x [N] − x [0] + O[(N t) ] (4.121)
t 0
We do not derive (4.120) here, but rather the simpler RK2 method, yet the path to higher-
order RK methods becomes clear. We want to update ˙x = f (x) from t k to t k+1 = t k + t
by approximating (4.114). With the explicit Euler method, we neglect the time variation of
f (x(t)) to obtain the rule

x [k+1] − x [k] = ( t) f x [k] ≡ k (1) (4.122)
(1)
But, once we have computed k , we can use it to approximate the integrand of (4.114)
(1)
more accurately. We take the mid-point of the full explicit Euler step, x [k] + k /2, and
evaluate the function at this point
(1)

k (2) = ( t) f x [k] + k /2 ≈ ( t) f (x(t k + t/2)) (4.123)
We now form a linear approximation to f (x(t)) over the time step,
( t/2) − (t − t k ) (t − t k ) 2

f (x(t)) = f (x(t k )) + f (x(t k + t/2)) + O( t )
( t)/2 ( t)/2
(4.124)
(2)
which we write in terms of k (1) and k ,

( t/2) − (t − t k ) (2) (t − t k ) 2
(1)
f (x(t)) = k + k + O( t ) (4.125)
2
2
( t) /2 ( t) /2
We substitute this approximation into (4.114),
t k + t 1
' ( t/2) − (t − t k ) (t − t k )
2
x [k+1] − x [k] = k (1) + k (2) + O( t ) dt
( t) /2 ( t) /2
2
2
t k
(4.126)