sikino のすべての投稿

束縛条件下の運動 – 自由度がうまく落とせない運動

このページは
束縛条件下の運動 – ホロノミックな束縛と非保存力
の続きです。

拘束がある場合の2次元の運動


ラグランジュの方程式を解く際に適切な座標変換が見付からず自由度が落とせない場合を考えます。その時の運動方程式は

の形で止まります。式(75)において、未知の関数は\(x(t), y(t), \lambda(t)\)の3つであり、条件式は3つなので解くことは出来るはずです。
方針は、式(75)から\(\lambda(t)\)を消去すれば良いのです。

束縛条件\(f(x,y,t)=0\)が成立する点\(x=x_0, y=y_0\)周りで、微小時間\(\Delta t\)の変化を考えます。
点\((x_0, y_0)\)まわりでテーラー展開すれば、

を得ます。よって

が満たされる変化でなければなりません。すなわち、\(x,y\)の時間変化は式(77)をいつも満たしていなければならないのです。両辺を\(\Delta t\)で割って極限をとれば、

が成立することが分かります。式(75)を書きかえれば、

を得ます。式(79c)もまた時間変化しても成立し続けなければなりません。よって、左辺の時間微分を取れば、

を得ます。整理して

を得ます。ここで関係式

を用いました。式(79a), (79b)から未定乗数\(\lambda(t)\)を消去すれば

を得ます。よって、式(81)と式(84)から連立方程式

を立てることが出来ます。これをあらわに解けば、運動方程式

を得ます。ここで、式(86), (87)の右辺の\(-\frac{\partial V}{\partial x}, -\frac{\partial V}{\partial y}\)は\(x, y\)方向の保存力であり、右辺の残りの項は束縛条件によって決まる束縛力であることが分かります。

拘束がある場合の3次元の運動


スカラー関数\(f\)で表される拘束がある条件のもとで, 運動方程式

またはベクトル表記で

について考えます。ここで表記を略すために

の略記を用いました(\(y, z\)方向についても同様)。
また、\(\nabla\)はナブラ演算子で、

です。2次元の場合と同様に、式(88d)の時間微分から

が導けます。ここで、\(\hat{H}f\)は関数\(f\)のヘッセ行列であり、

を表します。式(92)は、式(81)を3次元に拡張した形になっています。

式(88a),(88b),(88c)から未定乗数を消去すると従属な3つの関係式

を得ることが出来ます。この3式で独立なのは2つです。
独立な式として式(94a)と式(94b)を選ぶことにし、式(92)と共に書けば、連立方程式

を得ます。この式をあらわに解けば、運動方程式

を得ることが出来ます。ここで、\(n_x, n_y, n_z\)は法線ベクトル\(\mathbf{n}\)の\(x,y,z\)成分です。偏微分ではないことに注意してください。

ベクトル表記であれば

と書くことが出来ます。

単振り子の例


ここでは単振り子を例に挙げ、式(86)に代入した時本当に振り子を表す方程式になっているのか、また束縛条件の選び方に依らず、同じ運動を記述しているのかを例に挙げます。
保存力は重力で

を考えます。

考える束縛条件は
\(
f(x,y,t)=x^2+y^2-l^2=0
\)


\(
f(x,y,t)=\sqrt{x^2+y^2}-l=0
\)

です。

束縛条件1


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

束縛条件2


束縛条件として

を考えます。この束縛条件に対する偏微分を計算して、

を得ます。これを式(86), (87)に代入すれば、

であり、整理すれば

を得ます。

考察


同じ束縛条件なので、式(105)と(109)の示す運動は同じになるはずです。
それを示すには、2式の差異である

を示せればよい、ということです。そのために拘束力に対して座標変換

を考えます。右辺と左辺をそれぞれ計算すれば、

となり、同じ運動を記述していることが分かります。

参考文献


6. 振り子 -物理学の見つけ方
9. 自由な座標 -物理学の見つけ方

束縛条件下の運動 – ホロノミックな束縛と非保存力

このページは
束縛条件下の運動 – ホロノミックな束縛 1
の続きです。

単振り子



ポテンシャルが存在し、ホロノミックな束縛が課されている代表的な問題である、単振り子を考えましょう。
単振り子は、束縛条件下の運動 – ホロノミックな束縛 1で考えた円に束縛されている問題に対してポテンシャルが加わった問題です。

座標系は以下のような極座標系を選びます。

通常の極座標系とは異なります。\(\theta=0\)で鉛直下方向にとる座標系にしたかったので、通常の極座標系とは変えています。

ホロノミックな束縛条件\(f(x,y,t)=0\)を含めたラグランジアンは、デカルト座標系から極座標系に変えて

と書くことが出来ます。注記しますが、ラグランジュの未定乗数は時間に依存していて構いません。
ラグランジアンの偏微分をそれぞれ計算すれば、

ですので、運動方程式と束縛条件

が得られます。連立微分方程式(31)の3本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の3つなので、方程式を立てて\(r(t),\theta(t)\)を求めることが出来ます。

では、伸び縮みしない振り子を考えましょう。束縛条件は

をして表現することが出来ます。束縛条件の一般化座標による偏微分は

ですので、運動方程式

を得ます。式(34)の第2式を変形すると

ですので、これはまさに振り子を表す運動方程式です。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
 \(m=3, l=1, g=9.8\)
\(\theta(0)=\frac{4\pi}{5}, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

支点が動く振り子



次に、束縛条件が2つある場合を考えましょう。その例題として、
振り子の支点が摩擦の無いレールの上に載っていることを考えます。
ラグランジアンは

です。2つのホロノミックな束縛条件

を考えます。すると、ホロノミックな束縛条件を含めたラグランジアンは、

と書くことが出来ます。
これから伸び縮みしない振り子を考えることを見越して、座標系

を考えます。この新たな座標系でのラグランジアンは、

と書けます。それぞれの偏微分である



を用いると、運動方程式

を得ます。
連立微分方程式(46)の4本の方程式と2本の束縛条件に対して、未知の関数は\(x(t), y(t), r(t), \theta(t), \lambda_0(t), \lambda_1(t)\)の6つなので、方程式を立ててそれぞれの関数を求めることが出来ます。

では、支点が直線状のレールの上に置かれている状況を考えましょう。
支点に関する束縛条件\(f_0\)と, 振り子の紐の長さが変わらないという束縛条件\(f_1\)を考えると

と書くことが出来ます。なので、

とそれぞれ計算できるので、運動方程式

を得ます。求めたい未知の関数は\(x(t), \theta(t)\)です。
式(51a)と(51b)から、\(x(t), \theta(t)\)に関する運動方程式

を得ます。変形すると

ですので、あらわに求めれば、

を得ます。
特に支点の質量が振り子の質量に比べて非常に重い場合、すなわち\(m_0 \gg m_1\)の場合、

になります。第2式はまさに単振り子の方程式となっています。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
 \(m_0=3, m_1=1, l=1, g=9.8\)
\(x(0)=0, x'(0)=0, \theta(0)=\frac{4\pi}{5}, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

支点が放物線上に束縛されている振り子



座標系は1つ前の問題と同じく、

にとることにします。すると運動方程式は式(46)で表されるので、運動方程式の導出過程を省略することが出来ます。
問題設定から、束縛条件は

と書くことが出来ます。なので、それぞれの偏微分は

と書けますので、運動方程式は

です。式(60a)と(60b)から\(\lambda_0(t)\)を消去すると、

を得ます。更に、束縛条件から座標\(y(t)\)の時間変化は\(x(t)\)を用いて

と書くことが出来ます。式(61)と式(60d)から、\(x(t), \theta(t)\)の運動方程式

を得ます。あらわに解くと式が長くなってしまうので、ここで止めておきます。

実際に解きますと、

というような振る舞いになります。

・パラメータ
 \(m_0=3, m_1=1, l=1, g=9.8\)
 \(x(0)=0, x'(0)=5, \theta(0)=0, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

単振り子と空気抵抗


束縛条件と、非保存力である空気抵抗が存在する場合を考えます
座標系は

です。ホロノミックな束縛条件を含めたラグランジアンは、

と書けます。

さて、非保存力である空気抵抗が存在することを考えましょう。
非保存力なのでラグランジアンで書くことは出来ませんので、個別に求める必要があります。
空気抵抗の力\(F_x, F_y\)はデカルト座標系において

で掛ける事を既知とします。今知りたいことは、式(64)で表される一般化座標\(r, \theta\)で式(67)がどのように表されるのか?ということです。
\(r, \theta\)において、一般化座標における空気抵抗\(Q_x, Q_y\)は、

で変換することが出来ます。具体的にそれぞれの力を示せば、

ですので\(Q_x, Q_y\)は、

と表せられます。非保存力がある場合、ラグランジュの方程式の右辺にそれを入れれば良いので、

と表されます。よって運動方程式

を得ます。
長さ\(l\)の伸び縮みしない紐に繋がれた振り子をを考えれば、束縛条件は

ですので、運動方程式

を得ます。

実際に解きますと

・パラメータ
 \(m=3, l=1, g=9.8, c_1=0.3, c_2=0.3\)
 \(\theta(0)=0.8\pi, \theta'(0)=0\)
・計算手法
 陰的ルンゲクッタ法

このページはここで終わりです。
続いて、自由度がうまく落とせない、うまい座標系が見付からない場合を考えます。

次のリンク
束縛条件下の運動 – 自由度がうまく落とせない運動

参考文献


田辺行人、品田正樹著『理・工基礎 解析力学』 裳華房

束縛条件下の運動 – ホロノミックな束縛

ホロノミックな束縛を受けている質点の運動を考えます。

ホロノミックな束縛とは、
”束縛条件が一般化座標と時間のみで決まる方程式によって解析的に表現できるもの”
です。空気抵抗のない単振り子などは、半径(一般化座標)が一定という条件のみで決まるので、ホロノミックな束縛に相当します。

ラグランジュの方程式を基本として考えていきます。

スタートはラグランジアン、ゴールは運動方程式を立てる事です。

2次元平面上に束縛されている質点


1つの質量を持った質点が、ある軌道上しか動けない運動に対する運動方程式を導きましょう。
ある軌道上しか動けないという条件は、時刻\(t\)における質点の位置\((x,y)\)が、関数\(f(x,y,t)=0\)を満たす動きしか許されないということで表現することが出来ます。

束縛条件を含めた実効的なラグランジアンは, 未定乗数\(\lambda(t)\)を用いて

と書けます。ラグランジュの方程式

に代入すれば運動方程式

を得ます。

円に束縛されている運動


運動が半径\(l\)の円に束縛されている場合、束縛条件は

として表現されます。もちろん、\(f(x,y,t)=\sqrt{x^2+y^2}-l=0\)としても構いません。運動は束縛条件の形には依らないからです。

式(3)に代入すると、運動方程式

を得ます。方程式(3),(4),(5)の3本の方程式に対して未知の関数は\(x(t),y(t),\lambda(t)\)の3つなので、方程式を立てて\(x(t),y(t)\)を求めることが出来ます。

式(5)の第1式に\(x\)、第2式に\(y\)を掛けて両式を引くと

です。変形すれば

を得ます。どういう解き方でも構わないのですが、半径\(l\)が変わらないことから極座標

を用いることを考えます。すると、式(7)は

と書けます。よって

ここで、\(\omega,\theta_0\)は初期条件から来まる定数です。よって、解

を得ます。ちなみに、\(\lambda(t)\)は

です。

角周波数\(\omega\)で等速回転する筒に入れられた質点が描く軌道



続いて、ホロノミックな束縛条件が時間\(t\)に依存する問題を考えてみましょう。
そこで、角周波数\(\omega\)で等速回転する筒に入れられた質点の軌道を考えます。

極座標表示の方が扱いやすいので、極座標で考えていきます。

ホロノミックな束縛条件を含めた実効的なラグランジアンは、デカルト座標系のラグランジアンから出発して

を得ます。ここで、\(\lambda(t)\)は未定乗数を表します。
計算していけば、運動方程式

を得ます。方程式(18a),(18b),(18c)の3本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の3つなので、方程式を立てて\(x(t),y(t)\)を求めることが出来ます。

角周波数\(\omega\)で等速回転する筒を表現する関数を考えてみましょう。
筒の上の任意の点\((x,y)\)は、

と書くことで表現することが出来ます。よって、束縛条件から

という関係式を得ます。式(21)を運動方程式(18)に代入すれば、

を得ます。式(22a)を解けば、

であるので、解

を得ます。ここで、\(A,B\)は初期条件から決まる定数です。

実際に解いてみると、こんなグラフが得られます。

・パラメータ
\(A=0.01, B=0, w=2\)

このページはここで終わりです。
続いて、ポテンシャルが加わる問題である単振り子と束縛条件を考えていきます。
次のリンク
束縛条件下の運動 – ホロノミックな束縛と非保存力

参考文献


田辺行人、品田正樹著『理・工基礎 解析力学』 裳華房

最速のクイックソート(Fortran)

クイックソートが速いのは分かっていますが、コーディングの仕方によって速度は変わります。
本稿では、ネット上で公開されているどのクイックソートが早いのか調べていきます。

最も早いクイックソートはNUMPACのプログラムでした。

比較プログラム


比較するプログラムはネット上で公開されているクイックソート+αです。
対象は、倍精度実数をソートするプログラム、です。

1. 再帰を用いたクイックソート(f90)

t-nissieの日記: 【電脳】Fortranで書いたクイックソート

2. 再帰を用いないクイックソート(f90)

[Fortran]再帰を使わないquicksortその2 -fortran66の日記

※ただし、私が割と変えていますので、オリジナルそのものではないということを注記しておきます。この変更によって、実行速度はほぼほぼ変わらないことは確認しています。

3. Netlibのクイックソートdsort.f

https://www.netlib.org/slatec/src/dsort.f

4. NUMPACのクイックソート(sortdk.f)

SORTPACK(SORTxK,SORTxy,SRTVxz) (スカラー又はベクトルデータの内部ソーティング)

5. ヒープソート

fortran90によるヒープソートとバブルソート -シキノート

以上の5つのプログラムを比較していきます。
メインプログラムはこちら↓

コンパイルは

gfortran -O3 (ソートのプログラム達) main.f90

で行いました。

結果


結果を示します。
横軸にデータ数、縦軸にソートに掛かった時間を示しました。
”ソートに掛かった時間”とは、同じデータ数でソートを100回繰り返した時の1回当たりの時間です。


それぞれ、
赤線:再帰有りクイックソート
青線:再帰無しクイックソート
緑線:Netlibのクイックソート
紫線:NUMPACのクイックソート
黒線:ヒープソート
です。
最も早いのがNUMPAC, 最も遅いのがヒープソートだと分かりました。

ヒープソートもクイックソートも大体\(O(n\log n)\)ですので、グラフの傾きは両者でほとんど変わりません。再帰有りの方が遅い、という結果は面白いです。
先入観で”再帰は遅い”と考えていたのですが、そんなことはありませんでした。

続いて、同じデータを対数で見てみると以下の通りです。

走査したデータサイズの範囲では変化は有りません。更に大規模になれば、違いが見えてくるかもしれません。

ケプラー問題に対する陽的解法と陰的解法

2体のケプラー問題を数値的に解きます。
ここでは、適切な変換をして求めるのではなく、刻み幅制御されたプログラムで無理やり計算します。

そして、
陽的解法であるルンゲ=クッタ=フェールベルグ法

陰的解法であるルンゲ=クッタ=ガウス=ルジャンドル法
の離心率に対する計算量の違いを調べてみます。

Kepler問題

二次元で二体の問題を考えます。運動方程式は
\(
\begin{align}
\frac{d^2 x}{dt^2}&=-\frac{x}{(x^2+y^2)^{3/2}} \\
\frac{d^2 y}{dt^2}&=-\frac{y}{(x^2+y^2)^{3/2}}
\end{align}
\)
\(t=[0,20],~~x(0)=1-e,~~x'(0)=0,~~y(0)=0,~~y'(0)=\sqrt{\frac{1+e}{1-e}}\)
です。ここで、\(e\)は離心率を表します。

この問題に対する解は良く知られていて、
\(
\displaystyle (x+e)^2+\frac{y^2}{1-e^2}=1
\)

で表され、
\(
0\le e\lt 1
\)

の時、楕円となります。

関数の評価回数の離心率の依存性


楕円の軌道を持つ範囲において、
計算は離心率\(e\)が1に近づくほど難しくなります。
なぜなら、原点付近を通過する際に、質点の導関数の変化が大きくなるからです。

使用したプログラムの説明は
陽的解法はhttp://slpr.sakura.ne.jp/qp/runge-kutta-ex/
陰的解法はhttp://slpr.sakura.ne.jp/qp/irkgl-program/
です。実際のプログラムも置いておきます。

離心率を変化させたときの軌道はこんな感じです。

さて、計算速度の評価ですが、関数が何回評価されたかで比較します。
念の為、陰的解法はLU分解もあるので単純な比較は難しいことを注記しておきます。

図の一番上は関数の評価回数の離心率依存性を表しています。縦軸は評価回数、横軸は\(1-e\)を表します。真ん中の図は一番上と同じですが、軸を対数にとっています。一番下は、\(t=20\)まで計算した時に\(t=0\)の初期エネルギーとの相対誤差\(|(E(t=20)-E(t=0))/E(t=0)|\)を表しています。

特徴的な振る舞いは、離心率に対して陰的解法の評価回数は線形の振る舞い、陽的解法は指数で振る舞っている点です。
これは、質点が原点の近くを通るような問題の際に違いが際立つ事を示しています。
また、エネルギーの保存に関しても陰的解法の方が良いことが分かるでしょう。

一応注意しておきますが、ここでいう陰的解法はルンゲ=クッタ=ガウス=ルジャンドル法の振る舞いです。一般的な陰的解法については話していないことに注意してください。

プログラム


きっかけ

きっかけとしてはtwitterで流れてきまして、やってみよう、と思いました。

陰的ルンゲ=クッタ法の高速化

陰的ルンゲ=クッタ法の基本的な考えは
陰的ルンゲ=クッタ法
にて説明しました。
基本的には
\(
y(x+h)=y(x)+\Delta y
\)

の右辺\(y(x)+\Delta y\)を計算します。
しかし、陰的ルンゲ=クッタ法の方法を数値計算で行おうとすると望ましくない部分が現れます。
それは、

  1. \(y(x)=O(h^0),~~\Delta y=O(h^1)\)を同時に扱わなければならず、桁落ちが激しい
  2. 関数の評価回数が多い
  3. ヤコビアン、LU分解の計算コストが非常に高い

という点です。
簡易ニュートン法を用いる事を前提にしておくと、上の問題は若干解決することが出来ます。

本稿では陰的ルンゲ=クッタ法を発展させ、数値計算的に陰的ルンゲ=クッタ法のアルゴリズムを工夫し、どのように計算量を減らすか?に焦点を当てます。

目次

  1. 計算方法
  2. プログラム
  3. プログラムの評価
  4. 結論
  5. RADAU5の使い方
  6. 参考文献

注意
[2,5]の同著者の新しい論文では、本稿の計算方法([1]に従う方法)ではなく、
もう一段階変換してからニュートン法を利用しています。
正直な所、私自身が追いきれなかっとのと、複素数が入ってくるプログラムでしたので、[1]の方法で止めておきます。
2つの計算方法の収束の早さや精度を比較を[2]で行っていますので、気になる方はそちらをご覧ください。

計算方法


解きたい問題は\(N_{\text{eq}}\)本の連立一次微分方程式

です。これを\(s\)段のルンゲ=クッタ法で解くことを考えます。
座標や添え字は以下のように決めました。

すると次のステップの値は

と求められます。\(z_{np}^{[i]}\)を求める事が問題になります。ここで、\(d_j,~~(j=1,2,\cdots,s)\)は既知で

のように、Butcherテーブルから求められます。
具体的に、3段6次の陰的公式であるGauss-Legendreの場合、

と求められます。
\(z_{np}^{[i]}\)の具体的な形は

で計算されます。\(z_{np}^{[i]}\)をベクトルとして表したいので、

と変換します。
この非線形方程式はニュートン法によって求める事が出来て、以下の通り、\(k\)回の繰り返しで解が収束するまで計算されます。

ここで、解が収束した結果が求めたい\(z_{np}^{[i]}\)となります。すなわち、

です。行列\(J’\)は

と求められます\((n,m=1,2,\cdots,N_{\text{eq}}, p,q=1,2,\cdots,s)\)。ここで、\(\delta_{i,j}\)はクロネッカーのデルタ、\(a_{i,j}\)はButcherテーブルの値、\(J\)はヤコビアンで、

と計算されます。また、ベクトル\(\mathbf{w}_{np}^{(k)}\)は、

と定義します。

初期値の推定


初期値は\(i=0\)の初めのステップでは

として求め、それ以降では過去の結果を多項式補間して求めるとかなりよい精度です[1]。


ニュートン法の停止


ニュートン法は以下の条件が満たされた時、終了します。

まず、\(i=0\)の初めのステップでは必ず\(k\ge 2\)まで行い、
その後、以下の条件が満たされた時、終了します。

ここで、\(\eta_k^{[i]}\)は

で定義され、特に\(\theta_k\)は収束割合(convergion rate)と呼ばれます。
また、\(\kappa=0.1\sim 0.01,~~Ntol=\text{min}(0.03,\sqrt{Rtol})\)です[5]。ここで、\(Rtol\)は相対誤差を意味します。
実際に組んでみますと、理論が違うせいなのか、\(Ntol=\text{min}(0.03,\sqrt{Rtol})\)では十分なほど収束はしませんでした。なので、実際に組む時には\(Ntol=10^{-9}\)にしてしまっています。もしかしたら、\(10^{-9}\)では足らず、もっと必要かもしれません。

\(i\ge 1\)のステップでは\(k=1\)の時は前ステップの\(\eta\)の値を使い、


で判定します。ここで、倍精度演算ならば\(Uround=5\times 10^{-16}\)です。

\(i\ge 1,~~k\ge 2\)では

に従い、計算します。

ニュートン法の繰り返しの\(k\ge 2\)において、どこかで

が満たされる様であれば、刻み幅が大きすぎて収束しないことを意味します。なので、刻み幅を小さくする必要があります。

誤差判定


区間\(i\)の\(n\)番目の方程式の解の誤差は、以下の連立方程式を解いて得られます[1]。

ここで、\(\gamma_0\)はButcherテーブルの行列Aの実固有値で、ガウス=ルジャンドル陰的ルンゲクッタの3段6次であれば、
\(
\gamma_0=0.215314423116112178244733530380696
\)

を得ます。そして、上の方程式を解いた後に、計算結果を棄却するか判定するために、量

を計算します。
もしも\(||\text{err}^{[i]}||\lt 1\)であれば、その刻み幅で計算した結果を採用します。
しかし、\(||\text{err}^{[i]}||\ge 1\)であれば、以下の連立方程式を解きます([1]の”Hump”問題を参照)。

上を計算し、\(||\text{err}^{[i]}||\)を再計算した結果、もしも\(||\text{err}^{[i]}||\lt 1\)であれば、その刻み幅で計算した結果を採用し、そうでなければ刻み幅を次の節に従って変更します。

刻み幅制御


刻み幅制御をするためには, 2つ新しい刻み幅の推定値である\(h_{i+1}^{(1)},h_{i+1}^{(2)}\)を計算します。

その結果、刻み幅が小さくなるか大きくなるかに従って、どちらの刻み幅を採用するか決定します[2]。

プログラム


Fortran90で書いたプログラムはこちら。LU分解と連立方程式を解くため、LAPACKを使います。
それにリンクしてコンパイル、実行をしてください。
モジュールを使用していますが、これは関数の呼び出し回数を計測するためだけにグローバル変数として使っているため、消してもプログラムに何の影響もありません。

追記)
色々計算してみました。その結果、\(tol=10^{-8}\)より小さい値は使わない方が良さそうです。どうもこれ以上の精度にしてしまうと誤差の溜まり具合が増えてしまう感じがします。

ある配列で定義したグリッド上の値

計算結果を、ある配列で定義したグリッド上で欲しい場合、以下のプログラムで行うことが出来ます。
下に載せたプログラムと、上のコードの中にあるサブルーチン(grk, irkgl, dirk6, discrete_h, Jacobian)を一緒にコンパイルしてください。

終点だけの結果が欲しい場合

終点だけの結果が欲しい場合、不要なwork配列などは省略できます。
下に載せたプログラムと、上のコードの中にあるサブルーチン(grk, irkgl, dirk6, discrete_h, Jacobian)を一緒にコンパイルしてください。

収束判定の余地


上のプログラムは収束判定を[1]と同じにしているため、過剰に評価しているパラメータになっているかもしれません。

その余地としては、計算回数を減らすために重要な順に

  1. ニュートン法の収束判定 … Ntol
  2. 刻み幅の安全係数 … fac
  3. 刻み幅の離散化 … discrete_h内のimax
  4. ヤコビアンの更新条件 … if(Newt.le.2.or.theta.lt.1d-3)Jup=0の箇所

です。現状のプログラムでは安全のために過剰評価気味にしています。

4倍精度とGECP


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

また、倍精度、4倍精度に対応しているので、それらのプログラムを使って
陰的ルンゲクッタ法を書いたものが以下のものです。
irkgl_dge.f90
irkgl_qge.f90

倍精度 d
4倍精度 q

プログラムの評価



5つの問題について、評価した結果を上に載せます。評価の良し悪しは、
連立方程式の右辺の関数が何回評価されたかによって決めました。
比較対象は

  1. 陰的解法である本稿のプログラム(自作)
  2. 陽的解法であるルンゲ=クッタ=フェールベルグの4,5次(自作)
  3. 陰的解法であるRADAU IIAに基づくプログラム[3]

です。横軸に要求した精度、縦軸に実際に評価された回数を載せました。
まず、自作の陽的解法、陰的解法を比べますと、硬くない方程式であるEq.1,2,3では
陽的解法の方が10倍近く早いです。
しかし、硬い方程式であるEq.4,5は10倍から1000倍ほど陰的公式の方が早いという結果が得られました。望み通り、陰的解法が動いていることが確認できます。

…さて、専門家が書いたRADAU5ですが、めちゃくちゃ早いです。硬い方程式であるEq.4,5でも、自作したやつの1/10の回数で大体終わっています。
しかも、硬くない方程式であるEq.1,2,3ですら、自作の陽的解法よりも少ない回数で計算を終えています。本当にどういう事なんでしょうね…。
上には上がいるものです。

結論


ちゃんとした陰的解法が欲しいのであれば、自作せず、専門家のプログラムを使いましょう。

RADAU5の使い方


RADAU5のFORTRANプログラムはhttps://www.unige.ch/~hairer/software.htmlにあります。
使い方に手間取ったので、どうやって使うのかメモしておきます。

  1. http://www.unige.ch/~hairer/testset/testset.htmlに移動し、Van der Pol方程式のメインプログラムをダウンロードする。場所は
    ・VDPOL … driver for RADAU5,
    と書かれている所をクリックするとhttp://www.unige.ch/~hairer/testset/stiff/vdpol/driver_radau5.fに飛ぶ。これを保存し、driver_radau5.fという名前で保存する。
  2. https://www.unige.ch/~hairer/software.htmlに移動
  3. Stiff Differential Equations and Differential-Algebraic Problems
    という項目の
    ・RADAU5
    ・DC_DECSOL
    ・DECSOL
    のリンク先のプログラムをダウンロード。それぞれradau5.f, dc_decsol.f, decsol.fという名前で保存する。
  4. 合計4つのプログラムをダウンロードしたら、コンパイルを
    gfortran decsol.f dc_decsol.f radau5.f driver_radau5.f

    でコンパイルし、実行する。

参考文献


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

[2]E.~Hairer and G.~Wanner. Stiff diferential equations solved by Radau methods, J. Comput. Appl. Math., 111:93-111, 1999.

[3]E. Hairer, Fortran and Matlab Codes https://www.unige.ch/~hairer/software.html

[4]10. 常微分方程式 (2)https://www.ktech.biz/jp/archives/1003, K Technologies Sites (2014)

[5]Nicola Guglielmi and E.~Hairer, User’s Guide for code RADAU5 – Version 2.1 (packed in “radar5-v2.1.tar”) http://www.unige.ch/~hairer/software.html, 2005

グリーン関数による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法を用いるためには逆行列

を求めなければなりません。\(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

-1のn乗の計算速度

fortran90で、
\(a=(-1)^n, ~(n=0,1,\cdots, M)\)を計算する早いアルゴリズムは

a=1
if(mod(n,2).eq.1)a=-1

です。


物理、数学をやっていると至る所で\((-1)^n\)を見受けます。
これを数値計算する際に早い計算方法は何でしょうか?

試してみるのは以下のものです。


\(
(-1)^n
\)


\(
(-1)^n
\)
※①を倍精度型で計算


\(
e^{i\pi n}
\)
※\(i\)は虚数単位、\(\pi\)は円周率


\(
\begin{eqnarray}
\left\{
\begin{aligned}
& 1 ,~~(n=0,\mbox{偶数})\\
& -1 ,~~(n=\mbox{奇数})\\
\end{aligned}
\right.
\end{eqnarray}
\)


\(
a=-a
\)で逐次計算

使ったプログラムは以下のものです。

2次元時間依存しないシュレーディンガー方程式の数値解法

2次元の時間依存しないシュレーディンガー方程式を変分原理に基づいて解きます。

対象とする問題は以下の時間依存しないシュレーディンガー方程式です。

ここで、\(H_0\)は数値計算で用いる基底関数のハミルトニアン、\(V(\mathbf{r})\)はその他のポテンシャルです。

この系の固有関数\(\phi(\mathbf{r})\)を

として\(H_0\)の固有関数で展開します。ここで、\(\varphi_{\mathbf{r}}\)は

の固有値問題の解である、固有値\(E’_n\)に属する固有関数です。また、\(\varphi_{\mathbf{r}}\)は

として規格直交化されています。

式(2)を式(1)に代入した後、左から\(\varphi^*_{\mathbf{r}}\)を掛けて全空間で積分すれば

が得られます。ここで、\(V_{m,n}\)を

と置きました。式(7)を別の表現をすれば

と書くことが出来ます。なので、左辺のハミルトニアンの行列形式を対角化すれば、求めたい系の固有値、固有ベクトルが得られるわけです。固有ベクトルが得られたら式(2)に従って波動関数を構成すれば良いのです。

これが一般的な変分原理による解法です。
以降は2次元の問題に限り考えていきます。

2次元の場合


2次元の場合かつ\(H_0\)はxだけ含む部分とyだけ含む部分とで分離することが出来る場合を考えます。

この場合、固有関数は変数分離することが出来るので、2次元の基底関数を1次元系の基底関数の直積をとして考える方法を取ることが出来ます。すなわち、x,y方向の量子数\(n_x,n_y\)を用いて量子数\(n\)と以下のように関係して構成される、と考えるのです。

ここで、\(n\)と\(n_x,n_y\)には1対1の関係があれば順番は関係ないことに注意してください。それぞれの基底関数は

を満たします。もちろん、1次元系の基底関数がそれぞれ正規直交化されていれば二次元の場合でもその関係は保たれ、

が成立します。この基底関数で式(7)を書き換えてやれば

となります。

基底関数の対応


この問題はある長径、短径で指定される楕円内に存在する格子点の個数を数える問題と同じです。

nの順番を記録しておきさえすればどんな順番でも構いません。適当に求めればよいです。
ただし、ここではシュレーディンガー方程式を解いた結果、エネルギーの低い状態から順に欲しいので、\((n_x, n_y)\)の組み合わせでエネルギーが低い順に決定していきます。

今、エネルギーは\(x,y\)方向ともに単調に増加します。
\((n_x,n_y)\)の考えられる組み合わせで、
一番小さいエネルギーは\((1,1)\)なので、まず\(n=1\)は\((1,1)\)に対応させます。
次に小さいエネルギーで考えられるのは\((1,2),(2,1)\)のどちらかです。
ここで重要なのは現在より前に出てきた\((n_x,n_y)\)の組に1を足した組み合わせしか候補にあがりません。
なので、探索するのはこの条件を満たすものだけで充分です。

これをプログラムしたものは、下のtar.gz内に”qnumber”という名前でサブルーチンに入っています。

プログラム


こちらです。
http://slpr.sakura.ne.jp/qp/supplement_data/highprecision_tise2d.tar.gz

解凍すると、中にintegrate.f90 main.f90というファイルがはいっています。
Lapackを使うのでそれにリンクしてコンパイルしてください。intel®fortranならば、

ifort -mkl integrate.f90 main.f90

で良いと思います。
実行すると、

$ ./a.out
 Solving TISE...
     10 /    100
     20 /    100
     30 /    100
     40 /    100
     50 /    100
     60 /    100
     70 /    100
     80 /    100
     90 /    100
    100 /    100
     1.807[CPU sec]
                    1   1.0000086267300243    
                    2   2.0001412369101184    
                    3   2.0001696393476047    
                    4   3.0004659233630600    
                    5   3.0006690578629835    
                    6   3.0013176072105336    
                    7   4.0029644525394392    
                    8   4.0032488578583107    
                    9   4.0091516049507794
$

という結果が得られるでしょう。ここで、** / 100という数字はハミルトニアン行列要素の計算の経過を表しています。100は行列が\(100\times 100\)の行列であるという意味です。
1 1.0000086267300243
は、対角化後の固有エネルギーの一番低い状態から順に出力しています。
デフォルトでは二次元の調和振動子ですので、
\(
E_n=(n_x+1/2)+(n_y+1/2)
\)

の順になります。一番低いエネルギーから順に\(1,2,2,3,3,3,4,4…\)が解析解となります。

具体的な計算例


三角形型のポテンシャル

\(
\begin{equation}
V(x,y)=\left\{
\begin{array}{c}
0,~~~(y\lt x \cap x\lt a \cap y\gt -a) \\
10,~~~(otherwise)
\end{array}
\right.
\end{equation}
\)
ここで、\(a=3\).

2つの四角形のポテンシャル

\(
\begin{equation}
V(x,y)=\left\{
\begin{array}{c}
0,~~~(y\lt (x+b) \cap y \gt -(x+b) \cap y\lt -(x+b)+a \cap y \gt (x+b)-a) \\
0,~~~(y\gt (x-b) \cap y \lt -(x-b) \cap y\gt -(x-b)+a \cap y \lt (x-b)+a) \\
10,~~~(otherwise)
\end{array}
\right.
\end{equation}
\)
ここで、\(a=3\)