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

Question

我试图写一个算法来解决非线性 ODE

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 的示例，但对于我的公式来说，它不稳定并且会爆炸（对于RK4 和欧拉，t > 0.65）。因此是我的实现不正确，还是我没有看到另一个错误？

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）与精确解进行比较。