feji.me

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

October 27, 2020

単振り子のシミュレーションを行う。

ld2θdt2=gsin(θ)(1)l\frac{d^2\theta}{dt^2} = -g\sin(\theta) \tag{1}

前回とは異なり、運動方程式を無次元化するため、まずは各変数の次元を書き表すと下の表のようになることがわかる。

変数 次元
ll LL
θ\theta なし
tt TT
gg LT2LT^{-2}

ここではLLの基本単位をllとする。するとggよりT=LgT = \sqrt{\frac{L}{g}}からTTの基本単位はlg\sqrt{\frac{l}{g}}となる。θ\thetaは無次元量である。

無次元化した量をそれぞれチルダ記号を使って表記すると下のようになる。

l~=l/l=1t~=t/l/g(2)\tilde{l}=l/l = 1\\ \tag{2} \tilde{t}=t/\sqrt{l/g}\\

以上の2式を(1)式に代入する。

ll~d2θdt~2=m~gsin(θ)ll~gld2θdt~2=gsin(θ)d2θdt~2=sin(θ)(3)l\tilde{l}\frac{d^2\theta}{d\tilde{t}^2} = -\tilde{m}g\sin(\theta)\\ \tag{3} l\tilde{l}\frac{g}{l}\frac{d^2\theta}{d\tilde{t}^2} = -g\sin(\theta)\\ \frac{d^2\theta}{d\tilde{t}^2} = -\sin(\theta)

よって運動方程式を無次元化できた。改めて書き直しておく。

d2θdt~2=sin(θ)(4)\frac{d^2\theta}{d\tilde{t}^2} = -\sin(\theta) \tag{4}

さて、本題はこれからで(4)式をもとにシミュレーションを行う。

環境はPhyton3.8モジュールscipy

(4)式は二階常微分方程式なので数値計算を行うには一階常微分方程式を2つ用意し、それらを連立させて解いていく流れだ。私も詳しくは知らない。

なのでここでは以下のようにおけば

dθdt~=ωd2θdt~2=dωdt~=sin(θ)(5)\frac{d\theta}{d\tilde{t}}=\omega\\ \tag{5} \frac{d^2\theta}{d\tilde{t}^2} = \frac{d\omega}{d\tilde{t}} = -\sin(\theta)

となる。この2式を連立させる。

まずは具体的なコードを示す。

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def pend(t, y):
    theta, omega = y
    dydt = [omega, -np.sin(theta)]
    return dydt

y0 = (np.pi/3, 0.0)
sol = solve_ivp(pend, (0,10), y0, t_eval=t)
plt.plot(sol.t, sol.y[0], 'b', label=r'$\theta(t)$')
plt.legend(loc='best')
plt.xlabel(r'$\tilde{t}=t\ /\ \overline{t}$', fontsize='20')
plt.ylabel(r'$\theta$', fontsize='20')
plt.grid()
plt.show()

初期条件はθ(0)=π/3,ω=0\theta(0) = \pi/3, \omega=0である。

single pengulum