RK4法および特定の式

Question

のみ作動オイラー法は、私はそのための解決策

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

有する非線形ODE

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

を（一つの可能​​）を解くためのアルゴリズムを記述しようとしました私はPythonでオイラー法とRK4法の両方を実装しました。 、それはrosettacodeから例えば正常に動作

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() 

よう

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

はこのように、私のコードは今見え溶液で

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

：エラーチェックのために私はrosettacode上の例を使用しました私の処方では不安定で爆発します（t> 0.65、RK4とオイラーの両方）。したがって、私の実装は間違っているのでしょうか、それとも私には見えない別のエラーがありますか？

Answer 1

あなたの方程式の厳密解を検索：

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

私が見つけた：Cは初期条件に依存して、任意の定数である

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

。 t -> -Cについてはr(t) -> infinityであることが分かる。

私はどの初期条件を使用するのか分かりませんが、数値解を計算する際にこの特異点を満たしている可能性があります。

更新日： 初期条件はr(0)=1なので、定数CはC = ln(1/2) ~ -0.693です。数値解がいくつかの理由でクラッシュする理由を説明できます。0.30

更新日： コードを検証するには、たとえば0<=t<0.6で計算された数値解を正確な解と単純に比較できます。