トップ > プログラミング > Python(SciPy)で単振り子

Python(SciPy)で単振り子

プログラミング プログラミング-Python 物理学 物理学-力学

POINT

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

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

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

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

運動方程式

単振り子の運動方程式を立て，無次元化を行います．ここではかいつまんだ説明をします．詳細を知りたい場合は次の記事を参照してください．：無次元化が必要な理由と方法〜数値計算の疑問 - Notes_JP

単振り子の運動方程式は

\begin{aligned}
m\frac{\mathrm{d}^2}{\mathrm{d}t^2} (l\theta)
=-mg\sin\theta
\end{aligned}
で与えられます．

$\tau = t\sqrt{g/l}$として無次元化すると

\begin{aligned}
\frac{\mathrm{d}^2\theta}{\mathrm{d}\tau^2}
=-\sin\theta
\end{aligned}
となります．以下，$\dot{\theta}=\mathrm{d}\theta/\mathrm{d}\tau$, $\ddot{\theta}=\mathrm{d}^2\theta/\mathrm{d}\tau^2$と表すことにします．

厳密解

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

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

すると，上のエネルギー保存則から

\begin{aligned}
\frac{1}{2}\dot{\theta}^2 -\cos\theta = - \cos\theta_{\mathrm{max}}
\end{aligned}
なので，$\dot{\theta}>0$の範囲を考えれば
\begin{aligned}
\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{aligned}
がわかります（$\cos\theta=1-2\sin^2(\theta/2)$）．ここで，$\sin\xi=\sin(\theta/2)\bigl/\sin(\theta_{\mathrm{max}}/2)$と変数変換をします．両辺を$\xi$で微分することで得られる
\begin{aligned}
\cos\xi
=\cos(\theta/2)\cdot \frac{1}{2} \frac{\mathrm{d}\theta}{\mathrm{d}\xi}
\biggl/\sin(\theta_{\mathrm{max}}/2)
\end{aligned}
から
\begin{aligned}
\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{aligned}
です．したがって，$\theta(0)=0$とするとき
\begin{aligned}
\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{aligned}
となります．ここで，$F$は第1種不完全楕円積分と呼ばれます．

楕円積分$F$は，楕円関数$\mathrm{sn\,}$と$\mathrm{sn\,} \bigl(\tau, \sin(\theta_{\mathrm{max}}/2)\bigr)=\sin\xi$の関係があるので，厳密解が

\begin{aligned}
\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{aligned}
であることがわかりました．

【補足】
振動の周期を$T$とすれば$\theta(T/4)=\theta_{\mathrm{max}}$なので，

\begin{aligned}
\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{aligned}
となります．$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{aligned}
\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{aligned}
となります（$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{aligned}
\dot{\theta}
= \mathrm{sgn}\bigl(\dot{\theta}\bigr) \sqrt{ 2(\cos\theta - \cos\theta_{\mathrm{max}}) }
\end{aligned}
です．ここで， $\mathrm{sgn}$は
\begin{aligned}
\mathrm{sgn} (x)
&=
\begin{cases}
\,+1 &(x>0)\\
\,0 &(x=0) \\
\,-1 &(x<0)
\end{cases}
\end{aligned}
で定義される符号関数です．初期条件$\theta(0)=\theta_{\mathrm{max}}$の運動では，周期性から1/2周期ごとに$\dot{\theta}$の符号が変わることから
\begin{aligned}
\mathrm{sgn}\bigl(\dot{\theta}\bigr)
&=\mathrm{sgn} \bigl[ \mathrm{mod}(t,T) - T/2 \bigr]
\end{aligned}
と計算することができます（$T$は周期）．ここで$\mathrm{mod}(t,T)$は，時刻$t$を周期$T$で割ったときの剰余です．

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

単振り子の数値解(odeint)と厳密解
単振り子の数値解(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.5.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.5.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()