Page 198 - Numerical Methods for Chemical Engineering

Accuracy and stability of single-step methods 187

a simulation of duration t sim , the number of time steps taken is t sim /( t). If at each of
2
these steps, we introduce an additional local error proportional to ( t) , then the global
error, the net difference between the numerical trajectory and the true trajectory x(t), is
2
proportional to [t sim /( t)]( t) = t sim ( t). As the global error is proportional to the first
power of t, (4.150) is said to be a first-order accurate integration method. We usually wish
to use integration methods of high orders of accuracy, as then a reduction in t yields a
dramatic decrease in the global error. Note that the most accurate choice of θ is θ = 1/2,
2
the Crank–Nicholson method, for which the ( t) terms nearly cancel out.

Absolute stability of an integration method

Accuracy is not the only consideration when choosing an integration method. Usually, a
property of more practical importance is numerical stability. We examine here questions of
stability for the linear test equation

˙ x = Ax (4.156)
where we assume that the matrix A is diagonalizable. We wish to perturb the system by a
small amount ε away from its steady state x s = 0 and then compare the resulting numerical
At
response with a time step t to the true response x(t) = e ε.As A is assumed to be
diagonalizable, we can write
[0]
[0]
[0]
ε = x [0] = c w [1] + c w [2] +· · · + c w [N] c [0] ∈ C
1 2 N j
(4.157)
Aw [ j] = λ j w [ j] λ j = a j + ib j
The true response of the system to this perturbation is
N
At [0] λ j t [ j]
x(t) = e ε = c e w (4.158)
j
j=1
If x s = 0 is stable, Re(λ j ) < 0 for all j, and the true response x(t) returns to x s = 0 for all
perturbations. Thus, if we find that the numerical response does not behave similarly, we
obviously have a problem.

Definition Absolute stability. Let us generate a sequence of values x [k] that are meant to
approximate the response of the system ˙ x = Ax at times t k = k( t) for k = 0, 1, 2, ....
We expand each member of the sequence as a linear combination of the eigenvectors of
[k]
[ j]
A, x [k] = N c w . We say that the rule generating this sequence is absolutely stable if,
j=1 j
[k+1] [k]
for every eigenvalue of A with Re(λ j ) < 0, |c |≤|c |. That is, the numerical response
j j
of the system does not grow along the stable manifold of A.
From (4.150), the numerical response for ˙x = Ax is generated by the rule
[k+1] [k] [k] [k+1]
x = x + ( t) (1 − θ)Ax + θ Ax (4.159)
We expand the old and new states in eigenvectors of A,
N N
[k] [ j] [k+1] [k+1] [ j]
[k]
x = c w x = c w (4.160)
j j
j=1 j=1