Page 194 - Numerical methods for chemical engineering

180 4 Initial value problems

make an initial guess of x [k+1] followed by iterative solution of an implicit rule, is known
as a predictor/corrector method.
Let x [k+1,q] be our current estimate of x [k+1] at the qth Newton iteration, with the function
and Jacobian values f [k+1,q] and J [k+1,q] . We approximate f (x; Θ) in the vicinity of x [k+1,q]
as

f (x; Θ) ≈ f [k+1,q] + J [k+1,q] x − x [k+1,q] (4.140)

Substituting (4.140) into (4.128) for f (x [k+1] ; Θ), we have

α −1 x [k+1] + α 0 x [k] ≈ ( t)β −1 f [k+1,q] + J [k+1,q] x [k+1] − x [k+1,q] + U [k] (4.141)

This yields the following linear algebraic system for x [k+1, q+1] ≈ x [k+1] ,

A [k+1,q] [k+1,q+1] = b [k+1,q]
x
A [k+1,q] = α −1 I − ( t)β −1 J [k+1,q] (4.142)

b [k+1,q] =−α 0 x [k] + ( t)β −1 f [k+1,q] − J [k+1,q] [k+1,q] + U [k]
x
These Newton iterations are repeated until (4.128) is satisfied, yielding x [k+1] . The entire
process then is repeated for the next time step.
Clearly, an implicit method requires significantly more work per time step than an explicit
method; however, some tricks are available to reduce the computational burden. When t
is small, the changes in x from one iteration to the next are slight, as the difference between
x [k] and x [k+1] is proportional to t. Thus, A = α −1 I − ( t)β −1 J varies little and may be
held constant, at least temporarily. LU factorization allows us to solve a large fraction of
the systems of (4.142) without the need for elimination each time.

Stiffness and the choice of integration method

We might wonder why we bother to discuss implicit methods at all, if explicit methods allow
much easier and faster updates at each time step. The reason is that for many problems,
explicit methods may require, for reasons of numerical stability, the use of very small time
steps. Then, for a given overall simulation time, the number of updates performed with
an explicit method may be so much greater than with an implicit method that the implicit
method becomes favored.
Implicit methods are favored for IVPs that are stiff, in which the condition number of the
Jacobian matrix (the ratio of the largest and smallest eigenvalue moduli) is very large. Stiff
systems are by no means rare and so we must be prepared to use both explicit and implicit
methods. As a general rule, if we have no reason to expect that a problem is stiff, we first try
an explicit method. If it fails, i.e., it keeps running with no end in sight, we try an implicit
method. A more detailed treatment of stiffness and a comparison of the numerical stability
properties of explicit and implicit methods are provided following the demonstration of the
MATLAB ODE solvers.