Page 204 - Numerical methods for chemical engineering

190 4 Initial value problems

Euler method becomes unstable for t > 2 resulting in spurious oscillations and divergence.
By contrast, the A-stable implicit Euler method remains well-behaved, though not accurate,
even for large time steps.

Time step restrictions for stiff systems

We can now see why explicit methods have difficulty with stiff systems. For ˙ x =
[ j]
Ax, Aw [ j] = λ j w , (4.159) yields the numerical response
N
[0] k [ j] 1 − ω j (1 − θ)
[k]
x = c µ w µ j = ω j =−( t)λ j (4.168)
j j
1 + θω j
j=1
For absolute stability, we must have |µ j |≤ 1 for every stable eigenvalue with Re(ω j ) > 0;
however, when θ< 1/2, we have an upper limit on the allowable time step (Figure 4.11). If
λ max is the largest eigenvalue magnitude, for absolute stability we must have
ω c = 2 for explicit Euler (4.169)
ω max = ( t)λ max ≤ ω c
Thus, the time step must be smaller than ( t) c = ω c /λ max . But, inspection of (4.158) shows
that the time required for x(t) to relax back to the steady state x s = 0 is governed by the
[0] λ j t
smallest eigenvalue λ min ,as c j (t) = c e . Thus, we shall have to continue the simulation
j
until t end ≈ λ −1 . To maintain absolute stability, the number of time steps has to be greater
min
than
λ −1 (λ max /λ min ) κ(A)
t end min
N min ≈ = = = (4.170)
( t) c ω c /λ max ω c ω c
κ(A) is the condition number of A. Therefore, for stiff systems with high condition numbers,
explicit methods are forced to take a very large number of time steps, requiring many
function evaluations and much CPU time. The restriction (4.169) is not present for A-stable
integrators.

Error rejection

There is another benefit that accrues from absolute stability. Let us start a simulation at
[0]
[ j]
an initial state slightly perturbed from the previous one at x [0] = N (c [0] + δc )w ,to
j
j
j=1
generate the numerical trajectory
N
[k] [0] [0] k [ j]
x = c + δc µ w (4.171)
j j j
j=1
As in (4.168), µ j is the growth coefficient. The difference between the original trajectory
(4.168) and the perturbed trajectory (4.171) is
N
[k] [0] k [ j]
δx = δc µ w (4.172)
j j
j=1
When all |µ j | < 1, the effect of this perturbation decays to zero. This property, error
rejection, is important and desirable in numerical simulation, as then the round-off errors