常微分方程式の数値解法 3

今回の講義も常微分方程式の数値解法のベクトル編です．

例題として，\(d=4\) の場合の次の微分方程式（Keplerの2体問題）を考えましょう： \begin{equation} \frac{\mathrm d}{\mathrm dt} q_1 = p_1, \quad \frac{\mathrm d}{\mathrm dt} q_2 = p_2, \quad \frac{\mathrm d}{\mathrm dt} p_1 = -\frac{q_1}{(q_1^2 + q_2)^{3/2}}, \quad \frac{\mathrm d}{\mathrm dt} p_2 = -\frac{q_2}{(q_1^2 + q_2)^{3/2}}. \end{equation} これまでの表記との対応は \(y = [q_1,q_2,p_1,p_2]^\top\) です． よく知られているように，\(q=(q_1,q_2)\) は，一体を原点に固定したときのもう一体の位置ベクトルを，\(p = (p_1,p_2)\) は速度ベクトルを表します． また， \begin{equation} H(q,p) = \frac{p_1^2+p_2^2}{2} - \frac{1}{\sqrt{q_1^2+q_2^2}} \end{equation} とおくと，これは保存量になります（エネルギーやハミルトニアンと呼ばれます）． すなわち，微分方程式の解に沿って \begin{equation} \frac{\mathrm d}{\mathrm dt} H(q(t),p(t)) = 0 \end{equation} が成り立ちます． また， \begin{equation} L(q,p) = q_1 p_2 - q_2p_1 \end{equation} という保存量も持ちます．この量は角運動量と呼ばれ，これが保存量であることは，Keplerの第2法則に他なりません．

今回の講義では，Kepler方程式をいくつかの解法で数値計算した様子を比較してみます．

Symplectic Euler法

前回の講義で学んだsymplectic Euler法は次のようになります． Kepler問題を \begin{equation} \frac{\mathrm d}{\mathrm dt} q = p, \quad \frac{\mathrm d}{\mathrm dt}p = g(q) \end{equation} と表すことにすると，symplectic Euler法の更新式は次のようになります： \begin{align} q_{n+1} &= q_{n} + h p_{n}, \cr p_{n+1} &= p_{n} + h g(q_{n+1}). \end{align} Symplectic Euler法は1次精度の解法です．

Störmer-Verlet法

Symplectic Euler法を2次精度に拡張したStörmer-Verlet法の更新式は次のようになります： \begin{align} q_{n+1/2} &= q_{n} + \frac{h}{2} p_{n}, \cr p_{n+1} &= p_{n} + h g(q_{n+1/2}), \cr q_{n+1} &= q_{n+1/2} + \frac{h}{2} p_{n+1}. \end{align} Störmer-Verlet法は2次精度の解法です．

Störmer-Verlet法は，天文学の文脈で，1907年にStörmerにより提案された解法です． 今日に至るまで，天文学の文脈で広く利用されています（より正確に言えば，Störmer-Verlet法を高精度に拡張した解法がよく利用されています）． 一方で，分子動力学計算の文脈で，1967年にVerletが（Störmerとは独立に）本質的に同じ解法を提案したという経緯があります． そのため，分野によっては，Störmer法といったりVerlet法といったりすることもあります．

なお，この解法を歴史上はじめて用いたのはStörmerではありません． さらに200年以上遡り，NewtonがPrincipiaの中で，ケプラーの第2法則を幾何的に理解するために用いたのが初出です． そのため，数値解析学者の中には，Newton-Störmer-Verlet法と呼ぶ人もいます．

陰的な解法

今日の講義では，陰的な解法も扱います． ここでは簡単のため，二つの解法をプログラミングしてみます．

• 陰的Euler法
  \begin{equation} y_{n+1} = y_n + h f(y_{n+1}) \end{equation}
• 陰的中点則
  \begin{equation} y_{n+1} = y_n + h f\Big(\frac{y_{n+1}+y_n}{2}\Big) \end{equation}

どちらも，時間発展を計算するためには，連立非線形方程式を解く必要があります． 実用上は，不動点反復法・簡易Newton法・Newton法などがよく利用されますが，ここでは，パッケージの利用の練習も兼ねて，NLsolveパッケージを利用して数値計算することにします．

今日の課題

今日の課題は，Störmer-Verlet法と陰的中点則のプログラムを完成させることです． それぞれ，Symplectic Euler法と陰的Euler法のプログラミングが参考になります． これらを完成させたうえで，数値解に対して，エネルギーと角運動量の振る舞いを観察していきます．

サンプルファイルはこちら