【力学】単振り子の厳密解の初等関数による高精度な近似

振り子の解は高校で $\theta=\theta_0 \sin(\omega t +\phi)$ のようにならったかと思います。これは1次の精度での近似解です。振り子は厳密に解くこともできます。しかし、その解はヤコビのsn関数という特殊関数や周期に第一種完全楕円積分などが用いられ、やや複雑です。この記事ではこれを精度をある程度損なわず、そこそこに良い初等関数で構成される近似解を記述してみたいと思います。この記事は僕が高校の課題研究で扱っていたもので、当時は慣性抵抗と粘性抵抗までを含めた近似解を導出し、数値計算や実験との比較も行っていました。

無抵抗単振り子の厳密解

単振り子の初期条件は角 $\displaystyle \theta_0$ から静かに離すものとします。このラグランジアンは、
$\displaystyle L=\frac{1}{2}ml^2{\dot{\theta}}^2-mgl\left(1-\cos\theta\right)$
なので、オイラー-ラグランジュの運動方程式を立てると、
$\displaystyle \ddot{\theta}+\omega_{0}^2\sin\theta=0$
$\displaystyle \omega_{0}^2\equiv\sqrt{\frac{g}{l} }$
$\displaystyle dt=\frac{1}{\sqrt{2\omega_0}}\frac{d\theta}{\sqrt{\cos\theta-\cos\theta_0}}$
$\displaystyle \sin\phi\equiv\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}$
$\displaystyle \omega_{0}t=\pm\int_{\frac{\pi}{2}}^{\phi} \frac{d\phi}{\sqrt{1-\sin^2{\frac{\theta_0}{2}}\sin\phi}}$
と表されるので、積分を実行することで解が得られます。既によく知られた厳密な解析では、ヤコビのsn関数と定義することで、 $\displaystyle \theta$ の関数が表されます。ヤコビのsn関数を $\displaystyle \mathrm{sn}(u,k)$ とすると、次のように表されます。ただし、 $\displaystyle K(k)$ は第一種完全楕円積分といいます。
$\displaystyle \theta=2\sin^{-1}\left[\sin{\left(\frac{\theta_0}{2}\right)}\mathrm{sn}\left(\omega_0t+K\left(\sin{\frac{\theta_0}{2}}\right),\sin{\frac{\theta_0}{2}}\right)\right]$
$\displaystyle K(k)\equiv\int_{0}^{\frac{\pi}{2}}\frac{d\phi}{\sqrt{1-k^2\sin^2\phi}}$
そして、ヤコビのsn関数や第一種完全楕円積分は初等関数で表せないということが証明されています。(参考文献[1])

近似手法

$\displaystyle \omega_0t=\pm\int_{\frac{\pi}{2}}^{\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}}\frac{d\phi}{\sqrt{1-\sin^2{\frac{\theta_0}{2}}\sin^2{\phi}}}$
ここで、
$\displaystyle 1\gg\sin^2{\frac{\theta_0}{2}}\sin^2{\phi}$
であれば、
$\displaystyle \sqrt{1-\sin^2{\frac{\theta_0}{2}}\sin^2{\phi}}\approx1-\frac{1}{2}\sin^2{\frac{\theta_0}{2}}\sin^2{\phi}$
という近似が適用できます。
$\displaystyle \omega_0t\approx\pm\int_{\frac{\pi}{2}}^{\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}}\frac{d\phi}{1-\frac{1}{2}\sin^2{\frac{\theta_0}{2}}\sin^2{\phi}}$
$\displaystyle =\pm\int_{\frac{\pi}{2}}^{\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}}\frac{d\phi}{\cos^2{\phi}+\sin^2{\phi}-\frac{1}{2}\sin^2{\frac{\theta_0}{2}}\sin^2{\phi}}$
$\displaystyle =\pm\int_{\frac{\pi}{2}}^{\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}}\frac{d\phi}{\cos^2{\phi}+\left(1-\frac{1}{2}\sin^2{\frac{\theta_0}{2}}\right)\sin^2{\phi}}$
$\displaystyle =\pm\int_{\frac{\pi}{2}}^{\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}}{\frac{\frac{1}{\cos^2{\phi}}}{1+\left(1-\frac{1}{2}\sin^2{\frac{\theta_0}{2}}\right)\tan^2{\phi}}d\phi}$
次の変数変換をします。
$\displaystyle \sqrt{1-\frac{1}{2}\sin^2{\frac{\theta_0}{2}}}=\sqrt{\frac{\cos{\theta_0}+3}{4}}\equiv a\left(\theta_0\right)$
$\displaystyle a\left(\theta_0\right)\tan{\phi\equiv}\tan{\Phi}$
$\displaystyle d\phi{=\frac{\frac{\cos^2{\phi}}{\cos^2{\Phi}}}{a\left(\theta_0\right)}d\Phi}$
そうすれば、
$\displaystyle \omega_0t\approx\pm\int_{\frac{\pi}{2}}^{\tan^{-1}\left(a\left(\theta_0\right)\tan\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}\right)}{\frac{\frac{1}{\cos^2{\phi}}\frac{\cos^2{\phi}}{\cos^2{\Phi}}}{\left(1+\tan^2{\Phi}\right)a\left(\theta_0\right)}d\Phi}$
$\displaystyle =\pm\frac{1}{a\left(\theta_0\right)}\int_{\frac{\pi}{2}}^{\tan^{-1}\left(a\left(\theta_0\right)\tan\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}\right)}d\Phi$
$\displaystyle =\pm\frac{1}{a\left(\theta_0\right)}\left\{\tan^{-1}\left(a\left(\theta_0\right)\tan\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}\right)-\frac{\pi}{2}\right\}$
$\displaystyle a\left(\theta_0\right)\omega_0t\pm\frac{\pi}{2}\approx\pm\tan^{-1}\left(a\left(\theta_0\right)\tan\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}\right)$
ここで、
$\displaystyle \sin^{-1}\cos\left(a\left(\theta_0\right)\omega_0t\right)\approx\tan^{-1}\left(a\left(\theta_0\right)\tan\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}\right)$
とします。これによって、 $\displaystyle \pm$ を式中から無くすことができます。
$\displaystyle \tan\sin^{-1}\cos\left(a\left(\theta_0\right)\omega_0t\right)\approx a\left(\theta_0\right)\tan\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}$
$\displaystyle \tan^{-1}\left\{\frac{1}{a\left(\theta_0\right)}\tan\sin^{-1}\cos\left(a\left(\theta_0\right)\omega_0t\right)\right\}\approx\sin^{-1}\frac{\sin{\frac{\theta}{2}}}{\sin{\frac{\theta_0}{2}}}$
$\displaystyle \sin{\left(\frac{\theta_0}{2}\right)}{\rm \sin\tan}^{-1}\left\{\frac{1}{a\left(\theta_0\right)}\tan\sin^{-1}\cos\left(a\left(\theta_0\right)\omega_0t\right)\right\}\approx\sin{\frac{\theta}{2}}$
$\displaystyle \theta\approx2{\sin^{-1}}\left[\sin{\left(\frac{\theta_0}{2}\right)}{\rm \sin\tan}^{-1}\left\{\frac{1}{a\left(\theta_0\right)}\tan\sin^{-1}\cos\left(a\left(\theta_0\right)\omega_0t\right)\right\}\right]$
ただし、 $\displaystyle \tan$ が定義されない引数の点では、 $\displaystyle \theta$ が連続となるような値とします。
この振動の周期 $\displaystyle T$ は、
$\displaystyle T\approx\frac{4}{\omega_0}\int_{0}^{\frac{\pi}{2}}\frac{d\phi}{1-\frac{1}{2}\sin^2{\frac{\theta_0}{2}}\sin^2{\phi}}$
$\displaystyle =\frac{4}{a\left(\theta_0\right)\omega_0}\int_{0}^{\frac{\pi}{2}}d\phi$
$\displaystyle =\frac{2\pi}{a\left(\theta_0\right)\omega_0}$
また、第一種完全楕円積分は次のように近似されます。
$\displaystyle K\left(\sin{\frac{\theta_0}{2}}\right)\approx\frac{\pi}{2}\frac{1}{a\left(\theta_0\right)}$
この近似解における $\displaystyle a\left(\theta_0\right)$ は、周期の振依存性を補正しているとわかります。したがって、
$\displaystyle T\approx\frac{2\pi}{a\left(\theta_0\right)\omega_0}$
の精度が向上するような改良された $\displaystyle a\left(\theta_0\right)$ であれば、より精度の高い $\displaystyle \theta$ の解を表せられることが期待されます。
次に示す算術幾何平均と楕円積分の関係は第一種完全楕円積分の数値計算に用いられ、収束が早いことが知られています。
$\displaystyle I\left(a_0,b_0\right)\equiv\int_{0}^{\frac{\pi}{2}}\frac{d\vartheta}{\sqrt{{a_0}^2\cos^2\vartheta+{b_0}^2\sin^2\vartheta}}=I\left(\frac{a_0+b_0}{2},\sqrt{a_0b_0}\right)$
$\displaystyle a_{n+1}=\frac{a_n+b_n}{2},b_{n+1}=\sqrt{a_nb_n}$
$\displaystyle I\left(a_0,b_0\right)=I\left(a_\infty,a_\infty\right)=I\left(b_\infty,b_\infty\right)=\frac{1}{a_\infty}\int_{0}^{\frac{\pi}{2}}d\vartheta=\frac{\pi}{2a_\infty}$
$\displaystyle K\left(\sin{\frac{\theta_0}{2}}\right)=I\left(1, \cos\frac{\theta_0}{2}\right)\equiv\frac{\pi}{2a}$
$\displaystyle a$ は1と $\displaystyle \cos\frac{\theta_0}{2}$ の算術幾何平均とします。これを有限回で打ち切って利用します。
1回目
$\displaystyle \frac{1+\cos\frac{\theta_0}{2}}{2},\ \sqrt{\cos\frac{\theta_0}{2}}$
2回目
$\displaystyle \frac{\frac{1+\cos\frac{\theta_0}{2}}{2}+\sqrt{\cos\frac{\theta_0}{2}}}{2}=\left(\frac{1+\sqrt{\cos\frac{\theta_0}{2}}}{2}\right)^2,\sqrt{\frac{1+\cos\frac{\theta_0}{2}}{2}\sqrt{\cos\frac{\theta_0}{2}}}$
最後に、3回目の幾何平均をとると、
$\displaystyle a\approx\frac{1+\sqrt{\cos\frac{\theta_0}{2}}}{2}\sqrt{\sqrt{\frac{1+\cos\frac{\theta_0}{2}}{2}\sqrt{\cos\frac{\theta_0}{2}}}}$
さきほどの $\displaystyle a\left(\theta_0\right)$ は、
$\displaystyle a\left(\theta_0\right)\approx\frac{\pi}{2}\frac{1}{K\left(\sin{\frac{\theta_0}{2}}\right)}$
であったので、改良された $\displaystyle a\left(\theta_0\right)$ は次のように定められて、精度を改良できます。
$\displaystyle a\left(\theta_0\right)\equiv\frac{1+\sqrt{\cos\frac{\theta_0}{2}}}{2}\sqrt{\sqrt{\frac{1+\cos\frac{\theta_0}{2}}{2}\sqrt{\cos\frac{\theta_0}{2}}}}$
すなわち、ヤコビのsn関数を次のように近似していることとなります。
$\displaystyle \mathrm{sn}\left(u,k\right)=\sin(F^{-1}\left(u,k\right))$
$\displaystyle \approx{\rm \sin\tan}^{-1}\left\{\frac{1}{a\left(\theta_0\right)}\tan\sin^{-1}\sin\left(a\left(\theta_0\right)u\right)\right\}$

以下のグラフ*1は、 $\mathrm{sn}\left(u,k\right)$ とこの近似関数のプロットです。微小角で近似したものなので、母数が1に近いほど誤差は大きくなります。
$k=0.9$

$k=0.99$

$k=0.999$

$k=0.9999$

参考文献

[1] 松信伊, 理正夫, 竹内啓『初等関数の数値計算』（シリーズ新しい応用の数学（8））（1974）.

*1:

Jacobiのsn関数の初等関数のみによる近似関数

sn(x+K(k),k)≒sin(arctan(1/a*tan(arcsin(cos(a*x))))

a≡(1+sqrt(sqrt(1-k^2)))^2/4

添付画像順に母数k=0.9,0.99,0.999,0.9999 pic.twitter.com/xqLOLoJrOs
— フモフモさん (@fumofumobun) 2018年8月20日