计算科学 - 检查给定的微分方程是否是刚性的 - 吾爱随笔录

我必须确定以下微分方程是否僵硬：

y^{″} (t) = - 201 y^{'} - 200 y^{2} + 2, t \in [0, 20] .

$y''(t)=-201y'-200y^2 + 2, \quad t\in[0,20].$

可悲的是，我没有任何解决方案。所以，我所做的是实现显式和隐式欧拉，并查看各种步长的结果。

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]

对于 $N=100$ 步，如上面的代码，结果输出

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]

所以我们看到，显式欧拉发散了。（所有的nans）但隐含的没有。通过设置，我们还看到两者都得到相同的结果。[当步长足够小时，我们期望刚性微分方程可以使用显式方法求解] $N=100000$

但这似乎有点容易。你说这样好吗？你会怎么做？