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, -201*z - 200*z**2 + 2 ]) def explicit_euler(y0,t0,tEnd,N): y = zeros((2,N+1)) t = zeros(N+1) y[:,0] = y0 t = 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 = 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 yi, ti= implicit_euler(y0,t0,tEnd,N) print yi
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?