ルンゲクッタ四次法

Question

私は2つの方程式のシステムに対してルンゲクッタ四次法を実装しています。 T/Hがステップあるので

hは、セグメントの数です。

def cauchy(f1, f2, x10, x20, T, h): 
    x1 = [x10] 
    x2 = [x20] 

    for i in range(1, h): 
     k11 = f1((i-1)*T/h, x1[-1], x2[-1]) 
     k12 = f2((i-1)*T/h, x1[-1], x2[-1]) 
     k21 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12) 
     k22 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k11, x2[-1] + T/h/2*k12) 
     k31 = f1((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22) 
     k32 = f2((i-1)*T/h + T/h/2, x1[-1] + T/h/2*k21, x2[-1] + T/h/2*k22) 
     k41 = f1((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32) 
     k42 = f2((i-1)*T/h + T/h, x1[-1] + T/h*k31, x2[-1] + T/h*k32) 

     x1.append(x1[-1] + T/h/6*(k11 + 2*k21 + 2*k31 + k41)) 
     x2.append(x2[-1] + T/h/6*(k12 + 2*k22 + 2*k32 + k42)) 

    return x1, x2 

そこで私は、このシステム上でそれをテストしています：

def f1(t, x1, x2): 
    return x2 

def f2(t, x1, x2): 
    return -x1 

def true_x1(t): 
    return np.cos(t) + np.sin(t) 

def true_x2(t): 
    return np.cos(t) - np.sin(t) 

を私はまた異なる初期値と異なる機能でそれをテストした（正常に動作しているようだ。すべての作品ファイン）：

x10 = 1 
x20 = 1 
T = 1 
h = 10 

x1, x2 = cauchy(f1, f2, x10, x20, T, h) 
t = np.linspace(0, T, h) 

plt.xlabel('t') 
plt.ylabel('x1') 
plt.plot(t, true_x1(t), "blue", label="true_x1") 
plt.plot(t, x1, "red", label="approximation_x1") 
plt.legend(bbox_to_anchor=(0.97, 0.27)) 
plt.show() 

plt.xlabel('t') 
plt.ylabel('x2') 
plt.plot(t, true_x2(t), "blue", label="true_x2") 
plt.plot(t, x2, "red", label="approximation_x2") 
plt.legend(bbox_to_anchor=(0.97, 0.97)) 
plt.show() 

それから私は、エラーがO(step^4)のオーダーであるかどうかを確認したいので、私はこのようなステップと計算誤差を低減：

step = [] 
x1_error = [] 
x2_error = [] 
for segm in reversed(range(10, 1000)): 
    x1, x2 = cauchy(f1, f2, x10, x20, T, segm) 
    t = np.linspace(0, T, segm) 
    step.append(1/segm) 
    x1_error.append(np.linalg.norm(x1 - true_x1(t), np.inf)) 
    x2_error.append(np.linalg.norm(x2 - true_x2(t), np.inf)) 

そして、私はこの取得：

を

plt.plot(step, x1_error, label="x1_error") 
plt.plot(step, x2_error, label="x2_error") 
plt.legend() 

したがって、誤差はステップから線形です。 O(step^4)の順であるはずですから、これは本当に奇妙です。私が間違ったことを誰かに教えてもらえますか？

Answer 1

for i in range(1, h): 

これは1からh-1に反復します。最後のステップが欠落しているので、時間のx[h-1]と時刻Tの正確な解決の差はO(T/h)です。

したがってii+1からh手順について

for i in range(1,h+1): 

i-1 iから h手順について

又は

for i in range(h): 

を使用します。

---

また、np.linspace(0,1,4)は、最初は0で、最後はあなたが期待していたものを、おそらくされていない

array([ 0.  , 0.33333333, 0.66666667, 1.  ]) 

その結果、1ある4等間隔の数字を生成します。したがって、上記の補正では、両方の計算で同じ時点を使用することができます。

---

あなたがhまたはdtはステップサイズであり、Nはステップ数である彼らの通常の意味での文字を使用する場合は、あなたのコードに従うことが容易になるだろう。関数呼び出しでT/Nを繰り返し使用することを避けるには、ループh=T/Nまたはdt=T/Nの前に定義します。