質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Q&A

解決済

1回答

430閲覧

python ルンゲクッタ法

sshhoo

総合スコア15

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

0グッド

0クリップ

投稿2018/02/06 06:29

ルンゲクッタ法のコーディングに関して、よくわからないことがあったためご質問させて頂きました。

<やりたいこと>
四次のルンゲクッタ法を用いて一階の常微分方程式の数値解を求めたい。
また、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)")

イメージ説明

よろしくお願い致します。

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答1

0

ベストアンサー

2つ目のfor文でxの値を確認してみてください。rk()に意図した値を渡していますか?
また、ループ中x,x_2は意図した値になっていますか?

# このあたりで何かが必要? for t in t_points_2: print(x) # 意図した値? x_2=rk(x,t,f,dt_2) # このあたりで何かが必要? x_points_2.append(x_2)

投稿2018/02/06 06:36

編集2018/02/06 06:47
can110

総合スコア38234

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

sshhoo

2018/02/06 06:50

確認したところ 3.14159265297 が t の値に関わらず一様に与えられていました。 しかし、このプログラムでxを上記の値に定義したつもりはなく、関数run()内で、 x=0 と、初期値を定義したのにどうしてrk()に渡されないのかが疑問です。 何かややこしいコードを書いてしまったのでしょうか?
can110

2018/02/06 07:01

1個目のループではx=0から始まりxの値を更新していますね。 で、終了時点でxの値は3.14~になったわけですよね? で、2個目のループではrk(x,~)ではxの値を渡しているだけなのでxの値はそのまま変わらず。 それはおかしいのでは?ということです。
sshhoo

2018/02/06 07:31

ご指摘ありがとうございます。 xの値が一回目のfor文の結果を引きずってしまっていたのですね。 理解出来ました。 x_2=0 for t in t_points_2: if t==0: x=x_2 x2=rk(x,t,f,dt_2) x_points_2.append(x2) else: x=x2 x2=rk(x,t,f,dt_2) x_points_2.append(x2) if文をつけてみると、思い描いていた結果となりました。 お力添えして頂き、誠にありがとうございました。
can110

2018/02/06 07:42

う~ん、詳細未検証ですが、単純に x = 0 for ~   x = rk(x,~   x_points_2.append(x) でよいかと思います。
sshhoo

2018/02/06 07:53

試してみたところ、ご提案頂いたやり方の方がより短いコードで同じ結果が得られました。 アドバイスありがとうございます。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問