Initial Value Problems - Numerical Methods in Engineering with MATLAB - page 275

Numerical integration of stiff equations requires specialcare. The step size h
needed for stabilityis determinedbythe largest eigenvalue λ max ,evenif the terms

exp( − λ max x ) in the solutiondecay very rapidly and becomes insignificant as we move

away from the origin.
For example, consider the differentialequation 17

y + 1001 y +

=

1000 y 0 (7.16)

y , the equivalent first-order equations are Using y 1 = y and y 2 =

y 2

y =

− 1000 y 1 − 1001 y 2

In thiscase

0

− 1

=

1000 1001

The eigenvalues of are the roots of

=

− λ − 1

| − λ I | = 0

1000 1001 − λ

Expanding the determinant we get

− λ (1001 − λ ) + 1000 = 0

which has the solutions λ 1 = 1 and λ 2 = 1000. These equations are clearly stiff.Ac-

cording to Eq. (7.15) we wouldneed h 002 forEuler's method to bestable.

The Runge-Kuttamethodwould have approximately the samelimitation on the step
size.
When the problemis very stiff, the usual methodsofsolution,such as the Runge-
Kutta formulas, become impractical duetothevery small h required for stability. These
problems are best solvedwith methodsthat arespeciallydesigned for stiff equations.
Stiff problem solvers, which areoutside the scopeofthistext, have much better stabil-
ity characteristics; someofthemareevenunconditionally stable. However, the higher
degree of stability comes at a cost—the general rule isthatstability can be improved
onlybyreducing the order of the method (and thus increasing the truncation error).

< 2 /λ 2 = 0 .

EXAMPLE 7.7
(1) Show that the problem

19
4

y =− 10 y y (0) y − y (0) =− 9 = 0

17

Thisexample istaken from C.E. Pearson, Numerical Methods in Engineering and Science , van
Nostrand and Reinhold (1986).