グリーン関数による1次元の運動方程式の解法

グリーン関数を用いて1次元のニュートンの運動方程式を解きます。

  1. グリーン関数
  2. 重力下の質点
  3. 重力下のばね
    1. 1つ目の選び方
    2. 2つ目の選び方
  4. 注釈
  5. Thanks(フォントについて)

グリーン関数は非斉次微分方程式の特解を求める為に良く利用されます。
やっていることは、
微分方程式の満たす関数を、あたかも解いた風に書くこと
です。
なので実際、グリーン関数を使う、とは単なる式変形(微分方程式→積分方程式)をしているにすぎません。グリーン関数を本当に求め終わった時、はじめて解けた、と言えるのです。

物理現象では、グリーン関数が表現するのは「無限小時間に力が働く衝撃を受けた結果」を表現します。
どんな力も「複数の衝撃を積算したもの」と捉えることが出来ます。
衝撃に対する応答を考えてから、任意の力を衝撃の足し合わせと考え、問題を解きます。
すると、各々の力に対してその都度微分方程式を解いていくよりも
グリーン関数を求めた後の結果を足し合わせて任意の力に適応する、という考えの方が便利そうだと想像できます。

グリーン関数


微分方程式

をグリーン関数を用いて解くことを考えます。
ここで、左辺の\(\hat{D}\)は演算子を表し、例えば

のような形です。
グリーン関数で解くとは、解を(一般解)+(特解)で書くときの特解を与える手段です。
特解なので、何かしらの境界条件の下でグリーン関数を求める訳です。
物理の問題では、”因果律”と呼ばれる、結果は原因の後の時間に起こる、ということを仮定した境界条件を採用します。

全ての演算子を左辺に移してしまいましたが、右辺に時間に関する微分演算子が入っていても問題ありません。
式(1.1)の形式的な解は、グリーン関数\(G(t,t’)\)を用いて

と、式(1.3a),(1.3b)を満たす解を用いて書くことが出来ます。
ここで、\(\delta(x)\)はディラックのデルタ関数です。
実際に解(1.3)は式(1.1)を満たすことを確かめます。
式(1.3)に左から演算子\(\hat{D}\)を作用させます。
ここで、\(\hat{D}\)は\(t\)に関する演算子なので、積分と交換できます。よって

と変形でき、まさに解(1.3)は式(1.1)を満たす解として表されます。

重力下の質点


さて、質量\(m\)を持つ質点が、重力加速度\(g\)の重力下で運動する問題を考えます。運動方程式は

です。質量\(m\)で割って

を得ます。
もちろん、この問題の解は
\(
\displaystyle x(t)=-\frac{1}{2}gt^2+at+b
\)

です。この問題はグリーン関数を用いて解く必要は全くありませんが、解析解が簡単に求められる練習問題として良いでしょう。
グリーン関数は単なる形式ですが、一般的な解法であり、非常に強力です。
しかし、こんな簡単な問題も解けないのであればお話になりません。

さて、式(1.6)の形式的な解は

と書くことが出来ます。ここで、\(x_0(t), G(t,t’)\)は以下の微分方程式を満たします。

\(x_0(t)\)については簡単に解けてしまい、

です。ここで、\(\alpha, \beta\)は\(x(t)\)の初期条件によって決まる定数です。

続いてグリーン関数の微分方程式を解きましょう。
方程式の右辺を見ますと、\(t= t’\)だけで変なことが起こっています。
ですが、\(t\ne t’\)であれば右辺は完全にゼロです。
すなわち、\(G(t,t’)\)は\(t\ne t’\)で、微分方程式

を満たしています。これは簡単に解けてしまい、

です。ここで、\(A_0, A_1, B_0, B_1\)を適当な定数としておきました。
これらの係数間には、とある関係式が成り立っています。それを調べましょう。

デルタ関数を発散する点を含んで積分すると定数1になることを利用します。
式(1.18)の下の式の両辺を時刻\(t=t”\)周りで積分します。
すなわち、微小な時間\(\varepsilon\)を用いて、

を計算します。
もしも、たまたま\(t”=t’\)であれば、積分を実行して

を得ます。これは、時刻\(t=t’\)周りでは、グリーン関数の導関数は、右から極限と左から極限は一致しないことを言っており、不連続になると言っています。
\(t”\ne t’\)であれば、この極限は一致し、グリーン関数の導関数は連続であるということを言っています。

式(1.11)を代入すれば、係数間の関係式

を得ます。
式(1.7)のグリーン関数を含む項は特解ですので、不明な定数は存在しません。
今、4つの未知の定数\(A_0, A_1, B_0, B_1\)に対して式の数は2つなので、全ての係数を決めるには2つ条件式が足りません。
そこで、境界条件を加えます。この条件は式(1.7)の解の形を見て決められます。
運動方程式を解こうとする場合、通常は物理的な要請から、ある時刻の運動はその時刻より過去の運動のみから決定される、という因果律が課されます(※1)。
さて、因果律を仮定すると、事が起こる時刻\(t=t’\)より前には何も起こらない、すなわち\(G(t,t’)=0\)を意味しており、\(A_0=A_1=0\)となります。よって、

を得ます。よって、グリーン関数

を得ます。ここで、\(\theta(x)\)はヘヴィサイド関数を表します。
グリーン関数(t’=5に固定)を書いてみるとこんな感じ。

黒線がグリーン関数、赤線がその導関数です。導関数の右からと左からの極限の差がちょうどデルタ関数の積分、すなわち1になっていることが分かります。

これでグリーン関数が得られました。式(1.16)を使って式(1.7)の積分を計算すると、

を得ます。ここで、\(c\)は任意の値ではないことに注意してください。定まっている値です。

式(1.9)と式(1.17)の結果を代入すれば、運動方程式の解

を得ます。ここで、\(\beta+c\)を改めて\(\beta\)と書きました。
この結果は本当の解に一致します。

重力化のばね


続いて、鉛直方向に置かれたばねと重力の問題を考えましょう。
ばね定数を\(k\)と置けば、運動方程式は

です。質量\(m\)で割り、\(\omega=\sqrt{k/m}\)と置くと、

を得ます。
解は平衡点を中心とした振動運動です。平衡点\(x_0\)は力のつり合い
\(
-mg=kx_0~~~\to~~~x_0=-\frac{g}{\omega^2}
\)
から導くことが出来るので、
\(
\displaystyle x(t)=a\sin(\omega t +\varphi) -\frac{g}{\omega^2}
\)

です。これをグリーン関数を用いて導くことが出来れば目的達成です。

さて、グリーン関数の考え方として少なくとも2通りの方法が考えられます。
どちらの方法でも同じ解を与えるので、どちらをえらんでも構いません。
数学的に解きやすい方法を選ぶのが賢い方法です。

実際に試してみましょう。

1つ目の選び方

微分方程式を

として考えます。
解を

と考えるのです。\(x_0(t), G(t,t’)\)はそれぞれ

の微分方程式を満たします。式(1.9), 式(1.16)の結果を用いれば、容易に求められて、

です。式(1.22)の積分部分は

なので、解

を得ます。
さて、あらわな形ではないものの、解が得られました。
式(1.27)を満たす\(x(t)\)を探さなければならないのですが、
正直なところ解を見つけられる気がしません。これは自己無撞着の方程式であり、
例えば散乱問題等のように、グリーン関数を得るのが難しい時などでこの形式にとどめておく方法がとられます。
あらわに求める事が出来ずとも、他の式へ代入が出来たり、摂動展開なんかに用いられるので、微分方程式の形式よりその後の計算が扱いやすいのです。
ただ、今の目的からそれてしまうので、ここではこれ以上は話しません。

2つ目の選び方

微分方程式を

のように選びます。すると、解は

であり、\(x_0(t), G(t,t’)\)はそれぞれ

の微分方程式を満たします。
\(x_0(t)\)はただのばねの問題なので、

という解を得ます。一方、グリーン関数も\(t=t’\)を境に場合分けすると

を得ます。時刻\(t=t”\)周りの微小時間で積分すれば、

を得て、たまたま\(t”=t’\)であれば、

の関係式を得ます。ここで、グリーン関数自体はいかなる時刻でも右からと左から極限が一致するという事実をもちいました。
因果律を境界条件として使用すれば、係数は

であるので、グリーン関数は

と求められます。ここで、\(n\)は三角関数の周期性からくる任意の整数です。具体的なグリーン関数(t’=5に固定)の格好はこんな感じ。

黒線がグリーン関数、赤線がその導関数です。導関数の右からと左からの極限の差がちょうどデルタ関数の積分、すなわち1になっていることが分かります。

式(1.29)の積分を計算すると

と求められます。三角関数を時刻無限に飛ばした時の値ですが、拡張リーマン・ルベーグの補題から得られる事実を用います(※2)。これによると、

という事実が得られるので、式(1.37)の最後の項はゼロです。よって、

を得ます。
したがって、解

を得ます。これは求めたかった平衡点を中心に振動する運動です。

さて、1つ目の選び方で得た解(1.27)に解(1.42)を代入してみましょう。解(1.27)が本当にばねの問題なら、式は満たされるはずです。
式(1.27)の積分は

と計算できます。ここで、\(d\)を任意の定数ではない定数と置きました。
式(1.27)に代入すると、(RHS:right-hand side, 右辺)

を得ますが、今、\(\alpha’, \beta’\)は自由に決めることが出来ます。
もしも\(\alpha’=0,~~\beta’+d=0\)となるように定数を決めれば、解(1.27)は式(1.42)で表される解と同じで、微分方程式(1.20)の解となっていることが分かります。

注釈

※1) 因果律によるグリーン関数の境界条件は、妥当っぽく思われますが、人間が勝手にそうだろう、として経験的に与えている条件にすぎません。導くことは出来ない境界条件です。数学的には特解が何でも1つ与えられればいいので、時刻が未来であろうが特定の時間であろうが何でも構わないのです。
 今回の場合は解きたい微分方程式の右辺が定数(重力のみ)という非常に特別な場合ですが、一般的にはグリーン関数を用いた特解の被積分関数には、左辺の解である\(x(t)\)が含まれます。
この時、因果律を仮定しないと、『ある時刻の運動を計算するために、その時刻より未来の運動を知らなければならない』という理解しがたい状況が生まれます。言い換えれば、『現時刻の運動は、未来の振る舞いが決定されているので、どんな力を新たに加えようとも決して運動を変えることは出来ない』となってしまい、非物理的に感じてしまいます。
 もう一度言いますが、これは人間が「時間は必ず過去から未来に進む」という経験的に得られた事実を使用しているのであって、この根拠以外に因果律を採用する理由は存在しません。数学的に解こうとする際にはどんな条件でも構わないので、因果律を境界条件としても特解として正しい答えが得られるので良いのです。
 位置に対する微分方程式である場合は、例えば位置が無限大でゼロである条件だったり、ある点で位置と微分が決まっていたり、漸近で特定の形を持っていなければならない、という条件が当てはめられます。

(※2)これいいんですかね、使って。あんまり考えてはいません。もしも時刻が負の∞ではなかったらと思うと不安で仕方ありません。都合がよかったので使いましたが…ダメな気がしないでもないです。

Thanks

アイキャッチ画像のフォントはmagicringを使用しています。
配布ページ
http://inatsuka.com/extra/magicring/
配布元のサークル様
森の中の猫の小屋 http://inatsuka.com/

陰的ルンゲ=クッタ法

もし利用されるのであれば、実際のプログラムは
陰的ルンゲ=クッタ法の高速化
にして下さい。このページのプログラムは上のページのプログラムに比べて100倍くらい遅いです。


陰的ルンゲ=クッタ法とはどのように計算するのか? を主眼において説明します。

陽的ルンゲ=クッタ法は
ルンゲ=クッタ法の説明と刻み幅制御
4次ルンゲ=クッタ法
をご覧ください。

  1. 陽的解法と陰的解法
  2. 一階の微分方程式
  3. 複数の微分方程式
  4. Runge-Kutta-Gauss-Legendre法のButcherテーブル
  5. プログラム
  6. 刻み幅制御
  7. 参考文献

陽的解法と陰的解法


陽的ルンゲ=クッタ法のButcherテーブルは以下の通り

テーブルの右上と対角要素の値がゼロです。

しかし、陰的ルンゲ=クッタ法のButcherテーブルは

であり、すべての要素に値が入っています。

陽的解法は\(k_i\)を求めるためには\(k_{j\lt i}\)の要素が分かっていれば順繰りに計算できますが、
陰的解法は\(k_i\)を求めるために\(k_j\)の全ての要素が分かっていなければなりません。すなわち、陰的解法では1ステップ進めるために自己無撞着の連立方程式を解かなくてはいけないことを意味しています。

陰的解法の利点は、硬い方程式をスムーズに解くことが出来る点や、シンプレクティック性を持つ公式が存在する点です。

では、陰的公式の数値解法について述べて行きます。
急に任意の数の連立微分方程式を考えるのは複雑なので、方程式の数が一つだけの問題をまず考えて、その後一般化して考えます。

一階の微分方程式


問題

を初期条件\(y(x_0)=y_0\)の下で数値的に解き、\(y(x_0+h)\)の値を数値的に得ることを考えます。
\(s\)段の陰的ルンゲ=クッタ法を用いることを考えると、解は(1.2b)の形で書くことが出来ます。

ここで、係数\(k_i,~~(i=1,2,\cdots,s)\)は

を満たします。式(1.3)を式変形して、関数のゼロ点を探す問題に焼きなおします。すなわち、

の通り、関数\(g_i\)のゼロ点を探すのです。

関数\(g_i\)のゼロ点を与える式(1.4)の解を\(\mathbf{k}^{(\infty)}\)と書き、
解の近傍に存在する\(\mathbf{k}^{(0)}\)が刻み幅\(h\)のオーダーに収まる範囲を考えます。すなわち、

の範囲を考えます。\(\mathbf{k}^{(0)}\)を初期値とすれば、\(\mathbf{k}^{(0)}\)は\(h\)のリーディングオーダーの範囲で

と与えられます。ここで、関数\(y(x)\)もオーダー\(y(x_0+h)-y(x_0)\in O(h^1)\)の範囲にあることを仮定しています。

この想定の下ではNewton-Raphson法を用いることが出来て、

を繰り返して解が収束するまで計算を繰り返せば良く、初期値は\(\mathbf{k}^{(0)}\)です。
ここで略記した表記は

を意味します。また、行列Jは

で表される行列です。

\(J\)の行列要素を求める事はさておいて、Newton-Raphson法を用いるためには

解いて\(\mathbf{z}\)を求めなければなりません。\(J\)の逆行列を求める事は数値的にメリットは無いので、行列のLU分解からこの問題を連立方程式を解く問題に変えます。
式(1.12)のように変形します。

行列が下三角行列\(L\)と上三角行列\(L\)の積で書けるとします。すなわち、行列\(J\)を

の通りに分解するのです。連立方程式

を解くために、LU分解より、

です。\(U\mathbf{z}=\mathbf{w}\)と置けば、

をまず解いて、\(\mathbf{w}\)を得ます。その後、

を解くことで\(\mathbf{z}\)を得ます。

行列要素の計算

続いて、行列要素を計算することを考えます。

もう少し計算進めると、行列要素は

となります。ここで、\(\delta_{n,m}\)はKroneckerのデルタです。

簡易Newton法

行列\(J\)の要素をNewton法の繰り返し度に毎回計算するのは大変です。
なので、\(O(h^1)\)の範囲で関数の符号などが大きく変わらないと想定します。

となります。
すると、行列要素\(J\)は

と表されます。すなわち、Newton法の繰り返し毎に行列\(J\)を計算する必要はなく、

ということです。

複数の微分方程式


続いて以下のN本の微分方程式を解くことを考えます。

関数\(g\)のゼロ点を探す問題に焼きなおすと、

となります。
Newton法で繰り返すことを考えると、

です。それぞれ

を満たします。行列\(J\)は、

であり、\(\mathbf{k}\)の初期値は

です。行列要素\(J\)は

です。簡易Newton法では

として得られます。

ルンゲ=クッタ=ガウス=ルジャンドル法のButcherテーブル


ガウス=ルジャンドル求積法に基づくルンゲ=クッタ法のButcherテーブルは以下のものです。

2段4次

3段6次

このガウス=ルジャンドル求積法に基づく陰的ルンゲ=クッタ法は凄まじく素晴らしい特性を持ちます。

  • A安定である[4]
  • シンプレクティックである[3,4]
  • 対称である[4]
  • 硬い方程式に向く[1,3,4]
  • 刻み幅を途中で変更してもシンプレクティック性が保たれる[3]

至れり尽くせりの性質を持ちますね。とりあえずこれを使っておけば問題ないということです。

プログラム


Fortran90による、N本連立方程式を陰的ルンゲ=クッタ法で解くプログラムです。
簡易Newton法を用いています。

lapackを用いるので、mklかlapackにリンクしてコンパイルしてください。

3段6次 陰的Runge-Kutta

2段4次 陰的Runge-Kutta

さて、実際に解析解と刻み幅を変えて精度の比較をしましょう。
正しく実装できていれば刻み幅を1桁減らすと解析解との一致桁数が4桁、6桁合うようになるはずです。
横軸を刻み幅、縦軸を解析解との相対誤差にとったものが以下のグラフです。
問題はばねの問題で、\(x=100\)で比較しました。

確かに刻み幅を1桁減らすと解析解との一致桁数が4桁、6桁一致するようになっています。
上記計算は倍精度で行っています。なので、3段6次公式は倍精度演算にとってはオーバースペック気味です。
なぜなら、刻み幅が0.05辺りで桁あふれのような振る舞いが見えてしまっているからです。4倍精度演算になればちょうどいいくらいでしょう。

ここ以降で私が書いている部分はあまり良くない説明とプログラムです。
もし刻み幅制御のプログラムが欲しければ、
陰的ルンゲ=クッタ法の高速化
を参照してください。

刻み幅制御


3段6次の公式には別の次数が埋め込まれていますが、それは2次のオーダーです。ちょっと6次のオーダーに対して2次のオーダーを用いるのはもったいない気がします。

ここでは、良くない方法ですが、別の次数を得るために別の陰的公式を用いることを考えます。
別の陰的公式として、LobattoⅢCと呼ばれる3段4次の公式を採用しました。
端点の値が評価されるという性質があるからです。この性質が無いと不連続点をうまく検知できませんでした。
刻み幅制御のやり方は陽的Runge-Kutta-Fehlbergの方法と同じにしましたので、ここでは書きません。
詳しくはルンゲクッタ法の説明と刻み幅制御 -シキノートを参照してください。
陰的公式の正しい刻み幅制御の方法は文献[1]に載っていますので、本格的に知りたい方はそちらを参照してください。

刻み幅を制御した後、実際に返される値は3段6次の結果です。
懸念事項があります。3段6次のGauss-Legendre求積法に基づく公式はシンプレクティック積分ですが、
刻み幅を変えてしまうとシンプレクティック性が失われているのではないか?という事項です。
私が見聞きしている限り、シンプレクティック性を保つためには等間隔刻み幅をもちいなければならないはずです。

@adhara_mathphysさんによるコメントを戴きまして、Gauss-Legendre求積法に基づいた陰的Runge-Kutta公式は刻み幅の変化に依らずシンプレクティック性が保たれるようです。


またシンプレクティック性を保つ積分法は、全力学的エネルギーが正準変数の二次形式で書けるならどんなに時間発展してもエネルギーが保存されるという性質を持ちます。ばねの問題であればこの性質を満たし全力学的エネルギーは保存しますが、ケプラー問題は満たさないため、シンプレクティック性は保ちますが、全力学的エネルギーは保存しません。

3段6次と3段4次を利用した刻み幅制御陰的Runge-Kutta法

大規模計算にはメモリ等を考えていないので向いていません。ご注意ください。

3段6次:Gauss–Legendre methods(実際に返される値)
3段4次:LobattoⅢC (誤差を見積もるためだけに利用)

下のプログラムは
\(
\frac{dy}{dx}=y \cos(x),~~y(0)=1
\)

を解くものです。解析解は
\(
y(x)=e^{sin(x)}
\)

となります。


用途に合わせたその他プログラム

その他、必要と思われるプログラムをまとめて公開します。
4つに分けました。

プログラムの名づけ方

倍精度 d
4倍精度 q

LAPACK使用 l
GECP使用 g

openMP使用 m
不使用 n

自動間隔出力 a
等間隔出力 e

例えば、倍精度でLAPACKを用いて、openMPを使用して等間隔出力する場合のプログラム名は
main_dlme.f90という名前にしてあります。

デフォルトではばねの問題
\(
\begin{align}
\frac{d^2}{dt^2}y(t) = -y(t)
\end{align}
\)

を解く問題になっています。
main_dlma.f90
main_dlme.f90
main_dgma.f90
main_dgme.f90
main_qgma.f90
main_qgme.f90

GECPとは

LU分解と連立方程式の解法を、LAPACKを使わないで、
九州大学の渡部 善隆様が公開なさっているGECP(Gaussian Elimination with Complete Pivoting, 一般実行列に対する連立1次方程式の数値解を完全ピボット選択付き Gauss の消去法によって求める Fortran サブルーチン)
に置き換えたものです
http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/GECP/index.html
再配布可能なので、上記プログラムで必要な物には組み込んであります。

参考文献


[1]E. Hairer and G. Wanner, ‘Solving Ordinary Differential Equations II’ Springer, 1996

[2]ルンゲ=クッタ法のリスト

[3]J. M. Sanz-Serna, ‘Runge-Kutta schemes for Hamiltonian systems’, https://www.researchgate.net/publication/200175547_Runge-Kutta_schemes_for_Hamiltonian_systems

[4]遠藤 理平, 数値計算法概論 常微分方程式の解法, http://www.natural-science.or.jp/article/20120323225814.php

その他参照したページ


遠藤 理平, 数値計算法概論 常微分方程式の解法 http://www.natural-science.or.jp/article/20120323225814.php

田中正次, 三村和正, 山下茂, 陰的Runge-Kutta法の特性について(常微分方程式の数値解法) https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/100230/1/0643-4.pdf

渡部 善隆, 一般実行列に対する連立1次方程式の数値解を完全ピボット選択付き Gauss の消去法http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/GECP/index.html