I am reading through the very good Leveque book on numerical methods for ordinary and partial differential equations. In chapter 8 he shows a simple equation that is "stiff" in the ODE sense of the word.
$$
u' = \lambda(u – \cos{t}) – \sin{t}
$$
The solution is:
$$
u(t) = e^{\lambda(t – t_0)}(\eta – cos(t_0)) + \cos{t}
$$
So the ratio between $\lambda$ and the $\sin{t}$ will define the stiffness of the equation. That all makes sense. However, when I looked at the plot for this equation I could not understand what made it stiff? I always attribute stiffness with a lot of different time scales happening or a mixture of fast and slow frequencies. But the plot seems rather smooth, so why would the trapezoid rule have so much difficulty finding the solution?
That is my real question, based upon the plot below, what is making the solving stiff and causing the trapezoid method to bounce around to much–given that the curve is pretty smooth.
Best Answer
The $\sin(t)$ term acts like a driving force with frequency $1$. So indeed there are two time scales: one with speed $1$ (the driving force) and one with speed $\lambda$ (the part $u' = \lambda u$). If you want, you may augment the system with auxiliary equation $t' = 1$ to make it an autonomous system.
Another way to see the stiffness is to look at the stream lines. The way they turn almost at a right angle also indicates stiffness. The figure below plots streamlines of the field $(1, \lambda(u - \cos t) - \sin t)$
You can't tell that the problem is stiff by just looking at the solution (since it is always driven by slow processes away from the transition region), instead you need to inspect how rough is the trajectory field in the neighbourhood of the solution.
PS. I tried to reproduce the oscillation pattern of the trapezoidal method and I'm pretty sure that the solution plot and the numerical trajectory are obtained with different values of $\lambda$ ($\lambda \approx -5000$ for the method and $\lambda \sim -100$ for the solution). Below is the plot for $\Delta t = 0.2, \lambda = -30$ just to see how the trapezoidal method bounces in the ODE field.