はじめに
単振り子の運動方程式は未知の張力を考慮しなくて済む極座標によって記述されることがほとんどである。そのときに、三角関数が多用される。しかし、直交座標による運動方程式の叙述は、未知の力である張力を考慮に入れなければならないため複雑になる。そこで、今回は単振り子の運動方程式に直交座標系を導入し、Pythonの差分法で単振り子の運動をシミュレーションすることを目的とする。
条件設定
-y方向に重力が働いているxy座標系の原点に糸が吊り下がっていて末端におもり$m$が存在しているものとする。このとき、おもりに対する運動方程式は以下のように表すことができる。ただし、糸の長さを$R$、張力を$T$とする。$T$は未知数である。
\begin{equation}
\left\{ \,
\begin{aligned}
& a_x = -T \frac{x}{mR} \\
& a_y = -g-T \frac{y}{mR}\\
\end{aligned}
\right.
\end{equation}
一方で、円の方程式より以下のことが成立する。
y=-\sqrt{R^2-x^2} (y<0)
一方で以下の等式も成立する。
x^2+y^2=R^2
ゆえに、これを$t$で微分すると、
2x\ddot{x}+2y\ddot{y}=0
整理すると、
x\ddot{x}+y\ddot{y}=0
もう一度$t$で微分すると、
(\dot{x}^2+x\ddot{x})+(\dot{y}^2+y\ddot{y})=0
したがって、
\ddot{x}=-\frac{\dot{x}^2+\dot{y}^2+y\ddot{y}}{x}(=-T \frac{x}{mR})
以上の考察により、張力を位置、速度、加速度で表すと以下のようになる。
T=\frac{mR(v_x^2+v_y^2+y a_y)}{x^2}
これにより、おもりの直交座標表示に関する運動方程式が叙述されることとなる。
プログラム
ゆえに、以上のことから、差分法を用いて運動方程式を数値解析するプログラムを作成した。
python simple_pendulum_xy.py
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import japanize_matplotlib
import math
# アニメ
fig = plt.figure()
ims = []
# パラメータと定数
g=9.8
# 振り子の長さ
R=1.0
# 振り子の質量
m=1.0
#時刻の初期値
t=0
# 位置の初期値
x=0.99
# 速度の初期値
v_x=0.0
v_y=0.0
# 加速度の初期値
a_x=0.0
a_y=0.0
# 位置、速度、加速度のリスト
X=[]
Y=[]
V_X=[]
V_Y=[]
A_X=[]
A_Y=[]
#微小時間
dt =0.01
while(t<10):
# 束縛条件
y=-(R**2-x**2)**0.5
# 張力の大きさ
T=m*R*(v_x**2+v_y**2+y*a_y)/(x**2)
# 張力の大きさの上限
if T>10:
T=10
# 運動方程式
a_x= -T*x/(m*R)
a_y=-g-T*y/(m*R)
# 速度の更新
v_x=v_x+a_x*dt
v_y=v_y+a_y*dt
# 位置の更新
x = x + v_x*dt
y= y + v_y*dt
A_X.append(a_x)
A_Y.append(a_y)
V_X.append(v_x)
V_Y.append(v_y)
X.append(x)
Y.append(y)
# 時刻を更新
t = t + dt
# プロット
im = plt.plot([0,x],[0,y],marker='o',color='b')
ims.append(im)
# 10枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=100)
#保存
ani.save("simple_pendulum_xy.gif", writer="pillow")
plt.show()
これを実行すると以下のようなgifが出力される。
このように、上手く単振り子の振動を表現することができている。しかし、課題として$x=0$付近だと張力$T$が発散してしまうため、上限を設けることで簡易的に対応した。
まとめ
今回は、単振り子の運動方程式をxy座標系で叙述し、差分法によりPythonで数値解析し、シミュレーションすることを試みた。おもりが原点の真下にあるとき張力が発散するので、上限を設けて対処した。結果、単振り子の振動の様子を上手くシミュレーションすることができた。
参考文献