Page 206 - Numerical methods for chemical engineering

192 4 Initial value problems

Stiff stability of BDF methods

Above, we have studied in detail the stability properties of simple single-step methods of
the class (4.150). More generally, explicit methods always have time restrictions of the form
(4.169), and the more a method makes use of past data to predict the future, the more strict
the time step restrictions become. We consider here the stability properties of implicit multi-
step BDF methods
m h
[k+1] [k−n] [k+1]
α −1 x + α n x = ( t)β −1 f (4.179)
n=0
The coefficients are computed from (4.138) by substituting the Hermite interpolation poly-
nomial (4.134) into the update formula (4.130). An efficient implementation of BDF
methods is presented below in the discussion of DAE systems. For the test equation
[ j]
˙ x = Ax, Aw [ j] = λ j w , we again expand the states in the numerical trajectory in eigen-
c w , and define the growth coefficient for each mode as
vectors of A, x [k] = N j=1 j [k] [ j]
[k+1] [k] [k]
µ j = c /c . Absolute stability requires |µ j |≤ 1. In terms of the {c }, we write
j j
[k]
[ j]
−n
the states at other times as x [k−n] = N c µ w . Substitution into (4.179) with
j=1 j j
f [k+1] = Ax [k+1] yields
( (
m h N N
[k] −n [ j] [k] [ j]
α n c µ w = ( t)β −1 A c µ j w (4.180)
j
j
j
n=−1 j=1 j=1
Using Aw [ j] = λ j w [ j] and collecting terms, yields
N m h (
[k] −n [ j]
0 = c j α n µ j − β −1 ( t)λ j µ j w (4.181)
j=1 n=−1
Equating the sum in the braces to zero, defining again ω j =−( t)λ j , and multiplying by
µ , we obtain
m h
j
m h
m h −n m h +1
0 = α n µ − β −1 ω j µ (4.182)
j j
n=−1
For absolute stability, all growth coefficients that are roots of this stability polynomial must
have |µ j |≤ 1. Figure 4.13 shows the regions of absolute stability for BDF methods of
various orders p = m h + 1 in the complex plane of ω j . As we are concerned only with
Re(λ j ) < 0, a method is A-stable (absolutely stable for all t > 0) if the entire half-plane
Re(ω j ) > 0 lies within the region of absolute stability. This is the case for p = 1 and p = 2,
but not for p ≥ 3, as the methods become unstable for complex eigenvalues.
However, we see from Figure 4.13 that even though these higher-order methods are
not A-stable, they do remain absolutely stable when Re(ω j ) =−( t)Re(λ j ) 1. Thus,
when we use a time step much larger than the characteristic time [Re(λ j )] −1 of the mode,
[k+1] [k]
|c /c | < 1 and these modes decay to zero. Thus, these methods are said to have stiff
j j
decay.
From our discussion of the discretized time-dependent diffusion equation (4.174), we
haveseenthatthefastestmodeshavespatialwavelengthscomparabletothedistancebetween