ルンゲクッタ法 誤差 python

Question

ルンゲクッタ法を用いて微分方程式を解きたいのですが、ループ文がうまくつくれません。
具体的には関数RKで出た値'u0'をfor文で複数回繰り返す際に、値を更新していって欲しいのですが、コードを書くことが出来ずに同じ計算結果がリストに溜まってしまいます。

どのようなコードを書けばよいでしょうか？
初歩的なテクニックかもしれませんが、うまく書けません。
よろしくお願いいたします。

#各モジュールをインポート
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

#時間経過を定義
dt_1=0.0001

#時間の初期値を入力
t0=0

#時間の最大値
t_max=10
#試行回数
trials_1=np.arange(t0,t_max,dt_1)

#関数の値を入れるオブジェクトを生成
a_points_1=[]

#関数を定義
def u(t,x,x1):
    return -4*x
#x'の式
def y(t,x,x1):
    return x1

#runge-kuttaをtrials回実行
#結果をa_pointsに入れる
def rp():
    t,y0,u0=0,1,0
    for t in trials_1:
        x1=RK(y0,u0,dt_1)
        a_points_1.append(x1)
    print(a_points_1[:3])
        
#runge-kuttaの実行式
#k:y, l:u
def RK(y0,u0,dt):
    
    k1=y0*dt
    l1=u0*dt
    
    k2=(u0+l1*0.5)*dt
    l2=-4*(y0+k1*0.5)*dt
    
    k3=(u0+l2*0.5)*dt
    l3=-4*(y0+k2*0.5)*dt
    
    k4=(u0+l3)*dt
    l4=-4*(y0+k3)*dt
    
    k=(k1+2*k2+2*k3+k4)/6
    l=(l1+2*l2+2*l3+l4)/6
    
    y0=y0+k
    u0=u0+l
    
    return u0
#このu0の値をfor文で繰り返す際に更新していきたい
#現状では初期値のu0=0が常に代入されて、同じ計算結果が返ってくる
    
#プログラムを実行
rp()

Answer 1

デバッグに伴って、逆にわかりにくくなっているので、落ち着いて書き直したほうが早いかもしれません。
実装自体が目的と想定して、正しい値を出すであろうscipyの結果との比較を入れました。

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
%matplotlib inline

# initial values
dt = 0.001
t0 = 0.
t1 = 10.
x0 = np.array([0., 1.])

# time
N = int((t1-t0)/dt)
dt = (t1-t0)/N
trials = np.linspace(t0, t1, N)

# target function
def f(t, x):
    A = np.array([[0., 1.], [-1., 0.]])
    return np.dot(A, x)

# integrator
def gen_evolve(dt, f):
    def evolve(x, t):
        k1 = f(t, x)
        k2 = f(t + dt/2., x + k1*dt/2.)
        k3 = f(t + dt/2., x + k2*dt/2.)
        k4 = f(t + dt, x + k3*dt)
        k = (k1 + 2.*k2 + 2.*k3 + k4)*dt/6.
        x = x + k
        t = t + dt
        return x, t
    return evolve

# runge-kutta
def rp(evolve, x0, t0, t1):
    ans = []
    x = x0.copy()
    t = t0
    while t < t1:
        x, t = evolve(x, t)
        ans.append(x)
    return np.array(ans)

# self-defined integrator
ans0 = rp(gen_evolve(dt, f), x0, t0, t1)

# from scipy
r = scipy.integrate.ode(f)
r.set_integrator('dopri5')
r.set_initial_value(x0, t0)
ans1 = []
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    ans1.append(r.y)
ans1 = np.array(ans1)

# two methods should give same result
assert(np.linalg.norm(ans0-ans1)<0.1)

# plot
fig, ax = plt.subplots(ncols=3, dpi=200)
ax[0].set_aspect('equal')
ax[0].scatter(ans0[::100, 0], ans0[::100, 1], s=0.5, color='blue')
ax[1].set_aspect('equal')
ax[1].scatter(ans1[::100, 0], ans1[::100, 1], s=0.5, color='red')
ax[2].set_aspect('equal')
ax[2].scatter(ans0[::100, 0], ans0[::100, 1], s=0.5, color='blue')
ax[2].scatter(ans1[::100, 0], ans1[::100, 1], s=0.5, color='red')
plt.show()

---

Runge-Kuttaの式に関しては、

k1=y0*dt
l1=u0*dt

が多分正しくない気がします。

Comments

[sshhoo]

ご回答誠にありがとうございます。
ご指摘頂いた通り、自作したプログラムだと解析解との乖離が依然としてありましたので、大変助かりました。
お伺いしたいのですが、

# target function
def f(t, x):
A = np.array([[0., 1.], [-1., 0.]])
return np.dot(A, x)

Aに関して、なぜこのような行列が出てきたのかわからないのですが、解説して頂けないでしょうか？

お手数ですがよろしくお願い致します。

[mkgrei]

もとの問題は、

dx/dt = y
dy/dt = -4x

だったので、
v = [[x], [y]]として、

dv/dt = Av
A = [[0,1],
[-4,0]]

となります。
-4ではなく-1になっているのは手元で丸い円をとりあえず描こうとして変更し忘れたからです。
ぼーっとしていました。

[sshhoo]

ご解説ありがとうございます。
行列で少し躓きましたが、コードを参考にしながらようやく書き終えることが出来ました。
御礼申し上げます。

Answer 2

正直ルンゲクッタ法の事は名前しか覚えていません。
ですから、想像の範囲になりますが、こういう風に書けばいいんじゃないかなぁ。

def rp():
    y0, u0 = 1, 0
    for _ in trials_1:
        y0, u0 = RK(y0, u0)
        a_points_1.append(u0)

    print(a_points_1[:3])

def RK(y0, u0, dt=0.0001):
    ...

    return y0, u0

関数内での再代入は、関数外には一切影響しません。(グローバル宣言している場合を除く)
基本的に数値や文字列を関数に与えた場合、その値が変わることはあり得ません。

Answer 3

Pythonのintやfloatはimmutableなので、関数の中で値を更新されてもその外には反映されません。たとえば、numpy.arrayに変換してあげると、関数内で数値を更新させることができます。

t = np.array(0.0, dtype=np.float)
y0 = np.array(1.0, dtype=np.float)
u0 = np.array(0.0, dtype=np.float)

これをt,y0,u0=0,1,0の代わりに使ってください。
それから、値を更新する部分：

    y0=y0+k
    u0=u0+l

    return u0

を、次のようにしてください(LouiS0616さんのコメントを受けて修正)。

    y0 += k
    u0 += l

    return u0.copy()

以上を踏まえてコード全体を再掲

#各モジュールをインポート
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

#時間経過を定義
dt_1=0.0001

#時間の初期値を入力
t0=0

#時間の最大値
t_max=10
#試行回数
trials_1=np.arange(t0,t_max,dt_1)

#関数の値を入れるオブジェクトを生成
a_points_1=[]

#関数を定義
def u(t,x,x1):
    return -4*x
#x'の式
def y(t,x,x1):
    return x1

#runge-kuttaをtrials回実行
#結果をa_pointsに入れる
def rp():
    y0 = np.array(1.0, dtype=np.float)
    u0 = np.array(0.0, dtype=np.float)
    for t in trials_1:
        x1=RK(y0,u0,dt_1)
        a_points_1.append(x1)
    print(a_points_1[:3])

#runge-kuttaの実行式
#k:y, l:u
def RK(y0,u0,dt):

    k1=y0*dt
    l1=u0*dt

    k2=(u0+l1*0.5)*dt
    l2=-4*(y0+k1*0.5)*dt

    k3=(u0+l2*0.5)*dt
    l3=-4*(y0+k2*0.5)*dt

    k4=(u0+l3)*dt
    l4=-4*(y0+k3)*dt

    k=(k1+2*k2+2*k3+k4)/6
    l=(l1+2*l2+2*l3+l4)/6

    y0+=k
    u0+=l

    return u0.copy()  # もしくはfloat(u0)
#このu0の値をfor文で繰り返す際に更新していきたい
#現状では初期値のu0=0が常に代入されて、同じ計算結果が返ってくる

#プログラムを実行
rp()

Comments

[LouiS0616]

イミュータブルであることは関係しないのでは？再代入ですから。

---
def hoge(lst):
...lst = [1, 2, 3]

lst = []
hoge(lst)
print(lst) # []

[退会済みユーザー]

あー関数の中y0+=kとやる前提での対応法をかいておりました。そこも修正がいりますね。

[退会済みユーザー]

修正しました。これでどうでしょうか。

[退会済みユーザー]

イミュータブルに関するご指摘ですが、確かにa = a+1は再代入が発生してしいます。ただし、（objectのメソッドの定義次第ですが）これはa += 1と同じではなく、inplaceで行われるものがあります。今回の回答はnumpy.arrayがmutableであること、+=がinplaceで行われることを利用した回答です。そもそも関数の中でオブジェクトの状態を変えてしまうので、トリッキーと言えばトリッキーですし、デバッグしにくくなります。ただ、質問者さんは関数の中で値を更新したい、とのことでしたので、それを実現するための方法を提案してみました。

[LouiS0616]

迅速な対応感謝します。

うーん...
個人的には0次元配列、あるいは要素数1のコンテナは避けたい感があります。
ただこれに関しては単に好みの問題かもしれません。

[退会済みユーザー]

その点に関しては同意します。やるとしてもlistやdictに放り込んでしまうほうがマシな気がします。

[LouiS0616]

確かに文字通り『関数の中で値を更新したい』のならこのようになるでしょうね。

[sshhoo]

お二方ともご回答頂きまして誠にありがとうございます。

少し頭を整理してもう一度再考してみました。
やりたいこととしましては、関数RK内で初期値y0,u0を与えて関数内の計算を処理したのちに
y0=y0+k
u0=u0+l
として、値を更新した後に、更新した値を関数RKの初期値として用い同様の計算を行う行程を複数回実行して、微分方程式の数値解を求めたいという流れです。

関数RKから帰ってくる値をappendでリストに入れていますが、リストに入ってくる値は全て同じ値になってしまいます。

つまり、関数RKが永遠と初期値を使って同じ計算をしているのかなと考えております。

y0=y0+k
u0=u0+l
を
y_new=y0+k
u_new=u0+l
と書き直しますと

関数RKを実行するたびに、毎回初期値であるy0,u0を用いて計算するのではなく、最新のy_new,u_newの値を使って関数RKを動かしたいと考えております。

具体的には
for i in trials_1:
で(0<n<=i)
n回目の関数RKの施行時だと、RKの引数にはn-1回目で得られたy_new,u_newの値を用いて計算を行うプログラムを作りたいです。

質問力が貧弱であるため、うまく伝えられなかったと思い長々とコメントさせて頂きました。

質問に関してまた何かいい方法があれば、お返事をよろしくお願い申し上げます。

正直ルンゲクッタ法を誤って解釈している可能性も捨てきれません。
その際はご迷惑をおかけしたことを予め謝らせて下さい。

[退会済みユーザー]

念のためコード全体を再掲しました。今回の問題になっているのは、関数内での値の更新が、Pythonの組み込みのデータ(intやfloat)の性質上やりにくい点にあります。私のコードはnumpyライブラリに頼った方法ですが、LouiS0616のほうも見比べて使いやすいほうを選んでいただければいいんじゃないかと思います。なお、ルンゲタックの式があっているかは確認できていません。悪しからず・・・

[sshhoo]

大変丁寧にご回答頂きまして誠にありがとうございます。
自分なりに理解できたつもりでいます。
助かりました。
コメントで
a=a+1 と
a += 1
は同一ではないとおっしゃっていましたが、どのような違いがあるのですか？
初心者が理解できる範囲外である場合は、そう答えて頂いて結構です。

[退会済みユーザー]

最初に呼び出そうとするメソッドが違います。a+1はa.__add__を最初に、a+=1はa.__iadd__をそれぞれ最初に探しに行きます。ただし、aに__iadd__が無ければ、代わりにa=a.__add__(1)と再代入が発生する形で評価します。

[退会済みユーザー]

ためしに
(immutableであるintクラス)
int.__dict__を覗いてみると、__add__, __radd__はありますが、__iadd__はありません。

(mutableなnumpy.ndarrayクラス)
一方、numpy.ndarray.__dict__には__iadd__もあります。ちなみに、iaddのiはinplaceのiだと思います。

なお、immutableとはオブジェクトの中身の変更を許さない性質を意味するので、inplaceで作用するiaddがないのも当然ですね。

[sshhoo]

詳細なご回答誠にありがとうございます。
勉強になりました。

Answer 4

pythonは使った事がないけど確実に言える事があります。

return u0 していて x1=RK(y0,u0,dt_1) となっているなら
戻り値はx1に入っている訳でして、戻り値を u0 に反映させたいなら
以下のようにする必要があると思われます。

x1=RK(y0,u0,dt_1)
u0=x1 ←この行が足りない

それと引数に対して値を代入するのは基本的にまずいので
y0とu0をグローバル変数にして引数で渡さないようにするか、
以下のようにして別の変数に計算結果を入れて値を返す様にするべきです。

まずい：
y0=y0+k
u0=u0+l
return u0

よい：
y1=y0+k
u1=u0+l
return u1

お二方の回答を見て、そもそもの前提を見落としていた事に気が付きました、お恥ずかしい限りです。
この回答は無視して下さいますようお願いいたします。

Comments

[sshhoo]

ご協力並びにお心遣いありがとうございました。