ルンゲクッタ法のコーディングに関して、よくわからないことがあったためご質問させて頂きました。
<やりたいこと>
四次のルンゲクッタ法を用いて一階の常微分方程式の数値解を求めたい。
また、dtを変化させたときの数値解の値をグラフにして表示し、比較したい。
<質問内容>
関数’run’の中にある二つのfor文では同じことをさせているつもりなのですが、全く異なる値が返ってきます。
この原因を教えて頂きたいです。
python
1#各モジュールをインポート 2import numpy as np 3import matplotlib.pyplot as plt 4%matplotlib inline 5 6#今回の関数を定義 7def f(x,t): 8 return 4/(1+t**2) 9 10#t軸の値とdtを定義 11t_min = 0.0 12t_max = 1 13 14dt_1 = 0.1 15t_points_1 = np.arange(t_min,t_max,dt_1) 16 17dt_2=0.01 18t_points_2=np.arange(t_min,t_max,dt_2) 19 20#x軸の値を入れるオブジェクトを作成 21x_points_1 = [] 22x_points_2=[] 23 24#ルンゲクッタの関数 25def rk(x,t,f,dt): 26 k1 = dt*f(x,t) 27 k2 = dt*f(x+k1/2,t+dt/2) 28 k3 = dt*f(x+k2/2,t+dt/2) 29 k4 = dt*f(x+k3,t+dt) 30 x = x+(k1+2*k2+2*k3+k4)/6 31 return x 32 33#この関数内のfor文では同じことを実行しているつもりなのですが、後者から帰ってくる値が予想と大きく異なり、原因がわかりません。 34def run(): 35 x=0 36 for t in t_points_1: 37 k1 = dt_1*f(x,t) 38 k2 = dt_1*f(x+0.5*k1, t+0.5*dt_1) 39 k3 = dt_1*f (x+0.5*k2, t+0.5*dt_1) 40 k4 = dt_1*f(x+k3, t+dt_1) 41 x += (k1+2*k2+2*k3+k4)/6 42 x_points_1.append(x) 43 44 for t in t_points_2: 45 x_2=rk(x,t,f,dt_2) 46 x_points_2.append(x_2) 47 48plt.plot(t_points_1,x_points_1, label='dx/dt=4/(1+t^2)_1') 49plt.legend() 50plt.xlabel("t") 51plt.ylabel("x(t)") 52 53plt.plot(t_points_2,x_points_2, label='dx/dt=4/(1+t^2)_2') 54plt.legend() 55plt.xlabel("t") 56plt.ylabel("x(t)")
よろしくお願い致します。
回答1件
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2018/02/06 06:50
2018/02/06 07:01
2018/02/06 07:31
2018/02/06 07:42
2018/02/06 07:53