Checking if a given differential equation is stiff

i have to decide if the following differential equation is stiff: $$y''(t)=-201y'-200y^2 + 2, \quad t\in[0,20].$$

Sadly, I don't have any solutions. So, what I did was implementing explicit and implicit euler and look at the result for various step sizes.

from numpy import *
from scipy import optimize

rhs = lambda z: array([ z[1], -201*z[1] - 200*z[0]**2 + 2 ])

def explicit_euler(y0,t0,tEnd,N):
y = zeros((2,N+1))
t = zeros(N+1)

y[:,0] = y0
t[0] = t0

h = float(tEnd-t0)/N

for k in xrange(N):
y[:,k+1] = y[:,k] + h*rhs(y[:,k])
t[k+1] = t[k] + h

return y, t

def implicit_euler(y0,t0,tEnd,N):
y = zeros((2,N+1))
t = zeros(N+1)

y[:,0] = y0
t[0] = t0

h = float(tEnd-t0)/N

for k in xrange(N):
F = lambda a: -a + y[:,k] + h*rhs(0.5*(y[:,k]+a))
startvalue = y[:,k] + h*rhs(y[:,k])
y[:,k+1] = optimize.fsolve(F, startvalue)
t[k+1] = t[k] + h

return y, t

t0 = 0.
tEnd = 20.
y0 = array([1,0])
N = 100

ye, te = explicit_euler(y0,t0,tEnd,N)
print ye[0]
yi, ti= implicit_euler(y0,t0,tEnd,N)
print yi[0]

For $N=100$ steps, like in the code above, the result output

is:

[  1.00000000e+000   1.00000000e+000  -6.92000000e+000   2.95624000e+002
-1.19471120e+004  -2.31180176e+005  -1.13350513e+009  -3.84263356e+011
...
nan               nan               nan               nan
nan               nan               nan               nan
nan]
[ 1.          0.84122674  0.7135129   0.63248909  0.55673461  0.50831417
0.45792912  0.42604819  0.39009643  0.36765118  0.34076357  0.32415616
....
0.10461852  0.10443526  0.10425726  0.1040885   0.10392479  0.10376937
0.10361875  0.10347557  0.10333695  0.10320504  0.10307743]

So we see, that the explicit euler diverges. (all the nans) but the implicit doesn't. By setting $N=100000$ we also see that both get the same result. [We expect a stiff differential equation to be solvable using an explicit method when the stepsize is small enough]

But that seems a little to easy here. Would you say that's fine? How would you do it?

• I agree with the two answers that "stiffness" of an ODE is a very subtle concept. I also agree that Shampine has a lot of insight on the subject. So does Cleve Moler. I strongly suggest you take a look at this chapter of Moler's book, mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/…. Contrary to what was stated in one of the answers, he shows an example (page 19), due to Shampine, of a single ODE that exhibits stiffness. – Bill Greene Feb 1 '17 at 22:46
• Yeah, I know about the stability function (hope that's the proper term). Thanks for the clarification :) – xotix Feb 2 '17 at 10:20

There is no solid definition for stiff equations. I like Shampine's working definition the best: a differential equation is stiff if explicit methods are less computationally efficient than implicit methods. All of the other definitions about eigenvalues are purposefully vague because they cannot capture the fact that the "eigenvalues which matter" is dependent on how fast the methods are, and so in some cases what may seem "stiff" by eigenvalues is something you don't actually want to treat as "stiff" because explicit methods may step well enough that they work better.

Notice that then has a lot to do with the program you are using. Since you are using looping with NumPy, this is going to a be a very inefficient implementation. Because your implementation is inefficient, the fact that the kernal in fsolve is efficient may mean that, for this comparison, the implicit method will have less of a disadvantage due to implementation, and thus your threshold for stiffness would be lower than if you were using a heavily optimized code.

I encourage you to adopt this working definition and use the idea of stiffness and the eigenvalues estimates as a guideline for choosing methods, but double-check intuition with benchmarks.

• yeah, that's basically what's written in my lecture notes. We call a differential equation stiff, if an explicit method needs very small steps. It is solvable by a implicit with "normal" steps though. This is a basic introduction in the first year of my studies, so we don't really care too much about the most efficient implementations. Thanks for the clarification but my question is still open. :) – xotix Feb 2 '17 at 10:18
• With this definition, the way to check if an equation is stiff is simple: run both an explicit and an implicit method and see which one is faster. Multiple timescales serves as a guideline, but ultimately benchmarks are the only practical way to see. – Chris Rackauckas Feb 2 '17 at 13:00

The term "stiff" has many definitions, but a commonly understood one would be "there are two or more timescales associated with the process being modeled, and these timescales are very different".

In your case, you have a single ODE whose solution is a single function $y(t)$. As a consequence, it has only a single timescale, and so one would not consider the ODE stiff.

• I'm doing a basic lecture in numerical methods. So our definition probably differs. Our definition is literally "we call a differential equation stiff, if an explicit method needs very small steps to solve it." Which results in time consuming claculation although it might be very easy to calculate at some points. But thanks for the clarification, seems like we use a very basic definition then. – xotix Feb 2 '17 at 10:15
• @xotix: But what does "very small" actually mean? Small compared to what? Your comment seems to suggest that it is small compared to the length of time you want to integrate the ODE over, but that makes it relative to something you choose, not to the properties of the equation itself. – Wolfgang Bangerth Feb 3 '17 at 20:47