表計算ソフトで１次元調和振動子を解くためのスプレッドシート [ネコ騙し数学]
https://docs.google.com/spreadsheets/d/1X8VHH3I08SrWk7EK86SCBOUh3erGy6xJ2fLUvq4G_CQ/pubhtml
https://docs.google.com/spreadsheets/d/1X8VHH3I08SrWk7EK86SCBOUh3erGy6xJ2fLUvq4G_CQ/pub?output=ods
https://docs.google.com/spreadsheets/d/1JDmMdnqOV5peAmmOVgFGZjtWwTu6aoSmiMocmmM9oFM/pubhtml
https://docs.google.com/spreadsheets/d/1JDmMdnqOV5peAmmOVgFGZjtWwTu6aoSmiMocmmM9oFM/pub?output=ods
ネムネコがオレのPCに勝手にスプレッドシートをダウンロードした
と文句をつけられるのは嫌なので、スプレッドシートのダウンロードには注意するにゃ。
クリックするとダウンロードのフォルダーにods形式のスプレッドシートがダウンロードされるケロ。これをダウンロードすると、どのように計算しているかが分かると思うにゃ。
オイラー法とシンプレクティック法を用いた１次元調和振動子の解法 [ネコ騙し数学]
オイラー法とシンプレクティック法を用いた１次元調和振動子の解法
１次元調和振動子のハミルトニアン（力学的なエネルギー）は
である。
ここで、
であり、qはバネの自然長からの変位、pは運動量、（１）式のp²/2は運動エネルギー、q²/2はバネの位置エネルギー（ポテンシャルエネルギー）に相当する。
（２）式は、
と書き換えることが可能で、（２）式とバネの単振動に関するニュートンの運動方程式と同一のものである。
時刻t=0における初期条件を(q,p)=(1,0)とすれば、（２）、（３）式の解は
である。
１次精度のオイラー法を用いて解くには、
つまり、
と近似し、これを逐次繰り返して計算すればよい。
たとえば、t=0のとき、(q,p)=(1,0)であれば、Δt=0.1とすれば、
この計算結果をもとに、
といった具合に計算できる。
この計算は非常にシンプルな繰り返し計算だから、プログラムを作らずとも、表計算ソフトで簡単に計算が可能で、以下にt=0における初期条件を(q,p)=(1,0)としたものの計算結果を示す。
このように非常に粗い近似に基づく計算でありながら、比較的によく計算できていることがわかるだろう。時間の経過とともにqとpの厳密解と１次精度のオイラー法を用いた数値解との食い違いが大きくなるのは、近似計算を繰り返すたびに誤差が伝播し、蓄積するため。このため、オイラー法の誤差は雪だるま式に増加し、計算はやがて発散する。
次に、オイラー法を用い、Δt=0.1とし、t=10まで数値計算結果を(q,p)のグラフとして表したものを以下に示す。
解いた問題は、力学的エネルギー、つまり、ハミルトニアンHが一定であり、
(q,p)は単位円上に存在しなければならないのだが、時間の経過とともに軌道q²+p²=1からの逸脱が大きくなってゆく。
つまり、オイラー法を用いた近似計算法はエネルギー保存則を満たす解法ではなないということだ。
オイラー法のこの欠点を改良したものにシンプレクティック（Sympretic）法と呼ばれるものがある。
オイラー法で用いる計算式（４）を次のように改良したもの。
この僅かな改良によってエネルギー保存則からの逸脱がかなり改善され、オイラー法の宿命的な欠点である時間の経過に伴う誤差の伝播、蓄積が解消される。
以下に、表計算ソフトを用い、シンプレティック法を用いて解いた計算結果を示す。
比較参照のために、シンプレクティク法とオイラー法を用いて解いたqと厳密解をのグラフを以下に示す。
さらに、シンプレティック法による数値解を(q,p)のグラフで表したものを示す。
オイラー法とは異なり、シンプレクティック法を用いた数値解はかなり良くq²+p²=1の円周上に存在していることがわかるだろう。
シンプレクティク法とオイラー法を用いて解いたq,pを元に計算したハミルトニアンHと厳密な値1/2との相対誤差
の時刻による推移をグラフに示す。
このグラフが示すように、１次精度のシンプレクティク法は完全にエネルギー保存則を満たすわけでないが、振動をしつつも常にある誤差の範囲内に収まり、オイラー法のように指数関数的に誤差が増大することはない。
参考までに２次精度をのシンプレクティク方による計算結果も示してある。
数値計算ではよくあることなのだけれど、（４）式を少し変え（５）式にし計算するだけで計算結果が劇的に変化するのだから、数値計算は奥が深いよね。
そう思いませんか？
そして、少しでもこうした数値計算に興味を持ってもらえれば幸せです。
なお、これらの議論をより正確にするためには、ニュートン力学の進化形である解析力学などの知識が必要になるので、感覚的に理解できれば十分ではないでしょうか。
参考までに２次精度のシンプレクティック法で解いた計算結果は以下に示す。
