sikino のすべての投稿

Lagrange, プログラミングと数値計算, 物理学

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

コメントする

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

拘束がある場合の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}=\frac{\nabla f}{|\nabla f|}\)の\(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式の差異である

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

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

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

Fortran90によるプログラム

Fortran90で2次元平面を動く質点の束縛運動のプログラムを作ります。
そのプログラムは
https://slpr.sakura.ne.jp/qp/supplement_data/constraint.tar.gz
に置いてあります。

一応補足しておきますが、初期状態の位置は\(f(x, y, t)=0\)を満たすx,yでなければなりませんし、初期速度も\(f(x, y, t)=0\)の傾きと一致していなければなりません。
上の条件を満たさずとも何も警告は出さず、計算自体は行われますが、その計算は束縛運動ではなくなります。

上のプログラムの中で、変更すると思われるのは以下の2つです。
1. fp2d
2. fw2d

1. fp2dは、ポテンシャル\(V(x,y)\)を記述しており、デフォルトのプログラムでは、重力による位置エネルギー
\(
V(x,y)=mgy
\)
に設定しています。

2. fw2dは、束縛条件\(f(x,y, t)=0\)を記述しており、デフォルトのプログラムでは、
\(
f(x,y,t)=y – [2(cos(1.5*x)-1) + 0.1x^2]=0
\)
に設定しています。

もしも非保存力を入れたい場合は、grkの最後の方を変更してください。

数値計算では、
時間発展は刻み幅制御の陽的ルンゲクッタ、
gradient, Hesse行列、束縛条件に関する時間の二階微分は、Hyper-dual numbersを利用してほぼ厳密な値を得ています。

program main
  ! Author : sikino
  ! Date   : 2019/10/12 (yyyy/mm/dd)
  ! URL    : http://slpr.sakura.ne.jp/qp/constraint-condition3/
  use Parameters
  use Hyperdualmod
  implicit none
  integer::i,N,info,Nt
  double precision::t,h,tol,ta,tb
  double precision,allocatable::x(:),tgrid(:)
  external::grk
  double precision::E,Ek,Ep

  N=4
  allocate(x(1:N))

  ! Time range [ta, tb]
  ta = 0d0
  tb = 30d0
  ! Divide time range equally among Nt
  Nt = 201
  ! Tolerance of the RK method
  tol=1d-8

  ! Initial condition
  x(1)=0d0 ! x (ta)
  x(2)=4d0 ! x'(ta)
  x(3)=0d0 ! y (ta)
  x(4)=0d0 ! y'(ta)

  ! Time grid
  allocate(tgrid(1:Nt))
  do i=1,Nt
     tgrid(i) = (i-1)*(tb-ta)/dble(Nt-1)+ta
  enddo
 
  do i=2,Nt
     info = 0
     t = tgrid(i-1)
     h = (tgrid(i) - t)*0.1d0
     do while(info.le.0)
        call drkf45(grk,t,h,N,x,tgrid(i),info,tol)
     enddo
     
     ! Kinetic energy
     Ek = 0.5d0*mass*(x(2)**2+x(4)**2)
     ! Potential energy
     Ep = mass*g*x(3)
     ! Total energy
     E  = Ek + Ep
     write(10,'(8e25.10e3)')t,x(1),x(2),x(3),x(4),E,Ek,Ep
  enddo
 
  stop
end program main

subroutine fp2d(N,xH,fH)
  use Parameters, only:mass, g
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Potential

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  ! If conservative force,
  !  fH doesn't depend on the time.
  !-----------------------------
 
  fH = mass * g * y
 
  return
end subroutine fp2d

subroutine fw2d(N,xH,fH)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  type(Hyperdual),intent(in)::xH(1:N)
  type(Hyperdual),intent(out)::fH

  ! Wall of the function

  type(Hyperdual)::x,y,t
  x = xH(1)
  y = xH(2)
  t = xH(3)
  !-----------------

  fH = y - (2d0*(cos(1.5d0*x)-1d0) + 0.1d0*x*x)
!  fH = (y-5d0)**2 + x**2 - 25d0
 
  return
end subroutine fw2d

デフォルトのまま計算すると、fort.10に以下の出力が現れます。アニメーションは、同じファイルに入っているanime.pltで実行できます。

エネルギーは計算の範囲ではほとんどありません。全力学的エネルギーは、計算時間内でおおよそ5桁目が変化していました。

参考文献

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

Lagrange, 物理学

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

コメントする

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

単振り子


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

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

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

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

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

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

が得られます。連立微分方程式(31)の３本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の３つなので、方程式を立てて\(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\)
・計算手法
　陰的ルンゲクッタ法

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

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

参考文献

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

Lagrange, 物理学

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

コメントする

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

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

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

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

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

１つの質量を持った質点が、ある軌道上しか動けない運動に対する運動方程式を導きましょう。
ある軌道上しか動けないという条件は、時刻\(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)の３本の方程式に対して未知の関数は\(x(t),y(t),\lambda(t)\)の３つなので、方程式を立てて\(x(t),y(t)\)を求めることが出来ます。

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

です。変形すれば

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

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

と書けます。よって

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

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

です。

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


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

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

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

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

を得ます。方程式(18a),(18b),(18c)の３本の方程式に対して未知の関数は\(r(t),\theta(t),\lambda(t)\)の３つなので、方程式を立てて\(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その２ -fortran66の日記

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

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

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

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

ＳＯＲＴＰＡＣＫ（ＳＯＲＴｘＫ，ＳＯＲＴｘｙ，ＳＲＴＶｘｚ）　（スカラー又はベクトルデータの内部ソーティング）

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に近づくほど難しくなります。
なぜなら、原点付近を通過する際に、質点の導関数の変化が大きくなるからです。

使用したプログラムの説明は
陽的解法はhttps://slpr.sakura.ne.jp/qp/runge-kutta-ex/
陰的解法はhttps://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
４倍精度 q

プログラムの評価


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

  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次元の運動方程式の解法

12件のコメント

グリーン関数を用いて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通りの方法が考えられます。
どちらの方法でも同じ解を与えるので、どちらをえらんでも構いません。
数学的に解きやすい方法を選ぶのが賢い方法です。

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

１つ目の選び方

微分方程式を

として考えます。
解を

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

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

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

なので、解

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

２つ目の選び方

微分方程式を

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

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

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

Thanks

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

プログラミングと数値計算

-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”という名前でサブルーチンに入っています。

プログラム

こちらです。
https://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\)

プログラミングと数値計算, 物理学

1次元時間依存シュレーディンガー方程式の高精度な数値解法

2件のコメント

1次元時間依存シュレーディンガー方程式を確実に解くことを目的とします。

本稿の方針は
計算時間が掛かってもいいので、高精度に解く
ことを主眼においています。対象とする問題は、
ポテンシャルに空間的、時間的に不連続な点があっても高精度に解く
事を目的とします。

ポテンシャルが滑らかで、変なことが起こらないのならば、空間や時間を離散化して解く
クランク＝ニコルソン法による時間発展が簡単で優秀な方法です。

本稿では導出する過程では\(N\)次元の問題として扱っていますが、プログラムは1次元の問題のみです。

方針

空間に関する積分を適応型数値積分、時間に関する積分を適応刻み幅陽的ルンゲ＝クッタ法で行います。

時間依存しないシュレーディンガー方程式は適応型数値積分であるQUADPACKを使うことで求められます。
詳しくはこちらへ。

解法

原子単位系における時間依存シュレーディンガー方程式は

として書けます。
波動関数を

のように展開します。ここで、基底関数\(\varphi_n(\mathbf{r})\)はハミルトニアン\(H_0\)の固有関数で、以下の固有値問題

の固有値\(E_n\)に属する固有関数です。
この固有関数は以下の通り規格化されているとします。

関数で展開するのはハミルトニアンに含まれる空間に対する二階微分を解析的に行い、消す為です。
式(1)に代入すると、

となり、左から\(\varphi^*_m(\mathbf{r})\)を掛けて全空間で積分すると

という式が得られます。ここで、\(U_{m,n}(t)\)を

と書きました。

この\(U_{m,n}(t)\)は

の性質を満たします。すなわち、ポテンシャル\(V(\mathbf{r},t)\)が実数値関数であれば、\(U_{m,n}(t)\)はエルミート行列となります。行列形式で書くとすれば、

と表せられます(補遺1)。

ここでは行列形式で書いても見やすくなるだけで、今回は使わないので、式(6)を高精度に解くことに集中します。

高精度に解くために

この問題を解くために高精度に求めなければならない重要な部分は

  1. 空間に対する積分
    \(\displaystyle
    \int \varphi^*_m(\mathbf{r})f(\mathbf{r})\varphi_n(\mathbf{r})d\mathbf{r}
    \)
    と
  2. 時間に対する積分(連立微分方程式)
    \(\displaystyle
    \frac{d c_m(t;c_1,~\cdots,c_{N})}{dt}=F_m(t;c_1,~\cdots,c_{N})
    \)
    という箇所を、不連続箇所があっても高精度に求める手法を使えばよい訳です。

空間に関してはQUADPACKによる適応型の数値積分, 時間に関しては刻み幅制御のルンゲ=クッタ法を用いることで対処します。
各々の詳細は以下のページをご覧ください。
最速・高精度の数値積分
ルンゲクッタ法の説明と刻み幅制御

ここまで確実に解けると仮定すると、精度の不安は基底関数の数だけに依存します。
基底関数の数を積めば積むほど精度が高い計算が出来ますが、計算時間が莫大になります。
ここは計算機の性能と折り合いをつけるしかありません。

ここから考える数値計算上の基底関数\(\varphi_n(x)\)は計算区間の端でゼロになるsin基底関数で、区間[x_a,x_b]で定義され、
\(\displaystyle
\varphi_n(x)=\sqrt{\frac{2}{x_b-x_a}}\sin\left(n\pi\frac{x-x_a}{x_b-x_a}\right)
\)
です。これは、固有値問題
\(\displaystyle
-\frac{1}{2}\frac{d^2}{dx^2}\varphi_n(x)=E_n\varphi_n(x),~~E_n = \frac{1}{2}\left(\frac{n\pi}{x_b-x_a}\right)^2
\)
境界条件
\(\varphi_n(x_a)=\varphi_n(x_b)=0\)
の元での解となっています。

さて、定式化は終わりましたので、あとはプログラムして解くだけです。

莫大な計算時間

問題は計算時間が膨大にかかりすぎる事です。

計算時間を評価してみます。

・次元数を\(D\)
・1次元辺りの基底関数の数は\(N\)
・合計\(N^D\)元1次連立微分方程式を解く
・1つの連立微分方程式の右辺を計算するのに\(N^D\)個の積分
・適応刻み幅ルンゲクッタ法を\(N_{t}\)回繰り返す
・適応刻み幅ルンゲクッタ法は時間1ステップ辺り\(5\)回評価
・積分は空間の分点数\(N_{\mathbf{r}}^D\)回に比例
・ポテンシャルが対称行列であることを利用すれば\(1/2\)倍
・CPU個数を\(N_{\mathrm{cpu}}\)とし、完璧に並列計算が行われれば\(1/N_{\mathrm{cpu}}\)倍
・関数を1回呼び出すのに必要な時間を\(t_f\)

であるので、全計算時間（関数が評価される回数に比例する、と考える）は
\(\frac{5}{2N_{\mathrm{cpu}}}\times N^D_{\mathbf{r}}N_{t}N^{2D}t_f\)と求められます。

現実的な量で計算時間を見積もってみます。
10[原子単位]秒の時間発展を区間20[原子単位]メートルで計算することを考えます。

\(D=1,~ N_{\mathbf{r}}=1000,~N_{t}=2000,~N=50,~N_{\mathrm{cpu}}=1\)とすると、\(1.25\times 10^{10}\)回関数を呼び出さなければなりません。

手持ちのノートパソコン(1.9GHz、１コア使用)で、\(2.5\times 10^8\)回の関数\(f(x)=\sin(x)\)の呼び出しに掛かるCPU時間を測ると約7[cpu秒]かかりました。
ということは、関数を1回呼び出すのにかかる時間\(t_f\)は約\(3\times 10^{-8}\)秒です。

よって、\(1.25\times 10^{10}\)回の関数の呼び出しだけで約400秒、すなわち8分かかります。

現実には他の処理も入るので、8分は絶対にかかる、という事です。
仮にCPUが10個あれば40秒、CPU速度が倍程度の物を使えば20秒位になるでしょうか。
まぁまぁ受け入れられる時間になります。

2次元を計算しようと思ったら基底関数が二乗、かつ空間の分点も二乗になるので、CPUが10個,4GHz近くのCPUを使えても\(20\times 1000\times 50=1000000\)秒程度、すなわち11日です。

11日の計算は現実的ではありません。
二乗で効いてくる基底関数の数を減らすか、空間の積分の精度を犠牲にする、CPUを増やす等の対策をしなければなりません。

また、ここで示した時間は最低かかる時間ですので、実際に動かすのであればこれ以上は確実にかかるでしょう。関数の中の関数を呼び出す等の操作があれば倍々に増えていくのです。

計算量を減らす工夫

ポテンシャルが変わらない場合、\(U_{m,n}(t)\)の積分は、前の時刻と値が変わらないのでいちいち計算するのは無駄です。

数値的に”ポテンシャルが変わっていない”という情報を取り出すことが出来れば、この無駄を排除することが出来ます。

”ポテンシャル関数が変わっていない”というのはどうすれば評価できるのでしょうか？
私が考えたのは、Quadpackの自動積分を行う時に評価された引数の点における関数の値が全て一致すれば関数が変わっていない、と判断するという評価方法です。

実際にこれはうまく行きまして、ポテンシャルが一度だけ切り替わる、という問題に対して計算時間が
1085[CPU sec] → 25[CPU sec]
にまで減少しました。約1/50の計算時間になりました。

先ほどの11日で～と言っていたのはこの工夫をし無ければ、ですので、ポテンシャルの変化が殆どない問題に限れば1時間程度で計算が終わる、ということです。

もちろん、ポテンシャルが変化し続ける場合も試しましたが、変化はほぼありません。
83[CPU sec] → 84[CPU sec]
程度になりました。
この場合は2次元の例で出した通り、11日間掛かる推測です。

初期状態の準備

式(6)を解こうとしても初期条件が与えられなければ数値的に解くことは出来ません。
初期条件で良く使われるのは、1)解析的な解を与える場合と 2)系の固有状態を与える場合の２種類があります。

  1. 解析的な解を与える場合

    この場合、時間発展を考える場合、上記の時間発展方法ではsin関数基底に射影して係数\(c_n(t)\)を求めなければなりません。
    時刻\(t=t_a\)で初期状態\(\psi(x,t_a)\)が解析的に与えられた場合、係数\(c_n(t)\)は、

    として与えられます。

  2. 系の固有状態を与える場合

    この場合は時刻\(t_a\)のポテンシャル中の時間依存しないシュレーディンガー方程式
    \(\displaystyle
    \left[H_0 + V(\mathbf{r},t_c)\right]\phi_n(\mathbf{r})=E’_n\phi_n(\mathbf{r})
    \)
    を実際に解かなければなりません。
    そのためには、時間発展に使う基底と同じ基底\(\varphi_n(\mathbf{r})\)で展開し、対角化するのが賢い方法です。
    すなわち、

    として展開し、ハミルトニアンを行列表示にし、対角化して固有値、固有ベクトルを得ればいいのです。

計算の解析

計算結果を解析するにあたり、ある時刻の固有関数系の基底にどの位存在確率があるかという情報が良く使われます。

初期状態で利用した固有状態の存在確率を調べましょう。
今、求めたいのは時刻\(t_c\)のポテンシャル中の固有関数で展開した時の係数\(a_m(t)\)を求める事です。

ここで、\(\phi_n(\mathbf{r})\)は時刻\(t_c\)のポテンシャル中の固有関数であり、

を満たしています。また、\(\varphi_n(\mathbf{r})\)は数値計算を行う際の基底関数であり、

を満たしています。
任意の時刻\(t\)における波動関数を

として表したいので、\(a_n(t)\)について求めていけば、

という関係式を導くことが出来ます。

数値計算プログラム

プログラムは以下のリンクにあります。
https://slpr.sakura.ne.jp/qp/supplement_data/highprecision_1Dtdsetise.tar.gz
そのままのプログラムでは、以下の設定となっています。
プログラムの利用に関してはhttps://slpr.sakura.ne.jp/qp/about/をお読みください。

解凍すると、4つのファイル(main.f90, integrate.f90,anime.plt,anime_project.plt)が入っています。計算パラメータの設定は、main.f90に入っているので各自調整してください。

繰り返しますが、全て原子単位系です。デフォルトでは、

計算時間の範囲\([ta,tb]\) … \([-5, 20]\)秒
出力時間枚数\(Nt\) … \(100\)枚
1ステップ当りの時間積分の精度\(\mathrm{RKeps}\) … \(10^{-5}\)
1ステップ当りの時間の最小刻み幅\(\mathrm{hmin}\) … \(10^{-4}\)
※理想は最小刻み幅は倍精度ならば\(10^{-14}\)程度が良いと思いますが、時間変化するポテンシャルを計算しようとする時、いつまでも最小刻み幅で進む時間がたまにあります。こういう時間はさっさと抜け出たいので、それを考慮して設定できるようにしています。

位置空間の範囲\([xa,xb]\) … \([-10, 10]\)
出力時の分点数\(Nx\) … 200点
時間積分の精度\(\mathrm{QPeps}\) … \(10^{-6}\)

初期状態を時間依存しないシュレーディンガー方程式を解くか、解析解で与えるかiname … “tise”
inameが”tise”の場合、どの時刻\(t_{ini}\)の固有状態を使うか … \(t_a\)

基底関数の数\(N\) … \(30\)個
計算に用いるスレッド数\(N{\mathrm{cpu}}\) … 1つ
空間積分の分点数を強制的に増やすか否かQPforce … 0

初期状態波動関数… 時間依存しないシュレーディンガー方程式の基底状態

ポテンシャル\(
\begin{eqnarray}
V(x,t)=\left\{
\begin{aligned}
0.5 x^2 (t\lt 0)\\
0.1 x^2 (t\ge 0)
\end{aligned}
\right.
\end{eqnarray}\)

で解くプログラムとなっています。

lapackと、必要があればopemMPを使います。

ifort -mkl -openmp integrate.f90 main.f90
./a.out

もしくは、gfortranでlapackにリンクしてintegrate.f90 main.f90を一緒にコンパイル、実行してください。
実行後、以下のファイルを生成します。

  • wf_***.d
    波動関数の各時刻におけるデータです。時間順にwf_0.d(初期状態), wf_1.d, wf_2.d,…という風に指定された等間隔の時刻で出力します。アウトプットのファイル数はNtで指定した数です。
  • project_ratio.d
    各時刻でのsin関数基底への射影した際の式(2)の係数\(c_n(t)\)の絶対値二乗を出力します。最も高い状態が振幅を持ってしまう場合、基底数が足りていないことを示しています。\( |c_n(t)|^2\)の事です。
  • unitarity.d
    ユニタリー性を確認します。要は数値計算の発散とかで全計算区間内の存在確率が変化していないか見るためのものです。\(\sum_{n=1}^N |c_n(t)|^2\)の事です。
  • timerange.d
    時間グリッドを確認します。物理的な意味は無く、どの時間で計算負荷が大きいかを見るためのものです。
  • timeindex_correspondence.d
    時間インデックス（wf.***.dの***の値と、原子単位系の時刻の変換）を出力します。

もしも初期状態を時間依存しないシュレーディンガー方程式を解いて準備していた場合、以下のファイルを追加で出力します。

  • initial_eigenvalue.d
    初期状態のエネルギー固有値を出力します。式(14)の\(E’_n\)です。
  • project_tdse_to_tisebasis.d
    各時刻における初期状態の固有状態に射影した時の係数の絶対値の二乗を出力します。式(13)の係数\(a_n(t)\)の絶対値二乗である\(|a_n(t)|^2\)を出力します。

波動関数をgnuplot上でアニメーションとして見たい場合、同封してあるファイル”anime.plt”をgnuplot上でloadすれば良いです。

gnuplot> load "anime.plt"

また、gifアニメが欲しい場合、

gnuplot> load "anime_project.plt"

とすれば、波動関数と時間依存しないシュレーディンガー方程式を解いた時の基底関数への係数の絶対値二乗を共に出力します。

実行例

デフォルトの設定で計算を行います。
実行すると

$ ./a.out
0 / 100        1.0000000000
10 / 100        1.0000001064
20 / 100        1.0000001082
30 / 100        1.0000012466
40 / 100        1.0000079960
50 / 100        1.0000082264
60 / 100        1.0000114047
70 / 100        1.0000154242
80 / 100        1.0000154993
90 / 100        1.0000193449
     1.948[CPU sec]
$

という文が出力されます。
左の/100とは、今計算が100枚中枚まで終わった、ということを示しており、枚数は変数Ntに対応しています。
右の1に近い数字は、全空間の存在確率を表しており、計算が発散しないかどうかを確認できます。これはルンゲクッタ法の要求精度\(RKeps\)に大きく依存します。
最後のCPUsecは計算に要したCPU時間であり、1コアではおおよそ実時間と同一です。なので、この計算は約2秒で終わる、ということです。ちなみにCPUは1.9GHzの物です。

anime_project.pltによってgifアニメを生成すると、以下のgifアニメが得られます。

その他の出力されるデータは以下の通り。

計算例

ポテンシャルによる反射


※この上図の縦軸は数値計算上の基底である\(\varphi_n(x)\)への射影した時の係数の絶対値二乗です。

ラビ周波数


基底状態と第一励起状態の間の周波数でポテンシャルを揺らしています。

意地悪なポテンシャル


こんな変なポテンシャルでも計算はできます、というデモです。

補遺1

もし仮に、\(U_{m,n}(t)\)が時間によって変化しないのであれば、
\(\displaystyle
\frac{\partial \mathbf{x}}{\partial t}=\mathbf{A}\mathbf{x}
\)
の形の微分方程式であるので、指数関数の解を仮定して\(\mathbf{A}\)を対角化、その固有値固有ベクトルを求めて、線形結合で表現すればokです（線形代数II/連立線形微分方程式等にあるのでそちらを参照してください）。