MATH/APPM 4650
Class notes, April 8

Stiff differential equations

Definition (sort of):
A linear differential equation with constant coefficients is called stiff if one solution of the homogeneous part is of the form y(t) = e^-ct for some "large" positive number c.

How large is "large"? Depends on the context.
Usually you can get "dimensionless" parameters by rescaling your problem. Look for "natural" length/time scales. For example if I were studying the motion of a river 10 meters wide and flowing at a rate of 20 meters per second, my natural length scale might be 10 meters and time scale might be 0.5 seconds; then if I'm studying the river over a distance of 2 kilometers, I would call that 200 in dimensionless units and the velocity would be 1 in dimensionless units.

If you had a c that was larger than about 10 in dimensionless units, it would often be considered stiff. Certainly c larger than 100 would be stiff.

Example: y' = -20 y + 5t, y(0) = 0.5.
Actual solution: y(t) = -1/80 + t/4 + 41/80 e^-20t.
Euler's method with h=0.1 is shown on the graph.


Graphically this happens because h is too large, and the vector field keeps overshooting.



If h were chosen smaller, it would work out fine. The following is with h = 0.05.
If h were slightly larger, it would be far worse. The following is with h = 0.125 (note the scale).

---

The problem is that in general you don't know the true solution, and you don't know how small h should be. Maybe the solution blows up because that's what the equations do? Or maybe h needs to be 0.00001 for Euler's method to work and that's computationally unreasonable? Or maybe there's a bug in your code!

If you suspect your equation might be stiff (e.g., you see a large parameter in the equation, or you know physically there's a lot of friction or resistance), it's better to try a method that doesn't force h to be small.

Suppose we try a model equation y' = λ y. To make this a little bit more general, we'll allow λ to be complex.

Plug this into your method and see what h has to be for stability.

Example: If we are using the midpoint method, then we get w_k+1 = w_k + h f(t_k+h/2, w_k + h/2 f(t_k, w_k)
= w_k + h λ (w_k + h/2 λ w_k)
= (1 + hλ + (hλ)²/2) w_k. We can write w_k+1 = Q(hλ) w_k where in this case Q(z) = 1 + z + z²/2. Our big question is whether the solution of this equation is close to the solution of the actual equation. (Note: not convergence, since we're not taking the limit as h goes to zero! We're just taking h to be small but fixed.)

Notice the solution of the difference equation is w_k = Q(z)^k w₀.
This decays to zero if and only if |Q(z)| < 1.

The actual solution decays to zero if and only if Re(λ) < 0.

If we want the difference equation to decay whenever the actual solution decays, no matter what h is, then we want |Q(z)| < 1 whenever Re(z) < 0.

Definition:
The region of absolute stability for a one-step method is
In the case of the midpoint method, we can write Graphing the region where this is less than one, we get the following: This means that whenever h is small enough that hλ is in the green region, the difference equation acts like the differential equation. When h is too large, we end up in the gray region, and we get weird behavior.

The midpoint method is thus not suitable for stiff differential equations.

Definition: A method is A-stable if the region of absolute stability contains the entire left half-plane (i.e., all complex numbers x+iy with x<0).

The midpoint method is not A-stable. In general no explicit method is A-stable. Some implicit methods are A-stable however.

---

Example: Consider the one-step Adams-Moulton method (the trapezoid method) given by
w_k+1 = w_k + h/2 [ f(t_k, w_k) + f(t_k+1, w_k+1) ].
Plug in the test equation y' = λ. Solve for w_k+1.

We get w_k+1 = (1+hλ/2)/(1-hλ) w_k. Thus the stability function is Q(z) = (1+z/2)/(1-z/2) and the region of absolute stability is determined by |Q(z)|² < 1, which is which reduces to x < 0. Thus the region of absolute stability is exactly the left half-plane, so the one-step Adams-Moulton method is A-stable.

In general, one could try to analyze multistep methods by using the test equation y' = λy. The corresponding difference equation would then have more than one root, which we could call Q₁(hλ), . . ., Q_m(hλ). We'd want all of these to have absolute value less than one.

However, it is known that this is impossible for a general multistep method. Of the methods we've studied, the only ones that are A-stable are the implicit Euler method and the implicit trapezoid method.