バネマス系の運動をオイラー法で計算

$dy/dx = f(x,y)$の常微分方程式から解析的に$y(x)$が解ければ、どんな$x$における$y(x)$でも計算できますが、なかなかそうはいきません。運動方程式も微分方程式で表されているので、このシステムはどんな振る舞いをするの？という疑問に答えるにはやはり解析解が必要です。でもやっぱり解析解を求めることが難しいので、そんなときに役立つツールの一つが数値計算のオイラー法です。

オイラー法をバネ-マス系の運動方程式に適用してシステムの振る舞いを数値計算して、解析解との比較を行ってみます。

オイラー法とは
システムの運動方程式
運動方程式にオイラー法を適用
オイラー法で計算
解析解とオイラー法の比較
1. 計算誤差について

オイラー法とは

変数が一つである微分法手式を常微分方程式と呼びます。変数が$x$である常微分方程式 $dy(x)/dx = f(x,y(x))$の解が$y(x)$であるとします。

このとき、この$y(x)$の少し先の点$x+h$について表すと、テイラー展開を用いて(1)式で表されます。

$$
y(x+h) = y(x) + h \frac{d y(x) } {dx} + \frac{h^2}{2!} \frac{d^2 y(x)}{dx^2} + \cdots \tag{1}
$$

$dy(x)/dx = f(x,y(x))$の関係を(1)式に当てはめると(2)式になります。

$$
y(x+h) = y(x) + h f(x,y(x)) + \frac{h^2}{2!} \frac{d}{dx}f(x,y(x)) + \cdots \tag{2}
$$

(2)式から、(3)式のように$f(x,y(x))$の１階微分以上を利用しないで解$y(x)$の次の点$y(x+h)$を計算することをオイラー法といいます。

$$
y(x+h) = y(x) + h f(x,y(x)) \tag{3}
$$

一階微分以上を無視しているので、数値計算の誤差が大きいことが特徴です。ではこれを運動方程式に適用してみます。

このように$y(x+h)$の計算を$y(x)$と$x$だけで行う方法を一段式といいます。ちなみに $y(x+h)$の計算を$y(x), y(x-h), y(x – 2h), \cdots$と$x$を使い計算する方法は多段式です。

システムの運動方程式

オイラー法で計算する運動方程式は、減衰項が無いバネ-マス系に振動を加えたシステムです。

このシステムの運動方程式はニュートンの運動方程式から(4)式となります。

$$
f(t) = m \frac{d^2 x(t)}{d t^2} + k x(t) \tag{4}
$$

運動方程式にオイラー法を適用

オイラー法は$dy(x)/dx = f(x,y(x))$の形の式に適用するので(5)式を定義します。

$$
\frac{d x(t) }{ d t} = v(t) \tag{5}
$$

すると(4)式は(6)式になります。

$$
f(t) = m \frac{d v(t)}{d t} + k x(t) \tag{6}
$$
$dy(x)/dx = f(x,y(x))$ の形で書くと(7)式のようになります。

$$
\begin{eqnarray}
\frac{d x(t) }{ d t} &=& v(t) \\
\frac{d v(t)}{d t} &=& \frac{1}{m} \big(f(t)-kx(t) \big)
\end{eqnarray} \tag{7}
$$

(3)式のオイラー法を適用すると、(8)式となります。

$$
\begin{eqnarray}
x(t+h) &=& x(t) + h v(t) \\
v(t + h) &=& v(t) + \frac{h}{m} \big(f(t) -kx(t+h) \big)
\end{eqnarray} \tag{8}
$$

オイラー法で計算

パラメータを以下のようにしてPythonで計算させてその運動波形を確認します。外力$f(t)$は$F_0 \sin (\omega )t$の周期的な力です。

$x(0) [m]$	0.1
$v(0) [m/s]$	0
$k [N/m]$	5
$m [kg]$	0.1
$F_0 [N]$	2
$\omega [rad] $	12.56

import matplotlib.pyplot as plt
from math import sin,pi

# パラメータ
v0 = 0
x0 = 0.1

m = 0.1 # [kg]
k = 5 # [N/m]
F0 = 2 # [N]
f = 2 # [Hz]

omega = 2*pi*f

# 刻み幅
h = 0.02 # [sec]

x_array = [x0]
v_array = [v0]
t_array = [0]

#オイラー法で計算
for i in range(1,500):
    t = i * h
    
    x_next = x_array[i-1] + h * v_array[i-1]
    v_next = v_array[i-1] + h/m * (F0* sin(omega*t) - k * x_next)
    
    t_array.append(t)
    x_array.append(x_next)
    v_array.append(v_next)
    

#図示
plt.figure(figsize=(10,4))
plt.plot(t_array, x_array)
plt.xlabel("time[sec]")
plt.ylabel(r"$x(t)$[m]")
plt.grid(True)
plt.xticks([i for i in range(11)])
plt.show()

それっぽい計算結果になりましたね。ではこの結果を解析解と比べてみましょう。

解析解とオイラー法の比較

バネ-マス系の解析解は(9)式でした。

$$
x(t) = x_0 \cos \omega_0 t + \bigg(\frac{v_0}{\omega_o} – \frac{\omega}{\omega_0}\frac{F_0/m}{\omega_0^2 – \omega^2}\bigg) \sin \omega_0 t + \frac{F_0/m}{\omega_0^2 – \omega^2} \sin \omega t \tag{9}
$$

この式で計算した値とオイラー法の計算結果を一緒に描くと以下の図のように、ほとんど一致しています。青線がオイラー法、オレンジ線が解析解です。