ルンゲクッタとか懐かしい

参考:
常微分方程式における数値積分の基礎とPythonサンプルコード-MyEnigma

なんか大学1,2年の時にロボットアームの制御とか,目でみて誤差を体感しよう的な事をやったなあって.懐かしい.

参考記事をなぞってみる(odeintと統一的に書きたかったのでコードはちょっとだけ改変;ちなみにNumpyらしい書き方になっているのでかなり高速化されている).

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


def df(init=None, time=None):
    return np.cos(time)


def scipy_odeint(func, init, time):
    return odeint(func, init, time)


def euler_integration(func, init, time, dt=None):
    if dt is None:
        dt = np.diff(time).mean()
    t = np.arange(time.min(), time.max()+dt*2, dt)
    ty = func(time=t[:-1])
    y = np.concatenate(([init], np.cumsum(ty) * dt))
    return t, y


def runge_kutta(func, init, time, dt=None):
    if dt is None:
        dt = np.diff(time).mean()
    t = np.arange(time.min(), time.max()+dt*2, dt)
    s1 = func(time=t[:-1])
    s2 = func(time=t[:-1]+dt/2)
    s3 = func(time=t[:-1]+dt/2)
    s4 = func(time=t[:-1]+dt)
    ty = (s1 + 2*s2 + 2*s3 + s4) * dt / 6
    y = np.concatenate(([init], np.cumsum(ty)))
    return t, y


x = np.arange(0, np.pi*2.0, 0.1)
y = np.sin(x)
y_ode = scipy_odeint(df, 0, x)
xi1, yi1 = euler_integration(df, 0, x, 0.1)
xi2, yi2 = euler_integration(df, 0, x, 1.0)
xr1, yr1 = runge_kutta(df, 0, x, 1.0)
xr2, yr2 = runge_kutta(df, 0, x, 0.1)

fig = plt.figure(figsize=(19, 10))
ax = fig.add_subplot(111)
ax.plot(x, y, label='True')
ax.plot(x, y_ode, 'c--', label='odeint')
ax.plot(xi1, yi1, 'xb', label='euler_dx=0.1')
ax.plot(xi2, yi2, 'ob', label='euler_dx=1.0')
ax.plot(xr1, yr1, 'or', label='rk_dx=1.0')
ax.plot(xr2, yr2, 'xr', label='rk_dx=0.1')
ax.grid(True)
ax.legend()
plt.savefig('integration.png')

integration

広告
カテゴリー: 未分類 パーマリンク

ルンゲクッタとか懐かしい への1件のフィードバック

  1. ピンバック: 積分 | 粉末@それは風のように (日記)

コメントは受け付けていません。