[Pythonによる科学・技術計算]4次のルンゲ-クッタ法による1次元ニュートン方程式の解法

はじめに

常微分方程式の解法の一つである4次のルンゲ-クッタ法によるニュートン方程式の数値解法の例を挙げる。

---

内容

(1) [ウォーミングアップ] 4次のルンゲ-クッタ法による1階常微分方程式の解法。
(2) 調和振動子
(3) 減衰振動
(4) 強制振動

---

内容(1): 1階常微分方程式の解法

import numpy as np
import matplotlib.pyplot as plt
"""
4次のルンゲ-クッタ法による1階常微分方程式
The fourth-order Runge-Kutta method for the 1st order ordinary differencial equation

dx/dt = -x^3+sin(t)を解く
"""

def f(x,t):
    return -x**3 +np.sin(t)

a = 0.0
b = 10.0
N = 100
h = (b-a)/N

tpoints = np.arange(a,b,h)
xpoints = []
x = 0.0

for t in tpoints:
    xpoints.append(x)
    k1 = h*f(x,t)
    k2 = h*f(x+0.5*k1, t+0.5*h)
    k3 = h*f (x+0.5*k2, t+0.5*h)
    k4 = h*f(x+k3, t+h)
    x += (k1+2*k2+2*k3+k4)/6

plt.plot (tpoints, xpoints)
plt.xlabel("t")
plt.ylabel("x(t)")
plt.show()

---

内容(2):調和振動子

import numpy as np
import matplotlib.pyplot as plt
"""
4次のルンゲ-クッタ法による1次元ニュートン方程式の解法
The fourth-order Runge-Kutta method for the 1D Newton's equation
# 例: 線形復元力 (harmonic oscilation)
"""

def f(x, v, t):
    return -k*x  # 復元力

M = 1.0 # mass: 1 Kg
k = 1.0 #ばね定数
t0 = 0.0
t1 = 10.0

N = 50



del_t = (t1-t0)/N # time grid

tpoints = np.arange(t0, t1, del_t)
xpoints = []
vpoints = []

# initial condition
x0 = 0.0
v0 = 1.0

x, v = x0, v0


for t in tpoints:
    xpoints.append(x)
    vpoints.append(v)
    k1v =f(x, v, t)*del_t /M
    k1x = v * del_t

    k2v =  f(x+k1x/2, v+k1v/2, t+del_t/2 )*del_t /M
    k2x =(v+k1v/2)*del_t 

    k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t /M
    k3x =(v+k2v/2 ) *del_t 

    k4v = f(x+k3x, v+k3v, t+del_t )*del_t /M
    k4x = (v+k3v )*del_t 



    v += (k1v + 2 * k2v + 2* k3v + k4v)/6
    x += (k1x + 2 * k2x + 2* k3x + k4x)/6

plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t",  fontsize=24)
plt.ylabel("x(t)",  fontsize=24)

tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper left')

plt.show()


plt.plot (tpoints,vpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)

tt=np.arange(t0, t1, 0.01*del_t)
plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper left')

plt.show()

---

内容(3)減衰振動

import numpy as np
import matplotlib.pyplot as plt
"""
減衰振動
Damped oscillation

"""

def f(x, v, t):
    k=1.0 
    a=1.0
    return -k*x-a*v

M = 1.0 # mass: 1 Kg
t0 = 0.0
t1 = 10.0

N = 50



del_t = (t1-t0)/N # time grid

tpoints = np.arange(t0, t1, del_t)
xpoints = []
vpoints = []

# initial condition
x0 = 0.0
v0 = 1.0

x, v = x0, v0


for t in tpoints:
    xpoints.append(x)
    vpoints.append(v)
    k1v =f(x, v, t)*del_t /M
    k1x = v * del_t

    k2v =  f(x+k1/2, v+k1v/2, t+del_t/2 )*del_t /M
    k2x =(v+k1v/2)*del_t 

    k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t /M
    k3x =(v+k2v/2 ) *del_t 

    k4v = f(x+k3, v+k3v, t+del_t )*del_t /M
    k4x = (v+k3v )*del_t 



    v += (k1v + 2 * k2v + 2* k3v + k4v)/6
    x += (k1x + 2 * k2x + 2* k3x + k4x)/6

plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t",  fontsize=24)
plt.ylabel("x(t)",  fontsize=24)

#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper right')

plt.show()


plt.plot (tpoints,vpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)

#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper right')

plt.show()

---

内容(4)強制振動

import numpy as np
import matplotlib.pyplot as plt
"""
強制振動
Forced vibration

"""

def f(x, v, t):
    k=1.0 
    gamma=0.1
    return -k*x-2* gamma *v + 2*sin(t)

M = 1.0 # mass: 1 Kg
t0 = 0.0
t1 = 50.0

N = 1000



del_t = (t1-t0)/N # time grid

tpoints = np.arange(t0, t1, del_t)
xpoints = []
vpoints = []

# initial condition
x0 = 0.0
v0 = 1.0

x, v = x0, v0


for t in tpoints:
    xpoints.append(x)
    vpoints.append(v)
    k1v =f(x, v, t)*del_t /M
    k1x = v * del_t

    k2v =  f(x+k1/2, v+k1v/2, t+del_t/2 )*del_t /M
    k2x =(v+k1v/2)*del_t 

    k3v =f (x+k2x/2, v+k2v/2, t+del_t/2 )*del_t /M
    k3x =(v+k2v/2 ) *del_t 

    k4v = f(x+k3, v+k3v, t+del_t )*del_t /M
    k4x = (v+k3v )*del_t 



    v += (k1v + 2 * k2v + 2* k3v + k4v)/6
    x += (k1x + 2 * k2x + 2* k3x + k4x)/6

plt.plot (tpoints, xpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t",  fontsize=24)
plt.ylabel("x(t)",  fontsize=24)

#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.sin(tt), '-',label='Exact: Sin (t)',color = 'Green')
plt.legend(loc='upper right')

plt.show()


plt.plot (tpoints,vpoints, 'o',label='4th order Runge-Kutta')
plt.xlabel("t", fontsize=24)
plt.ylabel("v(t)", fontsize=24)

#tt=np.arange(t0, t1, 0.01*del_t)
#plt.plot (tt, np.cos(tt), '-',label='Exact: Cos (t)', color = 'red')
plt.legend(loc='upper right')

plt.show()

---

補遺

ここで行ったルンゲ-クッタ法では，系の保存量の一つである全エネルギーは時間ステップの増加に伴い数値的には保存されなくなる。これを改善する(シンプレクティック性を持たせる)解法の一つがベルレー法である(Qiita記事を参照)。