Python(SciPy)で単振り子

POINT

  • 単振り子の厳密解とPython(SciPy)の計算結果を比較する.
  • 厳密解の導出を解説する.

数値計算の妥当性を確認するために,2通りの方法

  1. 常微分方程式をSciPy(odeint, ode, solve_ivp)を用いて解いたものをプロットする方法
  2. 厳密解を楕円積分・楕円関数を用いてプロットする方法

で同じ結果が得られることを確認します.

運動方程式

単振り子の運動方程式を立て,無次元化を行います.ここではかいつまんだ説明をします.詳細を知りたい場合は次の記事を参照してください.:

単振り子の運動方程式は
\begin{align}
m\frac{\mathrm{d}^2}{\mathrm{d}t^2} (l\theta)
=-mg\sin\theta
\end{align}で与えられます.

$\tau = t\sqrt{g/l}$として無次元化すると
\begin{align}
\frac{\mathrm{d}^2\theta}{\mathrm{d}\tau^2}
=-\sin\theta
\end{align}となります.以下,$\dot{\theta}=\mathrm{d}\theta/\mathrm{d}\tau$, $\ddot{\theta}=\mathrm{d}^2\theta/\mathrm{d}\tau^2$と表すことにします.

厳密解

運動方程式の両辺に$\dot{\theta}$をかけると,
\begin{align}
\dot{\theta}\ddot{\theta}=-\dot{\theta}\sin\theta
\end{align}となります.したがって,
\begin{align}
\frac{\mathrm{d}}{\mathrm{d}\tau}\biggl( \frac{1}{2}\dot{\theta}^2 -\cos\theta\biggr)=0
\end{align}です.これは,エネルギー保存則を表しています.

$-\pi<\theta<\pi$の範囲の運動を考え,$\theta$の最大値を$\theta_{\mathrm{max}}$としましょう.このとき,振り子の速度は$0$つまり,$\dot{\theta}=0$となるはずです.

すると,上のエネルギー保存則から
\begin{align}
\frac{1}{2}\dot{\theta}^2 -\cos\theta = - \cos\theta_{\mathrm{max}}
\end{align}なので,$\dot{\theta}>0$の範囲を考えれば
\begin{align}
\frac{\mathrm{d}\tau}{\mathrm{d}\theta}
&=\frac{1}{ \sqrt{2(\cos\theta - \cos\theta_{\mathrm{max}} )} } \\
&=\frac{1}{2\sin(\theta_{\mathrm{max}}/2)}
\frac{1}{\sqrt{1-\biggl[\dfrac{\sin(\theta/2)}{\sin(\theta_{\mathrm{max}}/2)} \biggr]^2 }}
\end{align}がわかります($\cos\theta=1-2\sin^2(\theta/2)$).ここで,$\sin\xi=\sin(\theta/2)\bigl/\sin(\theta_{\mathrm{max}}/2)$と変数変換をします.両辺を$\xi$で微分することで得られる
\begin{align}
\cos\xi
=\cos(\theta/2)\cdot \frac{1}{2} \frac{\mathrm{d}\theta}{\mathrm{d}\xi}
\biggl/\sin(\theta_{\mathrm{max}}/2)
\end{align}から
\begin{align}
\frac{\mathrm{d}\theta}{\mathrm{d}\xi}
&=\frac{2 \sin(\theta_{\mathrm{max}}/2)\cos\xi }{\cos(\theta/2)} \\
&=\frac{2 \sin(\theta_{\mathrm{max}}/2)\sqrt{1-\sin^2\xi} }{\sqrt{1-\sin^2(\theta_{\mathrm{max}}/2)\sin^2\xi}}
\end{align}です.したがって,$\theta(0)=0$とするとき
\begin{align}
\tau
&=\int_0^\theta \frac{\mathrm{d}\tau}{\mathrm{d}\theta} \,\mathrm{d}\theta\\
&=\int_0^\xi \frac{\mathrm{d}\xi}{\sqrt{1-\sin^2(\theta_{\mathrm{max}}/2)\sin^2\xi}} \\
&=F\bigl(\xi, \sin(\theta_{\mathrm{max}}/2) \bigr)
\end{align}となります.ここで,$F$は第1種不完全楕円積分と呼ばれます.

楕円積分$F$は,楕円関数$\mathrm{sn\,}$と$\mathrm{sn\,} \bigl(\tau, \sin(\theta_{\mathrm{max}}/2)\bigr)=\sin\xi$の関係があるので,厳密解が
\begin{align}
\theta(\tau)
&=2\arcsin\Bigl[\sin(\theta_{\mathrm{max}}/2) \\
&\qquad\qquad\qquad \times\mathrm{sn} \bigl(\tau, \sin(\theta_{\mathrm{max}}/2)\bigr) \Bigr]
\end{align}であることがわかりました.

【補足】
振動の周期を$T$とすれば$\theta(T/4)=\theta_{\mathrm{max}}$なので,
\begin{align}
\frac{T}{4}
&=\int_0^{\theta_{\mathrm{max}}} \frac{\mathrm{d}\tau}{\mathrm{d}\theta} \,\mathrm{d}\theta\\
&=\int_0^{\pi/2} \frac{\mathrm{d}\xi}{\sqrt{1-\sin^2(\theta_{\mathrm{max}}/2)\sin^2\xi}} \\
&=F\bigl(\pi/2, \sin(\theta_{\mathrm{max}}/2) \bigr) \\
&=K\bigl(\sin(\theta_{\mathrm{max}}/2) \bigr)
\end{align}となります.$K$は第1種完全楕円積分と呼ばれます.

SciPyで解く

odeintを使って単振り子の常微分方程式を解くサンプルコードが公開されています.今回は,これをそのまま用います.パラメータを$b=0$, $c=1$としたものが,今回考える運動方程式です.


数値計算では$\theta(0)=\theta_{\mathrm{max}}$の状態から計算を始め,$\dot{\theta}<0$の運動($\theta$が小さくなっていく)を見ています.一方で,上で求めた厳密解は,$\theta(0)=0$かつ$\dot{\theta}>0$のものなので,比較できるように表式を書き換える必要があります.これには,

  1. $\dot{\theta}<0$の運動領域をみる
  2. 時間の原点を1/4周期ずらす

とすればよく,その表式は
\begin{align}
\theta(\tau)
&=2\arcsin\Bigl\{\sin(\theta_{\mathrm{max}}/2) \\
&\qquad \times\mathrm{sn} \bigl[-(\tau-T/4), \sin(\theta_{\mathrm{max}}/2)\bigr] \Bigr\}
\end{align}となります($T$は周期).上で見たように,周期$T$は完全楕円積分で表され,値は次を利用すれば取得できます.

また,ヤコビの楕円関数$\mathrm{sn}$は以下で求めることができます.


計算結果をプロットしたものが以下です.odeintを使って求めた$\theta$がtheta(t)_odeint,厳密解をプロットしたものがtheta(t)_exactです.

単振り子の数値解(odeint)と厳密解
単振り子の数値解(SciPy)と厳密解


次に,$\omega=\dot{\theta}$の厳密解について考えます.上のエネルギー保存則$\frac{1}{2}\dot{\theta}^2 -\cos\theta = - \cos\theta_{\mathrm{max}}$から
\begin{align}
\dot{\theta}
= \mathrm{sgn}\bigl(\dot{\theta}\bigr) \sqrt{ 2(\cos\theta - \cos\theta_{\mathrm{max}}) }
\end{align}です.ここで, $\mathrm{sgn}$は
\begin{align}
\mathrm{sgn} (x)
&=
\begin{cases}
\,+1 &(x>0)\\
\,0 &(x=0) \\
\,-1 &(x<0)
\end{cases}
\end{align}で定義される符号関数です.初期条件$\theta(0)=\theta_{\mathrm{max}}$の運動では,周期性から1/2周期ごとに$\dot{\theta}$の符号が変わることから
\begin{align}
\mathrm{sgn}\bigl(\dot{\theta}\bigr)
&=\mathrm{sgn} \bigl[ \mathrm{mod}(t,T) - T/2 \bigr]
\end{align}と計算することができます($T$は周期).ここで$\mathrm{mod}(t,T)$は,時刻$t$を周期$T$で割ったときの剰余です.

この計算により得られた厳密解と,odeintによる計算とを比較した結果が以下です:

単振り子の数値解(odeint)と厳密解
単振り子の数値解(odeint)と厳密解

参考文献/記事

  • $\S11$に単振り子の厳密解(周期)の説明があります.
    力学 (増訂第3版)   ランダウ=リフシッツ理論物理学教程

    力学 (増訂第3版) ランダウ=リフシッツ理論物理学教程

    • 作者: エリ・ランダウ,イェ・エム・リフシッツ,広重徹,水戸巌
    • 出版社/メーカー: 東京図書
    • 発売日: 1986/04
    • メディア: 単行本
    • 購入: 4人 クリック: 42回
    • この商品を含むブログ (22件) を見る
  • odeintのドキュメントでは,新しいコードはodeintではなくsolve_ivpを使うように書かれています.今回は,サンプルコードをそのまま使うことができたので,odeintを使用しました.


付録:コード

厳密解, odeint, ode, solve_ivpでグラフを作成するコードを載せておきます.すべての手法を比較した図は以下です:
単振り子(すべての手法を比較)
単振り子(すべての手法を比較)

厳密解

上で述べたことを忠実にコード化すれば厳密解が計算できます.関数はnumpyに用意されているので,比較的簡単にできます.実行すると,次のグラフが出力されます.
単振り子(厳密解)
単振り子(厳密解)

import numpy as np
from scipy.special import ellipj, ellipk
import matplotlib.pyplot as plt

#  initial conditions: theta(0) = pi - 0.1, omega(0) = 0
theta0 = np.pi - 0.1
# omega0 = 0.0
# y0 = [theta0, omega0]

# create a time array from t0...t_max/dt sampled at dt second steps
# https://matplotlib.org/examples/animation/double_pendulum_animated.html
t_min = 0
t_max = 20
dt = 0.1
t = np.arange(t_min, t_max, dt)

# exact solution
k = np.sin(theta0/2)
m = k**2
period = 4*ellipk(m)
print('T=', period)

t2 = t - ellipk(m)

sn = ellipj(-t2,m)[0]
exact_theta = 2*np.arcsin(k*sn)
exact_theta[0] = theta0 # 計算誤差でsqrtの中身がマイナスとなることの対策
sign = np.sign(np.mod(t,period) - period/2)
exact_omega = sign * np.sqrt( 2*(np.cos(exact_theta) - np.cos(theta0)) )

# plot
# https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.pyplot.plot.html
plt.plot(t, exact_theta, '.b', markersize=5, label='theta(t)_exact')
plt.plot(t, exact_omega, '.g', markersize=5, label='omega(t)_exact')
plt.legend(loc='best')

plt.rcParams['text.usetex'] = True
plt.xlabel(r'$\tau$')
plt.ylabel(r'$\theta$')

plt.grid()

plt.savefig('exact.jpg', bbox_inches='tight')

plt.show()

odeint

上で触れたとおり,scipy.integrate.odeint — SciPy v1.3.0 Reference Guideのサンプルコードを順に並べればよいだけです.次の図が出力されます:
単振り子(odeint)
単振り子(odeint)

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html
# theta''(t) + b*theta'(t) + c*sin(theta(t)) = 0

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

# theta'(t) = omega(t)
# omega'(t) = -b*omega(t) - c*sin(theta(t))
# Let y be the vector [theta, omega]
def pend_odeint(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

# parameter    
b = 0
c = 1

#  initial conditions: theta(0) = pi - 0.1, omega(0) = 0
theta0 = np.pi - 0.1
omega0 = 0.0
y0 = [theta0, omega0]     

# create a time array from t0...t_max/dt sampled at dt second steps
# https://matplotlib.org/examples/animation/double_pendulum_animated.html
t_min = 0
t_max = 20
dt = 0.1
t = np.arange(t_min, t_max, dt)
          
# odeint
sol_odeint = odeint(pend_odeint, y0, t, args=(b, c))

# plot
# https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.pyplot.plot.html
plt.plot(t, sol_odeint[:, 0], '.b', markersize=5, label='theta(t)_odeint')
plt.plot(t, sol_odeint[:, 1], '.g', markersize=5, label='omega(t)_odeint')
plt.legend(loc='best')

plt.rcParams['text.usetex'] = True
plt.xlabel(r'$\tau$')
plt.ylabel(r'$\theta$')

plt.grid()

plt.savefig('odeint.jpg', bbox_inches='tight')

plt.show()

ode

ドキュメントscipy.integrate.ode — SciPy v1.3.0 Reference GuideのExamplesを参考に,手探りで作成しています.実行結果は以下です:
単振り子(ode)
単振り子(ode)
【注意】
  • 解の出力〜プロットまでの方法が,ベストな方法かわかっていません.とりあえず,odeintと同じデータ構造になるようにしました.
  • 実行時に以下の警告が出ます.内容については,まだ理解できていません...

    UserWarning: dop853: step size becomes too small
    self.messages.get(istate, unexpected_istate_msg)))


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

# theta'(t) = omega(t)
# omega'(t) = -b*omega(t) - c*sin(theta(t))
# Let y be the vector [theta, omega]
def pend_ode(t, y, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt
    
# parameter    
b = 0
c = 1

#  initial conditions: theta(0) = pi - 0.1, omega(0) = 0
theta0 = np.pi - 0.1
omega0 = 0.0
y0 = [theta0, omega0]

# create a time array from t0..t_max/dt sampled at dt second steps
# https://matplotlib.org/examples/animation/double_pendulum_animated.html
t_min = 0
t_max = 20
dt = 0.1
t = np.arange(t_min, t_max, dt)

# ode
r = ode(pend_ode).set_integrator('dop853')
r.set_initial_value(y0, t_min).set_f_params(b, c)

for item in t:
    if item == t_min:
        sol_ode = r.integrate(item)
    else:
        sol_ode = np.vstack( (sol_ode, r.integrate(item)) )

# plot
# https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.pyplot.plot.html
plt.plot(t, sol_ode[:, 0], '.b', markersize=5, label='theta(t)_ode')
plt.plot(t, sol_ode[:, 1], '.g', markersize=5, label='omega(t)_ode')
plt.legend(loc='best')

plt.rcParams['text.usetex'] = True
plt.xlabel(r'$\tau$')

plt.grid()

plt.savefig('ode.jpg', bbox_inches='tight')

plt.show()

solve_ivp

いくつか注意が必要です.

出力結果は以下:

単振り子(solve_ivp)
単振り子(solve_ivp)

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

# theta'(t) = omega(t)
# omega'(t) = -b*omega(t) - c*sin(theta(t))
# Let y be the vector [theta, omega]
def pend_ode(t, y, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

# https://stackoverflow.com/questions/48245765/pass-args-for-solve-ivp-new-scipy-ode-api    
def pend_solve_ivp(t, y):
    return pend_ode(t, y, b, c)

# parameter    
b = 0
c = 1

#  initial conditions: theta(0) = pi - 0.1, omega(0) = 0
theta0 = np.pi - 0.1
omega0 = 0.0
y0 = [theta0, omega0]

# create a time array from t0..t_max/dt sampled at dt second steps
# https://matplotlib.org/examples/animation/double_pendulum_animated.html
t_min = 0
t_max = 20
dt = 0.1
t = np.arange(t_min, t_max, dt)

t_span = (t_min,t_max)

# solve_ivp
# https://stackoverflow.com/questions/56153628/using-scipys-solve-ivp-to-solve-non-linear-pendulum-motion
sol_solve_ivp = solve_ivp(pend_solve_ivp, t_span, y0, method='RK45', t_eval=t, rtol=1e-8)

# plot
# https://matplotlib.org/3.1.0/api/_as_gen/matplotlib.pyplot.plot.html
plt.plot(sol_solve_ivp.t, sol_solve_ivp.y[0], '.b', markersize=5, label='theta(t)_solve_ivp')
plt.plot(sol_solve_ivp.t, sol_solve_ivp.y[1], '.g', markersize=5, label='omega(t)_solve_ivp')
plt.legend(loc='best')

plt.rcParams['text.usetex'] = True
plt.xlabel(r'$\tau$')

plt.grid()

plt.savefig('solve_ivp.jpg', bbox_inches='tight')

plt.show()