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

はじめに

単振り子や2重振り子の運動方程式は、ラグランジアンを用いることで求めることができる。しかし、単振り子は周期的な単純な運動をする一方で、2重振り子は複雑な運動を行う。今回は、Pythonを用いて運動方程式を解析的に解くことで単振り子や2重振り子の運動の様子をシミュレーションすることを目的とする。ただし、シミュレーションに用いた運動方程式は以下のサイトを参考にしたので、導出等を知りたい方は参考にして欲しい。

単振り子

運動方程式

xy平面に質点$m_1$が存在していて、原点を支店として糸で繋げられている単振り子を考える。
重力は-y方向にかかっているとする。
質点の座標は糸の長さ$l_1$,重力方向との角度$\theta_1$を用いて以下のように考えることができる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & x_1 = l_1 sin \theta_1 \\
    & y_1 = -l_1 cos \theta_1 \\
    \end{aligned}
\right.
\end{equation}

一方で動径方向の運動方程式は以下のように表すことができる。

\frac{d^2 \theta_1}{dt^2}= -g sin \theta_1 

プログラム

これを元にプログラムを作成すると以下のようになる。

python simple_pendulum.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


n=100




# 振り子の長さ
l_1=2


#時刻の初期値
t=0
#角度の初期値と角速度の初期値
theta_1 = math.pi/2
d_theta_dt_1=0
dd_theta_dd_t_1=0

#微小時間
dt =0.1

while(t<10):
 
  dd_theta_dd_t_1 = -g*np.sin(theta_1)

 
  d_theta_dt_1 = d_theta_dt_1 + dd_theta_dd_t_1*dt
  theta_1 = theta_1 + d_theta_dt_1*dt


  t = t+dt
  
  

  x_1=l_1*np.sin(theta_1)
  y_1=-l_1*np.cos(theta_1)
  


  X_1= np.linspace(0,x_1,n)
  Y_1= np.linspace(0,y_1,n)


  im1=plt.plot(X_1,Y_1,'-')
  
  ims.append(im1)


# 10枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=100)
#保存
ani.save("furiko_simple.gif", writer="pillow")
plt.show()

このプログラムを実行すると以下のような動画が出力される。

考察

上記の動画のように単振り子は往復運動を示す。

2重振り子

運動方程式

運動方程式の導出は以下のサイトを参照して欲しい。本記事では、導出された結果を用いるのみとする。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & \ddot{\theta_1} = -\frac{(m_1+m_2)gl_1 sin\theta_1+m_2 l_1 l_2 (\ddot{\theta_2} cos(\theta_1-\theta_2)+\ddot {\theta_2^2} sin(\theta_1-\theta_2))}{(m_1+m_2)l_1^2}\\
    & \ddot{\theta_2} = -\frac{m_2 gl_1 sin\theta_2+m_2 l_1 l_2 (\ddot{\theta_1} cos(\theta_1-\theta_2)+\ddot \theta_1^2 sin(\theta_1-\theta_2))}{m_2 l_2^2} \\
    \end{aligned}
\right.
\end{equation}

一方で、質点$m_1$と質点$m_2$の位置は以下のようになる。

\begin{equation}
\left\{ \,
    \begin{aligned}
    & x_1 = l_1 sin \theta_1 \\
    & y_1 = -l_1 cos \theta_1 \\
    \end{aligned}
\right.
\end{equation}

\begin{equation}
\left\{ \,
    \begin{aligned}
    & x_2 = l_1 sin \theta_1+ l_2 sin \theta_2 \\
    & y_2 = -l_1 cos \theta_1- l_2 cos \theta_2 \\
    \end{aligned}
\right.
\end{equation}

プログラム

以上のことに留意してプログラムを書くと以下のようになる。

python double_pendulum.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


n=100

m_1=1.0
m_2=2.0

# 振り子の長さ
l_1=2
l_2=1

#時刻の初期値
t=0
#角度の初期値と角速度の初期値
theta_1 = math.pi/2
d_theta_dt_1=0
dd_theta_dd_t_1=0

theta_2 = 0
d_theta_dt_2=0
dd_theta_dd_t_2=0

#微小時間
dt =0.01

while(t<10):
 
  dd_theta_dd_t_1 = -((m_1+m_2)*g*l_1*np.sin(theta_1)+m_2*l_1*l_2*(dd_theta_dd_t_2*np.cos(theta_1-theta_2)+((d_theta_dt_2)**2)*np.sin(theta_1-theta_2)))/((m_1+m_2)*l_1**2)

  dd_theta_dd_t_2 = -(m_2*g*l_2*np.sin(theta_2)+m_2*l_1*l_2*(dd_theta_dd_t_1*np.cos(theta_1-theta_2)-((d_theta_dt_1)**2)*np.sin(theta_1-theta_2)))/(m_2*l_2**2)


  d_theta_dt_1 = d_theta_dt_1 + dd_theta_dd_t_1*dt
  theta_1 = theta_1 + d_theta_dt_1*dt

  d_theta_dt_2 = d_theta_dt_2 + dd_theta_dd_t_2*dt
  theta_2 = theta_2 + d_theta_dt_2*dt

  t = t+dt
  
  

  x_1=l_1*np.sin(theta_1)
  y_1=-l_1*np.cos(theta_1)
  
  x_2= l_1*np.sin(theta_1)+l_2*np.sin(theta_2)
  y_2= -l_1*np.cos(theta_1)-l_2*np.cos(theta_2)

  X_1= np.linspace(0,x_1,n)
  Y_1= np.linspace(0,y_1,n)

  X_2= np.linspace(x_1,x_2,n)
  Y_2= np.linspace(y_1,y_2,n)

  im1=plt.plot(X_1,Y_1,'-')
  im2=plt.plot(X_2,Y_2,'-')
  ims.append(im1+im2)


# 10枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=100)
#保存
ani.save("furiko_2tai_simple.gif", writer="pillow")
plt.show()

考察

このプログラムを実行すると以下のようになる。

このように、2重振り子はカオスな運動となる。しかし、やはり運動の様子を正確に記述するためには十分な計算量が必要となる。

まとめ

今回は、単振り子や2重振り子の運動方程式からPythonを用いて振り子のシミュレーションを行うことを目的としてプログラムを作成した。結果、微分方程式を解くために用いた差分法によるアルゴリズムは上手く動作しているようで、シミュレーションをすることができた。