Checking if a given differential equation is stiff

Question

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?

Answer 1

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.

Answer 2

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.