# What's the intuition behind “stiff systems”?

One of the mysteries that my numerics class left me with are "stiff systems". Our prof didn't really manage to explain what they are. Now that the class (and the test) are over, I'm still curious.

The most concise definition about "stiff systems" I've found so far is: "There is no universally accepted definition of stiffness". So, if there is none, I don't expect one.

I'd be thankful about information on both of the following points:

• From what kind of (real-world) problems do stiff systems arise? What's so "stiff" about them?
• What's the good thing - or the bad thing - that I'd worry about when I have to numerically solve such a real-world problem?
-
A tutorial on stiffness detection in Mathematica's NDSolve makes for an interesting and useful read. –  Sasha May 3 '12 at 22:09
The question was discussed nicely at scicomp.SE; personally, I stick with the working definition "it's stiff if explicit methods perform inefficiently on it". –  Guess who it is. May 4 '12 at 10:03

The way I always thought of it is this:

Imagine you have two movable objects linked by a spring. When one object moves, it stretches the spring, which makes the other object move. Now, you can draw up a system of differential equations which describe all this, and numerically integrate to figure out what the system does. (Actually, a system this trivial might even admit a closed-form solution.)

Now imagine taking the spring and replacing it with a steel bar.

Scientifically, a steel bar is just like a spring - it can be stretched and squashed. But practically, the actual deformation of the bar is going to be minute. (Unless you're dealing with really huge forces.)

Basically, what you're saying is that a steel bar is just a spring with an insanely large spring constant. In other words, the steel bar is a "stiff spring".

Now think about how you're going to solve that numerically. With a "normal" spring, when you move one end, the other end oscillates around a bit and eventually settles down at the natural length of the spring. When you move a steel bar, what happens in the real world is that the other end instantly reaches equilibrium. What happens in the simulation is that the huge spring constant causes the system to oscillate like crazy, and you have to take tiny, tiny steps to make it converge to equilibrium properly.

In short, stiff systems are a pain to solve numerically. In fact, while reading an nVidia (?) tutorial on numerical integration, they gave the following advise for dealing with stiff systems: "Try to make it un-stiff. Failing that, use implicit methods. Good luck with that..."

This isn't a "definition", but should give you some idea of the intuition behind the term. (Unless, of course, I'm completely wrong... In which case I'm sure somebody will say so!)

-
I just found out what stiff systems are. But modeling steel bar as spring(for moderate forces) is just bad model. So I would like to know if you can get stiff problem even your model is good. Like for the steel bar you would better model it as a constraint. –  tom May 3 '13 at 20:11

I can try to explain a simplified aspect of the problem with an example. Let's say that you are trying to numerically solve the following linear equation:

$$y'(t) = Ay(t)$$

where $y\in\mathbb{R}^n$ and $A$ is $n\times n$. If the eigenvalues and eigenvectors of $A$ are $\lambda_i$ and $v_i$ (assumed distinct with $Re(\lambda_i)<0$) then the exact solution is

$$y(t) = \sum_i a_i \exp(\lambda_i t) v_i$$

for constants $a_i$. Now consider solving it with the forward Euler method:

$$y(t+\delta t) = y(t) + \delta t A y(t) = (I + \delta t A) y(t)$$

For this to converge we must have $|1 + \lambda_i \delta t|<0$ for all $i$.

In the simple case where all eigenvalues are real, let the most negative have value $\lambda_i=-\kappa$. In this case we require

$$\kappa \delta t - 1 < 1 \quad\Rightarrow\quad \delta t < \frac{2}{\kappa}$$

to prevent an oscillating instability, so there is a maximum step size. But one of the eigenvalues may have a very small negative real part, $\lambda_i = -\epsilon$, so that it takes on $O(1/\epsilon\delta t)$ time steps to converge to equilibrium. Putting this together, you can see that the number of steps required for the system to converge to equilibrium is

$$O(\kappa/\epsilon)$$

i.e. it depends on the ratio of the absolute values of the largest and smallest eigenvalues. This ratio is often referred to as the stiffness ratio and a stiff system is often taken to be one which has a high stiffness ratio, for some suitable definition of "high".

-
Why do we require $$\kappa \delta t - 1 < 1 \quad\Rightarrow\quad \delta t < \frac{2}{\kappa}$$ to prevent oscillations? –  Ambesh May 3 '13 at 19:38

Imagine that the solution of the initial value problem that you are integrating, changes very slowly, but on the other side, you have nearby solutions that change very fast.

Your method of integration is approximate and is never exactly 'on top' of your solution, so small errors tend to magnify themselves. In front of this situation, algorithms with variable step size will try to make the step size very small to avoid this.

This, of course, will make the process of integration very slow in the best case. In the worse case, small steps of integration will increase dramatically the rounding errors of the method.

Note that a solution of a differential equation may be stiff for some values while not stiff for other values. A general indicator of stiffness is a big difference in the eigenvalues of the linearized problem. Though this does not always applies it gives you an idea. This means that your problem is scaled in different ways, so it requires small step size and large step size at the same time.

-