# Stiff differential equation where Runge-Kutta $4$th order method can be broken

Is there a stiff differential equation that cannot be solved by the Runge-Kutta 4th order method, but which has an analytical solution for testing?

How about "impractical"? Here's an example from Numerical Recipes in C, Second Ed., Sec. 16.6:

$$u'=998 u+1998 v$$

$$v'=-999 u -1999 v$$

where $u(0) = 1$ and $v(0) = 0$.

Then

$$u(x) = 2 e^{-x}- e^{-1000 x}$$

$$v(x) = -e^{-x} + e^{-1000 x}$$

Runge-Kutta 4th order would require a stepsize of less than $0.001$ for any accuracy, which is obviously aggravating for the one of the solutions which could be acquired with a far coarser step.

• Thanks, I like your equations and the book that was referenced, but if you calculate eigenvalues of the Jacobian, you will find $\lambda_1 = 2560$ and $\lambda_2 = -1559$. A system of differential equations is stiff, when $$S = \dfrac {max_{1 \leqslant i \leqslant n} \vert{Re(-\lambda_i)}\vert} {min_{1 \leqslant i \leqslant n} \vert{Re(-\lambda_i)}\vert} >> 1.$$ For this equations $S = 1.64$, please correct me if I wrong. – Warlock Apr 7 '13 at 8:29

Cleve Moler, in this note, gives an innocuous-looking DE he attributes to Larry Shampine that models flame propagation. The differential equation is

$$y^\prime=y^2-y^3$$

with initial condition $y(0)=\frac1{h}$, and integrated over the interval $[0,2h]$. The exact solution of this differential equation is

$$y(t)=\frac1{1+W((h-1)\exp(h-t-1))}$$

where $W(t)$ is the Lambert function, which is the inverse of the function $t\exp\,t$. (Whether the Lambert function is considered an analytical solution might well be up for debate, but I'm in the camp that considers it a closed form.)

For instance, in Mathematica, the following code throws a warning about possible stiffness:

With[{h = 40}, y /.
First @ NDSolve[{y'[t] == y[t]^2 - y[t]^3, y[0] == 1/h}, y,
{t, 0, 2 h}, Method -> {"ExplicitRungeKutta",
"DifferenceOrder" -> 4}]]

• But why is it a stiff problem, beside the fact that Matlab says so? – Ron Gordon Apr 6 '13 at 19:54
• @Ron: the working definition of "stiff" that I use is "...equations where certain implicit methods, in particular BDF, perform better, usually tremendously better, than explicit ones". Having also tried this out with my own implementations of explicit and implicit RK with stiffness detection, I agree with the assessment of Shampine and Moler. – J. M. is a poor mathematician Apr 6 '13 at 20:02
• I don't disagree with that assessment, it is just very nonobvious to me, especially given that I have a very dim grasp of how the Lambert W function behaves. Perhaps it is also that my numerical education came from, for better or worse, Numerical Recipes. There, the definition of "stiff" - and one that agrees with that given in the reference you provide - is when there are two dimension scales in a problem that are vastly different. I of course may be missing something, but this is why the stiffness here isn't clear to me. – Ron Gordon Apr 6 '13 at 20:08
• @Ron, maybe the asymptotics might make things clearer for you: $$W((h-1)\exp(h-t-1))=\exp(h-t-1)(h-1)-\exp(2 h-2 t-2) (h-1)^2+\frac32 \exp(3 h-3 t-3) (h-1)^3+\cdots$$ – J. M. is a poor mathematician Apr 7 '13 at 2:20
• Thank you for the pretty cool equation, but it is required an additional research. – Warlock Apr 8 '13 at 15:36