SciPyの微分で振り子をシミュレートする

こんにちは、アドカレ四日目です。

はじめに

PythonにはSciPyという数学処理ができるライブラリがあります。
SciPyではいろいろな種類の高度な数学処理ができます。
今回はSciPyの微分のサンプルとして、振り子をシミュレートしようと思います。

振り子のシミュレーション

前提

以降のコードは全てGoogle Colaboratoryで実行することを前提としており、ここに書かれた実行結果もまた、Google Colaboratoryで実行した結果です。
Google Colaboratoryはシミュレーションの結果をアニメーションとして再生したいため利用しています。

バージョン等に関しては、全て2022/12/04(記事執筆当時)のデフォルトを使用しています。

ライブラリのインポート

使用するライブラリをインポートします。
SciPy : シミュレーション時に微分を行う
NumPy : 実際にシミュレーションによって座標を求める
Matplotlib : シミュレーション結果を描画する

import numpy as np
from numpy import sin, cos, pi
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from matplotlib import rc
from IPython.display import HTML

定数の定義

ここでは、シミュレーションに使用する定数を定義します。
DURATIONとINTERVALについては、SciPyで微分する際に利用する値で、INTERVALは微分する際の最小の微小時間で、DURATIONは微分する期間です。
微小時間についての詳しい説明は、追加する！！！！！！！で行います

# 重力加速度
G  = 9.8
# 落下開始角度
INIT_RAD = pi/6
# 落下開始速度
INIT_V = 1
# 振り子の長さ
L = 1
# 微分に使う値
DURATION = 10
INTERVAL = 0.05

振り子の計算

振り子の公式

振り子を計算する公式はラグランジュ力学から

\frac{\delta\theta}{\delta t}=\omega
\\
\frac{\delta\omega}{\delta t}=-\frac{G}{L}\sin{\theta}

と微分の式で表せます。
Qiitaだと分数の長さが短いのなんで...

SciPyで微分をするためには、このまま式を関数にします。
引数の形は(計算に使用する値, 微小時間)にします。

def pendulum(state, t):
    theta, omega = state
    dydt = np.zeros_like(state)
    dydt[0] = omega
    dydt[1] = -(G/L)*sin(theta)
    return dydt

微分

SciPyで実際に微分をします。
SciPyでの微分の表現は離散的で、細かい値ずつ変数を変えていき、毎回その値を求めることによって微分を計算しています。
微分をするにはodeint関数を使います。引数はodeint(微分する式に当たる関数, 関数に使う変数, 微小に変わる値の配列)です。
微分する式は先ほど定義したpendulumです。
関数に使う変数は、今回は[θ, ω]としています。
微小時間の配列は、少しずつ値を増加させた配列を使います。今回は、定数で指定したように、0からDURATIONまでINTERVALずつ増加する配列を用意します。

# 初期値
state = [INIT_RAD, INIT_V/L]
# 0からDURATIONまでINTERVALずつ増加する
t = np.arange(0.0, DURATION + INTERVAL, INTERVAL)
# 微分実行
solv = odeint(pendulum, state, t)

極座標をxy座標に変換

solvは各t秒における角度と速度配列です。
なので、角度theta、大きさLの極座標として、三角関数で直交座標に変換します。

theta = solv[:, 0]
x = L * sin(theta)
y = - L * cos(theta)

結果を描画

matplotlibを使用して、シミュレーション結果を描画します。
また、アニメーションとして出力します。

fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-L, L), ylim=(-L, L))
ax.grid()
 
line, = ax.plot([], [], 'ro-', animated = True)

def update(i):
    next_x = [0, x[i]]
    next_y = [0, y[i]]
    line.set_data(next_x, next_y)
    return line,

anim = FuncAnimation(fig, update, frames=np.arange(0, len(t)), interval=INTERVAL*1000, blit=True)

rc('animation', html='jshtml')
anim

まとめ

実行結果