Stability and Stiff Equations

Perhaps the best way to motivate this concept is by looking at an example (taken from Iserles). Let

$$A=\left[\begin{array}{ccccccccc} %%1 2 3 4 5 6 7 -20 &10 & 0 &\c... ...: & & & &10 &-20 &10\\ 0 &\cdot & \cdot &\cdot &0 &10 &-20 \end{array}\right]$$

We solve

$${\bf Y}' = A {\bf Y}$$

$${\bf Y}(0) = {\bf I},$$ (57)

using AB2. To be specific, Figure (7) shows the Sup-norm of the approximate solution $y(t)$, for $M=10$, the size of the $M\times M$ matrix ${\bf A}$, for two and slightly different values of $h$, the step size. The dashed line corresponds to the approximate solution with $h=2.702703 \times 10^{-2}$ and the solid to $h=2.73972 \times 10^{-2}$. You can download the matlab code that was used in this example. Notice from the plots that the approximates are drastically different even though the difference in the values of $h$ is small.

Figure 7: Sup-norm of approximate solution $y(t_n)$ of (57) with $M=10$ using AB2. The dashed line corresponds to $h=2.702703 \times 10^{-2}$ and the solid line to $h=2.73972 \times 10^{-2}$

Why? after all, the difference in $h$ is so small. If you were to try AB of 3$^{rd}$ order you'd find it makes matters worse! If we had tried a BDF (low $p$, more on this later), the solution and the approximation would be reasonably close. $\Box$

Recall our analysis of Euler on the problem $y'=\lambda y$. One is tempted to conclude that a low-order method has poor approximating properties, for certain $\lambda, h$, as compared to a high order method.

But it is not the order of the method that caused the problem in the above example. Recall our analysis of trapezoidal scheme on $y'=\lambda y \rightarrow$ 2$^{nd}$ order method that showed the correct asymptotic behavior IRRESPECTIVE of $h$!

In summary: we need to understand the distinction between the ``order'' of the method and its stability, e.g. The trapezoidal has superior stability properties: in fact, we would find it to be stable independent of $h$! This does not mean that we can choose $h$ arbitrarily and expect the approximation and the solution to be close to each other: convergent and stable are not the same as accurate. However, any scheme which is consistent, convergent, and stable will be accurate if we take $h$ sufficiently small.

$\Box$

What's a stiff equation? No precise definition exists. Operationally,

$$\left\{\begin{array}{l} Y'=f(x,Y)\\ Y(x_0)=Y_0 \end{array}\right.$$

is ``STIFF'' if its numerical solution by some methods requires (perhaps in a portion of an interval) a significant depression of the step size in order to avoid instabilities.

Example: One way to assess qualitatively the stiffness of a system of equations is this: Take A a matrix of constant coefficients

$$\left\{\begin{array}{l} {\bf Y}'=A{\bf Y}\\ {\bf Y}(x_0)={\b... ...o between largest and smallest eigenvalue is huge!} \end{array}$$

Sometimes we see ``Stiffness-ratio'' as a way to ``quantify'' stiffness and is taken as ratio of the modules of largest to smallest eigenvalue of linearized system.

What's big? $10^3$ and above, perhaps.

Example: Kinetic Reactions have coupled systems with stiffness ratio $\approx 10^{17}$

Example: Bigbang (Einstein's General Theory) stiffness ratio of $\approx 10^{31}$.