RK4方法和欧拉方法只适用于某些公式

Question

我试着写一个算法来解非线性常微分方程

dr/dt = r(t)+r²(t)  

哪一个有（一种可能的）解决办法

r(t) = r²(t)/2+r³(t)/3  

为此，我在python中实现了euler方法和RK4方法。对于错误检查，我使用了rosettacode上的示例：

dT/dt = -k(T(t)-T0)

解决方案

T0 + (Ts − T0)*exp(−kt)

因此，我的代码现在看起来像

import numpy as np
from matplotlib import pyplot as plt

def test_func_1(t, x):
    return x*x

def test_func_1_sol(t, x):
    return x*x*x/3.0

def test_func_2_sol(TR, T0, k, t):
    return TR + (T0-TR)*np.exp(-0.07*t)

def rk4(func, dh, x0, t0):
    k1 = dh*func(t0, x0)
    k2 = dh*func(t0+dh*0.5, x0+0.5*k1)
    k3 = dh*func(t0+dh*0.5, x0+0.5*k2)
    k4 = dh*func(t0+dh, x0+k3)
    return x0+k1/6.0+k2/3.0+k3/3.0+k4/6.0

def euler(func, x0, t0, dh):
    return x0 + dh*func(t0, x0)

def rho_test(t0, rho0):
    return rho0 + rho0*rho0

def rho_sol(t0, rho0):
    return rho0*rho0*rho0/3.0+rho0*rho0/2.0

def euler2(f,y0,a,b,h):
    t,y = a,y0
    while t <= b:
        #print "%6.3f %6.5f" % (t,y)
        t += h
        y += h * f(t,y)

def newtoncooling(time, temp):
    return -0.07 * (temp - 20)

x0 = 100
x_vec_rk = []
x_vec_euler = []
x_vec_rk.append(x0)

h = 1e-3
for i in range(100000):
    x0 = rk4(newtoncooling, h, x0, i*h)
    x_vec_rk.append(x0)

x0 = 100
x_vec_euler.append(x0)
x_vec_sol = []
x_vec_sol.append(x0)

for i in range(100000):
    x0 = euler(newtoncooling, x0, 0, h)
    #print(i, x0)
    x_vec_euler.append(x0)
    x_vec_sol.append(test_func_2_sol(20, 100, 0, i*h))

euler2(newtoncooling, 0, 0, 1, 1e-4)

x_vec = np.linspace(0, 1, len(x_vec_euler))

plt.plot(x_vec, x_vec_euler, x_vec, x_vec_sol, x_vec, x_vec_rk)
plt.show()

#rho-function
x0 = 1
x_vec_rk = []
x_vec_euler = []
x_vec_rk.append(x0)

h = 1e-3
num_steps = 650
for i in range(num_steps):
    x0 = rk4(rho_test, h, x0, i*h)
    print "%6.3f %6.5f" % (i*h, x0)
    x_vec_rk.append(x0)

x0 = 1
x_vec_euler.append(x0)
x_vec_sol = []
x_vec_sol.append(x0)

for i in range(num_steps):
    x0 = euler(rho_test, x0, 0, h)
    print "%6.3f %6.5f" % (i*h, x0)
    x_vec_euler.append(x0)
    x_vec_sol.append(rho_sol(i*h, i*h+x0))

x_vec = np.linspace(0, num_steps*h, len(x_vec_euler))
plt.plot(x_vec, x_vec_euler, x_vec, x_vec_sol, x_vec, x_vec_rk)
plt.show()

对于rosettacode中的示例，它工作得很好，但是对于我的公式，它是不稳定的，并且会爆炸（对于t>；0.65，对于RK4和euler）。因此，我的实现是不正确的，还是我没有看到另一个错误？你知道吗

Answer 1

搜索方程的精确解：

dr/dt = r(t)+r²(t)

我发现：

r(t) = exp(C+t)/(1-exp(C+t))

其中C是依赖于初始条件的任意常数。可以看出，对于t -> -C，r(t) -> infinity。你知道吗

我不知道你用什么初始条件，但可能是你在计算数值解时遇到这个奇点。你知道吗

更新： 因为初始条件是r(0)=1，所以常数C是C = ln(1/2) ~ -0.693。它可以解释为什么你的数值解在某些t>；0.65时崩溃

更新： 为了验证你的代码，你可以简单的比较你的数值解，比如说0<=t<0.6和精确解。你知道吗