RK4-метод и метод Эйлера работает только для определенных формул

Question

Я попытался написать алгоритм для решения нелинейных ОДУ

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

, который имеет (один из возможных) раствор

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

Для этого Я реализовал как метод эйлера, так и метод RK4 в python. Для проверки ошибок я использовал пример на 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, так и для эйлера) для моей формулы. Это значит, что моя реализация неверна, или есть еще одна ошибка, которую я не вижу?

Answer 1

В поисках точного решения вашего уравнения:

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

Я нашел:

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

Где C произвольная постоянная, зависящая от начальных условий. Можно видеть, что для t -> -Cr(t) -> infinity.

Я не знаю, какое начальное условие вы используете, но может быть, вы встречаетесь с этой особенностью при вычислении численного решения.

ОБНОВЛЕНИЕ: Поскольку ваше начальное состояние - r(0)=1, константа C - C = ln(1/2) ~ -0.693. Это может объяснить, почему ваше числовое решение падает с некоторого t> 0,65

ОБНОВЛЕНИЕ: Для проверки кода вы можете просто сравнить свое числовое решение, рассчитанное для, скажем 0<=t<0.6 с точным решением.