sikino のすべての投稿

Schrödinger, Heisenberg, Interaction描像

Schrödinger描像、Heisenberg描像、Interaction描像というものがあります。

Schrödinger描像は波動関数による表現方法で、多くの場合で分かりやすい表現なため良く取り入れられます。
ただし、時と場合によってHeisenberg描像の方が分かりやすかったり、Interaction描像の方が最適、というときがあります。

表現方法が違うだけで同じものを表現します。

  1. Schrödinger描像(S描像)

    \(\hat{H}\)が時間依存しない場合、時刻\(t\)での波動関数は初期状態(\(t=0\))の波動関数\(\Psi_S(0)\)を用いて、
    \(
    \displaystyle \Psi_S(t)=e^{-i\frac{\hat{H}}{\hbar}t}{\Psi_s(0)}\ \ \ \ (1)
    \)

    と書くことができます。ここで添え字”\(S\)“はS描像であることを強調しています。
    \(\displaystyle \hat{U}(t)=e^{-i\frac{\hat{H}}{\hbar}t}\ \ \ (2)\)
    と表記されることが多く、時間発展演算子と呼ばれています。
    今、(2)の両辺を時間微分すると、\(\hat{U}(t)\)の満たす微分方程式
    \(
    \displaystyle i\hbar \frac{d\hat{U}(t)}{dt}=\hat{H}\hat{U}(t), \ \ \ \hat{U}(0)=1\ \ \ (3)
    \)

    が得られます。また、式(3)は実は\(\hat{H}\)が時間依存していても成立します。
    この意味は、式(3)が時間発展演算子のより一般的な記述方法であることを意味しています。

    今、Schrödinger描像で、とある演算子\(\hat{A}\)の期待値を考えます
    (\(\hat{A}_S\)は時間依存してもしなくてもok)。
    (例えば\(\hat{A}_S\)は位置空間で、位置演算子だったら\(\hat{x}=x\)、運動量演算子だったら\(\hat{p}=-i\hbar\frac{d}{dx}\)です。)
    すると、その期待値を表す関数は時間依存し、時刻\(t\)の\(\hat{A}\)の期待値\(\langle\hat{A}\rangle_t\)は
    \(
    \begin{align}
    \langle\hat{A}\rangle_t
    &\displaystyle =\langle{\Psi_S}(t)|\hat{A}_S|\Psi_S(t)\rangle \\
    &\displaystyle =\langle e^{-i\frac{\hat{H}}{\hbar}t} \Psi_S(0)|\hat{A}_S|e^{-i\frac{\hat{H}}{\hbar}t}\Psi_s(0)\rangle
    \end{align}
    \)
    と書けます。これがSchrödinger描像における演算子の期待値を表現しています。

  2. Heisenberg描像(H描像)

    時刻\(t\)でのHeisenberg描像の波動関数\(\Psi_H(t)\)の定義は
    \(
    \begin{align}
    \Psi_H(t)&=e^{i\frac{\hat{H}}{\hbar}t}\Psi_S(t) \\
    &=e^{i\frac{\hat{H}}{\hbar}t}e^{-i\frac{\hat{H}}{\hbar}t}\Psi_S(0) \\
    &=\Psi_S(0)
    \end{align}
    \)
    です。
    Heisenberg描像で、\(\hat{A}\)の期待値\(\langle\hat{A}\rangle_t\)を考えましょう。
    \(
    \begin{align}
    \langle\hat{A}\rangle_t&=\langle\Psi_S(t)|\hat{A}_S|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_S(t)|e^{-i\frac{\hat{H}}{\hbar}t}e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}e^{i\frac{\hat{H}}{\hbar}t}|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_H(0)|\hat{A}_H(t)|\Psi_H(0)\rangle
    \end{align}
    \)
    すなわち、S描像の演算子\(\hat{A}_S\)との間には、
    \(
    \hat{A}_H(t)=e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}
    \)

    という関係があるのです。

    では、H描像で演算子\(\hat{A}_H(t)\)の時間依存性はどうなるのでしょう。
    上の式を両辺微分して\(i\hbar\)を掛けると、
    \(
    \begin{align}
    i\hbar\frac{d}{dt}\hat{A}_H(t)&=
    -\hat{H}e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}
    +e^{i\frac{\hat{H}}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}}{\hbar}t}\hat{H} \\
    &=-\hat{H}\hat{A}_H(t)+\hat{A}_H(t)\hat{H} \\
    i\hbar\frac{d}{dt}\hat{A}_H(t)&=\left[\hat{A}_H(t), \hat{H}\right]
    \end{align}
    \)

    もしも、\(\hat{A}_S\)が時間依存する形だったら
    \(
    \displaystyle i\hbar\frac{d}{dt}\hat{A}_H(t)=\left[\hat{A}_H(t), \hat{H}\right]+i\hbar \frac{\partial \hat{A}_H(t)}{\partial t}
    \)
    と書けます。

  3. Interaction描像(I描像)

    前提がありまして、系のハミルトニアン\(\hat{H}\)を\(\hat{H}=\hat{H}_0+\hat{V}\)と分けます。
    そして、\(\hat{H}_0\)の時間発展は分かっているものとします。

    この時、時刻\(t\)でのInteraction描像の波動関数\(\Psi_I(t)\)の定義は
    \(
    \begin{align}
    \Psi_I(t) &= e^{i\frac{\hat{H_0}}{\hbar}t}\Psi_S(t) \\
    &(= e^{i\frac{\hat{H_0}}{\hbar}t}e^{-i\frac{\hat{H}}{\hbar}t}\Psi_S(0) )
    \end{align}
    \)
    です。
    このInteraction描像の波動関数の意味を簡単に説明すると、
    \(\Psi_I(t)\)は、ある時刻で、S描像の波動関数を\(e^{-i\frac{\hat{H}}{\hbar}t}\)に従って時間を進める方向に時間発展させる部分と、\(e^{i\frac{\hat{H}_0}{\hbar}t}\)に従って時間を戻す方向に時間発展させる部分の2つがある、ということです。
    詳しくは後ほど述べます。

    I描像での演算子\(\hat{A}\)の期待値\(\langle\hat{A}\rangle_t\)はどうなるのでしょう。
    \(
    \begin{align}
    \langle\hat{A}\rangle_t&=\langle\Psi_S(t)|\hat{A}_S|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_S(t)|e^{-i\frac{\hat{H}_0}{\hbar}t}e^{i\frac{\hat{H}_0}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}_0}{\hbar}t}e^{i\frac{\hat{H}_0}{\hbar}t}|\Psi_S(t)\rangle \\
    &\displaystyle =\langle\Psi_I(t)|\hat{A}_I(t)|\Psi_I(t)\rangle
    \end{align}
    \)
    と書けます。

    すなわち、I描像では演算子が時間依存し、波動関数も時間依存します。
    よって、
    Interaction描像での演算子\(\hat{A}_I(t)\)の時間依存を記述する方程式と、
    Interaction描像での波動関数\(\Psi_I(t)\)の時間依存を記述する方程式
    があり、それらを求めてみましょう。

    • 演算子\(\hat{A}_I(t)\)の時間依存

      \(
      \hat{A}_I(t)=e^{i\frac{\hat{H}_0}{\hbar}t}\hat{A}_S e^{-i\frac{\hat{H}_0}{\hbar}t}
      \)
      演算子の期待値は上式の通りであり、もしも\(\hat{H}_0\)が時間依存しない場合、Heisenberg描像で\(\hat{H}\to \hat{H}_0\)と置き換えたものと同じ形になります。すなわち、
      \(
      \displaystyle i\hbar\frac{d}{dt}\hat{A}_I(t)=\left[\hat{A}_I(t), \hat{H}_0\right]
      \)
      となります。もしも\(\hat{H}_0\)が時間依存する場合は、
      \(
      \displaystyle i\hbar\frac{d}{dt}\hat{A}_I(t)=\left[\hat{A}_I(t), \hat{H}_0\right]+i\hbar \frac{\partial \hat{A}_I(t)}{\partial t}
      \)

    では、I描像で波動関数\(\Psi_I(t)\)が満たす運動方程式はどうなるでしょう。
    I描像の波動関数の定義式に左から\(e^{-i\frac{\hat{H_0}}{\hbar}t}\)を掛けると、
    \(
    \Psi_S(t) = e^{-i\frac{\hat{H_0}}{\hbar}t}\Psi_I(t)
    \)
    であり、これをS描像のSchrödinger方程式に代入します。
    \(
    \begin{align}
    & i\hbar \frac{\partial}{\partial t}\left(e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)\right)
    =\hat{H}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t) \\
    & i\hbar \left[-\frac{i}{\hbar}\hat{H}_0 e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)+e^{-i\frac{\hat{H}_0}{\hbar}t}\frac{\partial \Psi_I(t)}{\partial t}\right]
    =\hat{H}_0e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)+\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t) \\
    & i\hbar e^{-i\frac{\hat{H}_0}{\hbar}t}\frac{\partial \Psi_I(t)}{\partial t}=\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t) \\
    \end{align}
    \)
    左から\(e^{i\frac{\hat{H}_0}{\hbar}t}\)を作用させると、
    \(
    \displaystyle i\hbar \frac{\partial \Psi_I(t)}{\partial t}=e^{i\frac{\hat{H}_0}{\hbar}t}\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}\Psi_I(t)
    \)
    ここで、Interaction描像の演算子\(\hat{H}_I(t)\)を
    \(
    \displaystyle \hat{H}_I(t)=e^{i\frac{\hat{H}_0}{\hbar}t}\hat{V}e^{-i\frac{\hat{H}_0}{\hbar}t}
    \)
    と置くと、波動関数\(\Psi_I(t)\)が満たす運動方程式
    \(
    \displaystyle i\hbar \frac{\partial \Psi_I(t)}{\partial t}=\hat{H}_I(t)\Psi_I(t)
    \)
    が得られます。

さて、定義が終わったところで、Heisenberg描像とInteraction描像の意味を見てみましょう。
簡単にそれぞれの描像を説明すると、

Heisenberg描像量子-古典対応を成す
Interaction描像既に分かっている部分を取り除く

という特徴があります。

まずはHeisenberg描像から。
Heisenberg描像でHeisenberg描像の演算子\(\hat{A}_H(t)\)の時間発展を記述する運動方程式は
\(
i\hbar\frac{d}{dt}\hat{A}_H(t)=\left[\hat{A}_H(t), \hat{H}\right] \cdots(A)
\)
です。
今、\(\hat{A}_H(t)\)が位置演算子\(\hat{q}\)と運動量演算子\(\hat{p}\)によって表記されていると考えます。
この時、Schrödinger描像で、それぞれの演算子はハミルトニアンと
\(
\begin{align}
\left[\hat{q}_S,\hat{H}_S\right] &= i\hbar \frac{\hat{p}_S}{m} \\
\left[\hat{p}_S,\hat{H}_S\right] &= -i\hbar \frac{\partial}{\partial \hat{q}_S}\hat{H}_S
\end{align}
\)
という関係があります。
Heisenberg描像でSchrödinger描像のハミルトニアンがどう記述されるのかを考えてみると、
\(
\begin{align}
\hat{H}_S &= \frac{\hat{p}_S^2}{2m}+\hat{V}_S \\
&= \frac{1}{2m}e^{-i\frac{\hat{H}_S}{\hbar}t}\hat{p}_H^2 e^{i\frac{\hat{H}_S}{\hbar}t}+e^{-i\frac{\hat{H}_S}{\hbar}t}\hat{V}_H^2 e^{i\frac{\hat{H}_S}{\hbar}t} \\
&= e^{-i\frac{\hat{H}_S}{\hbar}t}\left( \frac{\hat{p}_H^2}{2m}+\hat{V}_H\right) e^{i\frac{\hat{H}_S}{\hbar}t} \\
\hat{H}_H&= e^{i\frac{\hat{H}_S}{\hbar}t}\hat{H}_S e^{-i\frac{\hat{H}_S}{\hbar}t} \\
&=\hat{H}_S
\end{align}
\)
となり、結局
\(
\hat{H}_S=\hat{H}_H
\)

であることが導かれます。
Heisenberg描像での位置演算子\(\hat{q}\)と運動量演算子\(\hat{p}\)の交換関係は
\(
\begin{align}
\left[\hat{q}_H,\hat{H}_H\right] &= i\hbar \frac{\hat{p}_H}{m} \\
\left[\hat{p}_H,\hat{H}_H\right] &= -i\hbar \frac{\partial}{\partial \hat{q}_H}\hat{H}_H \\
\end{align}
\)
であることが導けるため、\(\hat{q}_H\)と\(\hat{p}_H\)に関して、運動方程式に代入して
\(
\begin{align}
i\hbar\frac{d}{dt}\hat{p}_H(t)&=\left[\hat{p}_H(t), \hat{H}_H\right] \\
&= -i\hbar \frac{\partial}{\partial \hat{q}_H}\hat{H}_H \\
\frac{d}{dt}\hat{p}_H(t)&=-\frac{\partial}{\partial \hat{q}_H}\hat{H}_H
\end{align}
\)
であり、
\(
\begin{align}
i\hbar\frac{d}{dt}\hat{q}_H(t)&=\left[\hat{q}_H(t), \hat{H}_H\right] \\
&= i\hbar \frac{\hat{p}_H}{m} \\
\frac{d}{dt}\hat{q}_H(t)&=\frac{\partial}{\partial \hat{p}_H}\hat{H}_H
\end{align}
\)
となります。
この二つの式は、古典力学のハミルトンの運動方程式に酷似していますよね。ここからHeisenberg描像が量子-古典対応を成す、という所以になるわけです。

注意として、Heisenberg描像は、量子-古典対応を成しますが、演算子それ自身を測定することなどできないのです。あくまで期待値と対応するだけです。
また、Heisenberg描像は古典力学への回帰を見出し、エーレンファストの定理と関連します。


続いて相互作用描像です。
系のハミルトニアン\(\hat{H}\)を\(\hat{H}=\hat{H}_0+\hat{V}\)と分けます。
そして、\(\hat{H}_0\)の時間発展は分かっている時、Interaction描像の波動関数\(\Psi_I(t)\)は
\(
\begin{align}
\Psi_I(t) = e^{i\frac{\hat{H_0}}{\hbar}t}e^{-i\frac{\hat{H}}{\hbar}t}\Psi_S(0)
\end{align}
\)
でした。じっくり見ていきましょう。
例を考えます。
\(
\displaystyle \hat{H}=\hat{H}_0+V(x)=-\frac{\hbar^2}{2m}\nabla^2+V(x)
\)
で、ポテンシャル\(V(x)\)が単調減少で図のような山なりの波形をしている場合を考えます。
Interaction描像イラスト_c
初期状態\(\Psi_S(0)\)があったとすると時間経過すると、その波動関数の形は広がりながら右へ移動します。
Interaction描像では、この、広がる、という動きと、右へ移動する、という動きが合わさったと考えます。
広がる動きは、ポテンシャルが無いときの波束の動きです。ポテンシャルがない時、右に動くことはしません。
右への動きは、ポテンシャルしか無いときの波束の動きです。質点を考えれば、ポテンシャルの山を転がり、右へ動きます。

Interaction描像の波動関数\(\Psi_I(t)\)は既に分かっている広がる動き(\(\hat{H}_0=-\frac{\hbar^2}{2m}\nabla^2\)の部分)を取り除き、波束が右に動く描像(\(V(x)\)の部分)だけを切り抜くのです。

言葉で言えば、Interaction描像は、
\(t’\)秒後の\(I\)描像における波動関数は、\(t’\)秒後の\(S\)描像における波動関数に対して、\(t’\)秒間の\(\hat{H}_0\)による逆方向の時間発展を行えばよい、ということです。

あくまでもイメージなので、細かな議論はご了承を。

参考文献

David J. Tannor著、山下晃一ほか訳『入門 量子ダイナミクス(上)』(2011), p.231~236

gnuplotのカラーマップ

配色(pm3dの色)について
gnuplot ver4.6で確かめています。

一番良いと思ったカラーマップは、matlabで使われている”jet”だと感じました。

set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

colorjet
参考:Matlab colorbar with Gnuplot

※一番下に黒または白を追加した以下の物も良いと思います。

set palette defined ( -1 '#000030', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set palette defined ( -1 '#ffffff', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

※正負をはっきりと見せたい時。

set palette defined ( 0 '#0fffee',1 '#0090ff', 2 '#000fff',3 '#000090',4 '#ffffff',5 '#7f0000', 6 '#ee0000', 7 '#ff7000', 8 '#ffee00')

set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#ffffff',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

追記)2023/01/13
matlabでは、最近parulaというカラーマップが採用されているようです。

set palette defined ( 0 '#352a87',1 '#2053d4',2 '#0d75dc',3 '#0c93d2',4 '#07a9c2',5 '#38b99e',6 '#7cbf7b',7 '#b7bd64',8 '#f1b94a',9 '#fad32a',10 '#f9fb0e')

グレースケールにした時にも画像を変更する必要もなく、このままで綺麗に出るようです。
また、\(f(x,y)=\sin(xy)\)をjetと比較してみましょう。

グレースケールとした場合は、以下のようになります。

確かにparulaの方が高い低いがはっきりと分かります。

カラーの場合は偏りが生じないjetの方が綺麗ですかね?
しかし、parulaの方が柔らかな印象です。
グレースケールで印刷されることが予想されるならば、parulaの方がいいですね。
印象が大きく変わるかと思います。好きな方をご使用ください。


配色について。
gnuplotでは色の指定方法として3種類あるそうです((2),(3),(4))。

お勧めのカラーの指定順に,
(1)外部サイトgnuplottingにあるカラーマップを利用する
(2)cubehelixで指定する
(3)definedで指定する
(4)rgbformulaeで指定する

方法を紹介します。

(1)gnuplottingにあるカラーマップの利用


一番早く、綺麗なカラーマップを得るためには”gnuplotting(http://www.gnuplotting.org/)”という外部サイトに紹介されているgnuplotのカラーマップを使うことです。
ここに紹介されているカラーマップは、https://github.com/Gnuplotting/gnuplot-palettesより画像を引用して、
overview((c)gnuplotting)です。

使い方
  1. Gnuplotting/gnuplot-palettes -githubに飛びます。
  2. 右上にある “Clone or download” をクリックしてすべてのファイルをダウンロード
  3. カラーマップを見て使いたい名前の “***.pal” ファイルをgnuplot上でload
    で使えます。線の色なども上書きされるので、もしも線の色だけを変えたければ、linestyleなどで変更する必要があります。

(2)cubehelixで指定する


gnuplot version4.6以上

set palette cubehelix start 0.5 cycles -1.5 saturation 1

この色の指定方法はD.A.Greenによって発表された、宇宙の、強度イメージを色付けするために作られた色の指定方法です[1]。

それぞれのパラメータが何を表すかは折りたたんでおきますので、気になる方は下のリンクを開いてください。

[adsense1]

cubehelixでの具体例

対象とするデータは、水素原子の(n,l,m)=(4,2,0)の解析解です。
※このデータが欲しい方はこちらHydrogen_xz_nlm420.txt(1.5MB)に置きましたので、コピペしてください。
Ctrl+aで全選択できます。

3Dの表示ではこういうグラフです。
hydrogennlm420

特に記述が無いものはgammaの値はデフォルトの1.5です。

1.hydro420_s0.5c-1.5s1g1.5_c

set palette cubehelix start 0.5 cycles -1.5 saturation 1

きちんとした裏付けの下カラーリングされているので、論文等々にお勧めです。


2.hydro420_s0.5c-1s3g1.5_c

set palette cubehelix start 0.5 cycles -1 saturation 3

もしくは

 set palette defined(0"#000000",1"#5b13b8",2"#1c60ff",3"#00b7d8",4"#01eb75",5"#55f02e",6"#ccd73b",7"#ffc68c",8"#ffd5e2",9"#ffffff")

3.hydro420_s0.5c-1s3g3_c

set palette cubehelix start 0.5 cycles -1 saturation 3
set palette gamma 3

もしくは↑に近い色で

set palette defined(0"#000000",0.1"#8e007a",2.1"#579dff",3.1"#2de4ff",4.1"#4bffab",5.1"#98fe7d",6.1"#e4eb8c",7.1"#ffe1c0",8.1"#ffe9f0",9.1"#ffffff")

4.hydro420_s1c0s2g1.5_c

set palette cubehelix start 1 cycles 0 saturation 2

5.hydro420_s2c0s3g1.5_c

set palette cubehelix start 2 cycles 0 saturation 3

6.hydro420_s2c1s2g1.5_c

set palette cubehelix start 2 cycles 1 saturation 2

7.hydro420_s3c0.5s3g1.5_c

set palette cubehelix start 3 cycles 0.5 saturation 3

8.hydro420_s1c5s3g1.5_c

set palette cubehelix start 1 cycles 5 saturation 3

9.hydro420_s1c0s2g1.5neg_c

set palette cubehelix start 1 cycles 0 saturation 2 negative

10.hydro420_s1c-1s2g1.5neg_c

set palette cubehelix start 1 cycles -1 saturation 2 negative

11.hydro420_s3c0.5s1g1.5neg_c

set palette cubehelix start 3 cycles 0.5 saturation 1 negative

12.hydro420_s3c-2s2g1.5neg_c

set palette cubehelix start 3 cycles -2 saturation 2 negative

13.hydro420_s1c5s3g1.5neg_c

set palette cubehelix start 1 cycles 5 saturation 3 negative

82.54.4 Cubehelix

(3)definedで指定する

set palette defined(-3"blue",0"white",5"red",10"yellow")

z軸が-3になったとき青色(blue)、0になったとき白(white),5になったとき赤(red)、10になったとき黄色 (yellow・になるような配色です。
zの範囲が-50 < z < 100だった場合は-3,0,5,10の比で表されます。
すなわち、z=-50で青、z=-16付近で白、z=42付近で赤、z=100で黄色となるように自動的に調整されます。
色は有名な色は上の例のように言葉で記述できますが、16進数のカラーコードによる表示もできます。

カラーマップ例

カラーマップ例を載せます。


set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

matlabで使われている”jet”のカラーマップです。
colorjet
参考:Matlab colorbar with Gnuplot

set palette defined(0"#ff0000",1"#ff8000",2"#ffff00",3"#80ff00",4"#00ff00",5"#00ff80",6"#00ffff",7"#0080ff",8"#0000ff",9"#8000ff",10"#ff00ff",11"#ff0080",12"#ff0000")

cyclic_c
色は色相環を元に作っています。
周期的な関数を出力するときに役に立ちます。例えば、リーマン面を角度に依存させて表示する際などに利用すると良いです。こんなように。
Riemann_angle


set palette defined(1"#0000ff",2"#0080ff",3"#00ffff",4"#00ff80",5"#00ff00",6"#80ff00",7"#ffff00",8"#ff8000",9"#ff0000")

cl2_c
色相環を元に作成したものです。


set palette defined (-9 "purple",-6 "blue",-3 "turquoise", 0 "#f5f5f5", 3 "gold", 6 "orange", 9 "red")

cwhite
※匿名希望さん提供


set palette defined(0"#ffffff",0.8"#00008b",1.8"#2ca9e1",3"#008000",4.2"#ffff00",5"#eb6101",5.5"#8b0000")

カラフルtr_c
白⇒青⇒水色⇒黄色⇒オレンジ⇒赤
細かい変化を見たい時にいい感じになります。


set palette defined(0"#000000",1"#0000ff",2"#1e90ff",3"#00ffff")

黒から青tr_c
黒⇒青⇒明るい青
不連続な関数を表示させる時にこのマップを使うといい感じになります。


set palette defined(0"#ffffff",1"#0000ff",2"#ff0000",3"#ffff00")

白青黄tr_c
白⇒青⇒赤⇒黄
上の2つが気に入らなかったりした時に。


set palette defined(0"#00008b",1"#2ca9e1",2"#38b48b",3.5"#ffff00",5"#eb6101",5.3"#c9171e")

カラフルatr_c
青⇒水色⇒黄色⇒オレンジ⇒赤
一番下が白じゃないほうがいいときに。

(4) rgbformulaeで指定する

set palette rgbformulae 33,13,10

”set palette rgbformulae” のあとの3つの整数で赤(R)/緑(G)/青(B)に対応しています。
整数で赤、緑、青をどのような変化具合をさせるのかということを指定します。
直線的に色を変化させるのか、2次、3次関数的なのかはたまたsin,cos的なのかを。どういった整数を使えばどういう関数に相当するのか?
これを知るには gnuplot上で
”show palette rgbformulae”
と打ってください。gnuplot ver.4.6では、以下の画像のように出力されます。
rgbformulae_c
0~36まで出てきます。指定は-36~36までできます。負の値は反転することを意味しています。

また、いくつかのrgbformulaeを使った例が ”help palette rgbformulae” で見れます。実際に入力するとこんな感じ。
helprgb_c
helpで現れるカラーマップはこんな感じです。
rgbformulae_help_c

また、今使っているカラーマップをrgbformulaeで記述をしたい場合、

show palette fit2rgbformulae

とやれば、今使っているカラーマップをrgbformulaeで表示する場合の3つの数字を出してくれます。
あくまで似た図です。

デフォルトのカラーマップを変更する

デフォルトのカラーマップを変更するには、ホームディレクトリにて、
.gnuplot
という名前のファイルを作り、その中に

set palette defined(0"#ffffff",0.8"#00008b",1.8"#2ca9e1",3"#008000",4.2"#ffff00",5"#eb6101",5.5"#8b0000")

と記述します。これでデフォルトのカラーマップが変更されます。

ちなみに、僕が使っている.gnuplotファイルは以下のものです。

set terminal wxt dashed noraise enhanced font 'Times New Roman,20'

#matlab jet
#set palette defined ( 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set palette defined ( -1 '#ffffff', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set linetype  1 lc rgb "black" lw 1
set linetype  2 lc rgb "red" lw 1
set linetype  3 lc rgb "blue" lw 1
set linetype  4 lc rgb "forest-green"
set linetype  5 lc rgb "magenta" lw 1
set linetype  6 lc rgb "#FFD900" lw 1
set linetype  7 lc rgb "cyan" lw 1
set linetype  8 lc rgb "#7F00FF" lw 1
set linetype  9 lc rgb "#FF7F00" lw 1
set linetype  10 lc rgb "#00FF7F" lw 1
set linetype  11 lc rgb "gray" lw 1
set linetype  12 lc rgb "gray40" lw 1

上の設定を適応すると、カラーマップはjetに黒を追加したもの、線は以下のようになります。

配色のあれこれ、gnuplot demoより。
カラフルに塗りたい http://gnuplot.sourceforge.net/demo/pm3dcolors.html
モノクロに塗りたい http://gnuplot.sourceforge.net/demo/pm3dgamma.html
点に色をつけたい http://gnuplot.sourceforge.net/demo/rgb_variable.html
線に色をつけたい http://gnuplot.sourceforge.net/demo/rainbow.html

http://ayapin-film.sakura.ne.jp/Gnuplot/Primer/Misc/colormap.html

[1]D. A. Green(2011), ‘A colour scheme for the display of astronomical intensity images’, Bulletin of the Astronomical Society of India, 39, 289.(2011BASI…39..289G at ADS.)

またDave Green’s `cubehelix’ colour schemeをご覧下さい。

[adsense2]

LinuxMintでアップデート失敗

僕の環境は、
VMware(R) Player バージョン7.1.0 build-2496824
に、LinuxMint 17.1 Rebeccaをインストールしている環境です。

sudo apt-get update
sudo apt-get dist-upgrade

をした後、再起動したあと、
LinuxMintの起動後にモードを選択する画面が表示され、

udevadm trigger is not permitted while udev is unconfigured.
udevadm settle is not permitted...
...
ALERT! /dev/disk/by-uuid/....-....-....-.... -.... does not exist.
Dropping to shell!
BusyBox v1.21.1 (Ubuntu 1:1.21.0-1ubuntu1) built-in shell (ash)
Enter 'help' for a list of built-in commands.

(initramfs) _

画像ではこんな感じになりました。
20150511-51639_b_c1
※uuidを晒しても安全か分からなかったため、念のため隠しています。

こうなってしまった後の場合の復帰方法を書きます。手順を書けば、

  1. liveCD からLinuxを起動
  2. “sudo fdisk -l” を入力し、ブートディスクを確認する, 僕の場合は /dev/sda1 でした。
  3. “sudo mkdir /media/newroot” を入力
  4. “sudo mount /dev/sda1 /media/newroot”, を入力, ブートディスクを作成したディレクトリに変更
  5. “sudo chroot /media/newroot” を入力
  6. “sudo apt-get update”
  7. “sudo apt-get dist-upgrade”
  8. 再起動で幸せに!

で治りました。
Thread: “udevadm trigger is not permitted while udev is unconfigured.”より。
ただし、6,7の手順ではすんなりアップデートはできず、apt-get update入力後失敗し、
>代わりにconfigure…を入力してください
と出て、こちらを改めて入力しました。
また、dist-upgradeもすんなりはできず、オプション追加して再入力してください、とあったのでその通りにしました。


詳しく書いていきます。

  1. liveCD からLinuxを起動

    必要なのは使っているLinuxの.isoファイルです。
    LiveCDはLinuxをインストールしたときに使用した,LinuxOSが入っているCD/DVD/USBのことです。
    ですが、物理的なCDドライブは特に必要ありません。.isoイメージファイルさえあればokです。
    LinuxMint17.1の場合はisoイメージファイルを
    Download Linux Mint 17.1 Rebecca
    よりダウンロードできます。日本の場合はno-codecs版です。
    ※ちなみに、VMwareでUSBからのブートは出来ません。Biosに項目があるかもしれませんが、できません。
    ダウンロードしたら、VMwareの設定で、起動時に仮想ドライブを起動するようにします。
    仮想マシンの設定を選択して、設定します。
    20150512-164029a_ca_c
    ↑を選択してから
    20150512-213544_c_c
    を選択してokを。

    そして、問題になっているLinuxを起動し、Bios画面を表示させ、CDromから起動するように変更します。
    起動の途中でEscを押せばBios画面に入れるのですが、一瞬で読み込みが終わり、起動画面は消えるので結構シビアです。
    確実にBios画面に入るためにはVMwareをCDブート、USBブートする/VMwareのBIOS画面を出すより、このようにすればいいそうです。

    仮想OSを作成したフォルダ内の「*.vmx」ファイルに、以下の1行を追加する事で、ESCキーを受け付ける時間を延ばす事が可能です。
    BIOS.BootDelay = “5000”
    単位はミリ秒。上記の場合、5秒間の猶予が出来ます。
    2~3秒に設定していると、意外にあっという間に通過します。お勧めは4~6秒くらいです。
    余談ですが、以下の一行を追加すると、毎回BIOS画面に入れます
    bios.forceSetupOnce = “TRUE”

    そうして起動します。するとGUIでLinuxMintが起動するかと思います。そうしたらターミナルを開いて次のステップに進みます。

  2. “sudo fdisk -l” を入力し、ブートディスクを確認する

    次に、今使われているブートディスクは何なのか、を調べます。

    sudo fdisk -l

    を入力するとこの画像のように表示されることでしょう。
    20150512-162004_c
    ここで、左から2列目、”Boot” とある欄を見ますとアスタリスクが付いていることが分かります。
    なので、これが自分のBootディスク、ということになります。
    上の画面の場合は/dev/sda1, ということです。

  3. “sudo mkdir /media/newroot” を入力

    説明はありません。

  4. “sudo mount /dev/sda1 /media/newroot” を入力

    説明はありません。

  5. “sudo chroot /media/newroot” を入力

    説明はありません。

  6. “sudo apt-get update”

    ここでちょっと躓きました。僕の場合だけかもしれませんが、エラーが出されました。
    > Err http://…
    みたいな文がずらずらと。ただ、よくよくエラー文を読むと、
    代わりにconfigure…を入力してください
    とあったのでこれを入力しました。そうしたら、とりあえずokとなりました。

  7. “sudo apt-get dist-upgrade”

    次にこれを入力します。ここでもエラーが出ましたが、エラー文に
    オプションを追加して入力してください、
    とあったのでその通りにしたらokとなりました。

  8. 再起動

    再起動すれば元通り!
    この時、Bios画面か、仮想CD/DVDドライブから読み込ませないようにしておかないと意味がないので外しておきましょう。

途中、sudo に関して、

sudo: unable to resolve host

が出てきました。一応解決するためにはこのリンク先通りにやれば解決します。
これが必要な手順か、は知りません。

カーネルが、という場合もあるようです。その時は手順6,7が違うものになるそうです。

参考先

Thread: “udevadm trigger is not permitted while udev is unconfigured.” このエラーのフォーラム
Download Linux Mint 17.1 RebeccaLinux Mint 17.1 Rebeccaのダウンロード
VMwareをCDブート、USBブートする/VMwareのBIOS画面を出すBIOS画面を確実に出すためには。
sudo: unable to resolve host が表示されたらsudoの解決

とある詳しい方。

演算子の種類と説明

\(\hat{A}\)がエルミート演算子ならば、\(e^{i\hat{A}}\)はユニタリー演算子である。

この文の意味が分かる人はこのページはいらないと思います。


前提として、一度勉強した人が思い出す、という体を想定しています。
詳しくブラケット表記であるとか、正しく知りたい人は、
J. J. Sakurai著 桜井 明夫訳『現代の量子力学〈上〉』 (物理学叢書) (1989)



David J.Tannor著 山下晃一訳 『入門 量子ダイナミクス 時間依存の量子力学を中心に(上)』 化学同人
などを参考にしてください。このページの参考先もこの2つです。

[adsense1]

エルミート演算子(Hermite Operator)


エルミート演算子は自己随伴演算子とも呼ばれます。

  • ディラックのブラ・ケット表記で考えます。
    演算子\(\hat{A}\)の状態\(u\)と状態vによる内積を\(\langle u|\hat{A}|v\rangle \)と、表現します。
    関数による表現では、
    \(
    \displaystyle \langle u|\hat{A}|v\rangle =\int u^* \hat{A} v dx
    \)

    となるわけです。この時、
    \(
    \langle u|\hat{A}|v\rangle^*=\langle v|\hat{A}^{\dagger}|u\rangle
    \)

    を関係を満たす演算子を随伴演算子と呼び、\(\hat{A}^{\dagger}\)と表現します(\(\hat{A}\)のエルミート共役を取る、とも言います)。

    状態による表記では、
    \(
    \displaystyle \hat{A}|v\rangle =\lambda |v\rangle
    \)

    に対してエルミート共役を取ると、
    \(
    \displaystyle \langle v|\hat{A}^{\dagger} =\lambda^* \langle v|
    \)

    と書けます。

    そして、たまたま\(\hat{A}^{\dagger}\)が\(\hat{A}\)に等しい場合、すなわち、
    \(
    \hat{A}^{\dagger}=\hat{A}
    \)

    を満たすとき、演算子\(\hat{A}\)はエルミート演算子(自己随伴演算子)だ、と呼びます。

  • エルミート演算子の持つ性質
    1. エルミート演算子の固有値は実数である。
      \(
      \hat{A}|v\rangle=\lambda |v\rangle
      \)

      左から\(v\)を掛けて、内積を取り、式変形します。式変形は2通り考えられて、

      \(
      \begin{eqnarray}
      \left\{
      \begin{aligned}
      \langle v|\hat{A}|v\rangle&=\langle v|\lambda |v\rangle=\lambda\langle v|v\rangle \\
      \langle v|\hat{A}|v\rangle&=\langle v|\hat{A}^{\dagger}|v\rangle=
      \langle v|\lambda^*|v\rangle=\lambda^*\langle v|v\rangle
      \end{aligned}
      \right .
      \end{eqnarray}
      \)

      となります。同じものから出発したので値は同じものになるはずです。なので
      \(\lambda=\lambda^*\)
      これを満たす\(\lambda\)は虚数部がゼロでなければなりません。
      よってエルミート演算子の固有値は実数である、となります。

    2. あるエルミート演算子が異なる2つの固有値を持つ場合、これらの固有値に対応する固有関数は互いに直交する。
      2つの固有値を\(\lambda_1, \lambda_2\)と書いて、
      それぞれの固有値に属する固有ベクトルを\(|v_1\rangle, |v_2\rangle\)と書くことにします。すなわち、
      \(
      \hat{A}|v_1\rangle=\lambda_1 |v_1\rangle \ ,\ \ \ \ \hat{A}|v_2\rangle=\lambda_2 |v_2\rangle
      \)

      であるとします。ここで\(\lambda_1\neq \lambda_2\)です。
      1番目の式に左から左から\(\langle v_2|\)を作用させて式変形します。式変形は2通り考えられて、
      \(
      \begin{eqnarray}
      \left\{
      \begin{aligned}
      \langle v_2|\hat{A} |v_1\rangle&=\langle v_2|\lambda|v_1\rangle=\lambda_1\langle v_2|v_1\rangle \\
      \langle v_2|\hat{A}|v_1\rangle&=\langle v_2|\hat{A}^{\dagger}|v_1\rangle=
      \lambda_2^*\langle v_2|v_1\rangle=\lambda_2\langle v_2|v_1\rangle
      \end{aligned}
      \right .
      \end{eqnarray}
      \)

      同じものから出発したので、
      \(
      \begin{align}
      \lambda_1\langle v_2|v_1\rangle=\lambda_2\langle v_2|v_1\rangle \\
      \rightarrow (\lambda_1-\lambda_2)\langle v_2|v_1\rangle=0
      \end{align}
      \)

      仮定より、\(\lambda_1\neq \lambda_2\)なので、\(\langle v_2|v_1\rangle=0\)になるほかありません。
      \(\langle v_2|v_1\rangle=0\)は内積がゼロ、すなわち直交である、と言っているので仮定は示されました。
  • ある状態\(\phi\)が、演算子\(A\)の固有値\(\lambda_i\)に属する固有状態\({v_i}\)の組で書かれるとき、
    すなわち、
    \(
    \displaystyle |\phi\rangle=\sum_ia_i|v_i\rangle
    \)

    で書かれるとき、
    \(
    \displaystyle \frac{\langle\phi|\hat{A}|\phi\rangle}{\langle\phi|\phi\rangle}
    =\sum_i p_i \lambda_i,\ \ \ p_i\equiv \frac{|a_i|^2}{\sum_i|a_i|^2}\ \ \cdots (a)
    \)

  • 量子力学の中心的教義は以下の3つを主張しています。
      (i)     全ての観測量にはエルミート演算子\(\hat{A}\)が結び付けられる。
      (ii)   \(\hat{A}\)に属する観測量の測定について、起こりうる結果は\(A\)の固有値\(\{A_i\}\)のみ。
      (iii)  系が状態\(|\phi\rangle=\sum_ia_i|v_i\rangle\)にあれば、
      \(\lambda_i\)の値を得る確率は\((a)\)で与えた\(p_i\)で与えられる。
      また、\(\displaystyle \sum_ip_i\lambda_i\)は\(\hat{A}\)を測定した平均値、もしくは期待値である。
  • おまけ
    随伴行列\(A^{\dagger}\)は、\(\hat{A}\)の転置行列の複素共役で与えられます。
    \(
    (A_{ij})^{\dagger}=A^*_{ji}
    \)

    特に、エルミート行列の場合、\(A^{\dagger}=A^*\)なので、
    \(
    A_{ij}=A^*_{ji}
    \)

    対角要素については\(A_{ii}=A^*_{ii}\)が成り立つので、エルミート行列の対角要素は実数でなければならないことがわかります。

逆演算子


\(\hat{A}\)の逆\(\hat{A}^{-1}\)を意味します。

  • 定義

    \(\hat{A}|u\rangle=|v\rangle\ \ \ \cdots (1)\)
    のとき、逆演算子\(\hat{A}^{-1}\)は
    \(|u\rangle=\hat{A}^{-1}|v\rangle\ \ \ \cdots (2)\)
    と定義されます。

  • 性質

    (1)の左から\(\hat{A}^{-1}\)を作用させると、
    \(\hat{A}^{-1}\hat{A}|u\rangle=\hat{A}^{-1}|v\rangle\ \ \ \cdots (1-1)\)
    (2)の左から\(\hat{A}\)を作用させると、
    \(\hat{A}\hat{A}^{-1}|u\rangle=\hat{A}|v\rangle\ \ \ \cdots (2-1)\)

    (1)と(2-1), (2)と(1-1)を見比べれば、明らかに
    \(\hat{A}\hat{A}^{-1}=\hat{A}^{-1}\hat{A}=\mathbf{1}\)
    となります。

  • 逆演算子の存在
    全ての演算子に逆演算子があるわけではありません。以下の通り、逆が存在しない場合があることを示せます。
    演算子\(\hat{A}\)が2つの異なる初期ベクトルを同じ終ベクトルに写す状況を考えます。式で表せば、
    \(
    \begin{eqnarray}
    \left\{
    \begin{aligned}
    \hat{A}|u_1\rangle&=|v\rangle \ \ \ (3)\\
    \hat{A}|u_2\rangle&=|v\rangle \ \ \ (4)
    \end{aligned}
    \right.
    \end{eqnarray}
    \)
    ここで、
    \(\hat{A}^{-1}|v\rangle\)を考えた時、それを\(|u_1\rangle\)か、\(|u_2\rangle\)かを決める術はありません。
    よって逆が存在しないことになります。
    また、 \((4)-(3)\)を行うと、
    \(
    (4)-(3)=\hat{A}(|u_1\rangle-|u_2\rangle)=\mathbf{0}
    \)
    となります。この意味は、\(\hat{A}\)は\(\mathbf{0}\)ではない、いくつかのベクトルを消去する、ということを表現しています。

[adsense2]

ユニタリー演算子


ユニタリー演算子は演算子\(\hat{A}\)の逆演算子\(\hat{A}^{-1}\)が
\(\hat{A}\)のエルミート共役\(\hat{A}^{\dagger}\)に等しいとき、その演算子はユニタリー演算子だ
、と定義されます。すなわち、
\(
\hat{A}^{-1}=\hat{A}^{\dagger}
\)
を満たすとき、と定義されます。
また、別の表現では
\(
\hat{A}^{\dagger}\hat{A}=\hat{A}\hat{A}^{\dagger}=\mathbf{1}
\)
という場合もありますが同じことです。
ユニタリー演算子はよく、\(\hat{U}\)という表記がされます。

  • ユニタリー演算子の性質
    1. ユニタリー演算子はノルムを保存する。
      ユニタリー演算子が状態\(|u\rangle\)に作用した場合を考えます。この時、ノルムは自身の内積を取ればいいので、
      \(\langle u|\hat{A}^{\dagger}\)と\(\hat{A}|u\rangle\)を作用させればノルムになります。故に、
      \(
      \langle u|\hat{A}^{\dagger}\hat{A}|u\rangle=\langle u|u\rangle
      \)
      であり、ノルムは変わりません。
    2. ユニタリー演算子の固有値は絶対値\(\mathbf{1}\)を必ず含む。
      あるユニタリー演算子を\(\hat{U}\)と書き、固有値問題を
      \(
      \hat{U}|v\rangle=\lambda |v\rangle\ \ \ (5)
      \)
      と書きます。

      両辺のエルミート共役をとって、
      \(
      \displaystyle \langle v|\hat{U}^{\dagger} =\lambda^* \langle v|
      \)

      ユニタリー演算子の性質を使うと、
      \(
      \displaystyle \langle v|\hat{U}^{-1} =\lambda^* \langle v|
      \)

      右から\(\hat{U}\)を作用させて、
      \(
      \begin{align}
      \displaystyle \langle v|\hat{U}^{-1}\hat{U} &=\lambda^* \langle v|\hat{U} \\
      \displaystyle \langle v|\hat{U} &=\frac{1}{\lambda^*}\langle v| \ \ \ \ (6)
      \end{align}
      \)

      (5)に左から\(\langle v|\)を作用させると、
      \(
      \langle v|\hat{U}|v\rangle=\langle v|\lambda |v\rangle=\lambda\langle v |v\rangle\ \ \ (7)
      \)

      であり、また、内積を以下のように変形し、(6)を使うと、
      \(
      \displaystyle \langle v|\hat{U}|v\rangle=\frac{1}{\lambda^*}\langle v|v\rangle\ \ \ (8)
      \)

      従って、式(7),(8)は同じものから出発したので等しいはずで、
      \(
      \begin{align}
      \lambda\langle v |v\rangle&=\frac{1}{\lambda^*}\langle v|v\rangle \\
      &\rightarrow |\lambda|^2=1
      &\rightarrow |\lambda|=1
      \end{align}
      \)
      となるため、ユニタリー演算子の固有値は必ず絶対値1を含みます
      ※これは、\(e^{i\theta},\ \ \ \theta\)は実数、であればいいと言っています。1である必要性はありません。

    3. 異なる固有値を持つユニタリー演算子の固有状態は直交する
      異なる2つの固有値を\(\lambda_1, \lambda_2\)と書き、
      それぞれの固有値に属する固有ベクトルを\(|v_1\rangle, |v_2\rangle\)と書くことにします。
      \(|v_2\rangle\)と\(|v_1\rangle\)による内積をそれぞれ考えると、
      \(
      \langle v_2|\hat{U}|v_1\rangle=\langle v_2|\lambda |v_1\rangle=\lambda_1\langle v_2 |v_1\rangle\ \ \ (9)
      \)

      と変形できるし、
      \(
      \langle v_2|\hat{U}|v_1\rangle=\frac{1}{\lambda_2^*}\langle v_2|v_1\rangle=\lambda_2\langle v_2|v_1\rangle \ \ \ (10)
      \)

      とも変形できます。最後の変形は2.の性質、絶対値1を持つことを利用しています。
      2つの固有値\(\lambda_1, \lambda_2\)は違う値を仮定したので、(9)=(10)が成り立つためには
      \(\langle v_2|v_1\rangle=0\)
      でなければなりません。よって、ユニタリー演算子の異なる固有値に属する固有状態は直交していなければなりません。

ユニタリー演算子とエルミート演算子


エルミート演算子からユニタリー演算子を作る方法を示します。
この方法は運動量演算子の導出でも用いるので、覚えておいて損はないかと思います。

\(\hat{A}\)がエルミート演算子ならば、\(e^{i\hat{A}}\)はユニタリー演算子である。
演算子\(e^{i\hat{A}}\)の逆行列を考えて、それがエルミート共役に等しくなるか、見てみます。
定数に対するエルミート共役は単なる複素共役、\(\hat{U}=\hat{U}^{\dagger}\)であることを利用すると、
\(
\begin{align}
\left(e^{i\hat{A}}\right)^{-1}&=e^{-i\hat{A}} \\
&=1+(-i\hat{A})+\frac{(-i\hat{A})^2}{2!}+\cdots \\
&=1+(i\hat{A})^{\dagger}+\frac{\left\{(i\hat{A})^{\dagger}\right\}^2}{2!}+\cdots \\
&=\left(e^{i\hat{A}}\right)^{\dagger}
\end{align}
\)
となり、ユニタリー演算子の定義\(\hat{U}^{-1}=\hat{U}^{\dagger}\)を満たしていることがわかります。

逆行列とユニタリー行列


行列表記と演算子表記は同じものです。
量子力学では、ユニタリー行列の逆行列に関心があります。数式ならば、
\(\hat{U}^{-1}_{ij}=\hat{U}^{\dagger}_{ij}=\hat{U}^{*}_{ji}\)
ということです。

ユニタリー行列の性質の一つに、ユニタリー行列の列ベクトルは正規直交ベクトルである、ことを示しましょう。
逆行列との積は\(\mathbf{1}\)に等しいので
\(
(\hat{U}^{-1}\hat{U})_{ik}=\delta_{ik}
\)
です。これを踏まえ、ユニタリー演算子の性質を使って、
\(
\begin{align}
(\hat{U}^{\dagger}\hat{U})_{ik}&=\sum_j\hat{U}_{ij}^{\dagger}\hat{U}_{jk} \\
&=\sum_j\hat{U}_{ji}^{*}\hat{U}_{jk}=\delta_{ik}
\end{align}
\)
という結果が得られます。最後の式から、行列\(U\)の2つの列ベクトルの内積、これが\(\delta_{ik}\)に等しいことがわかります。
故にユニタリー行列の列ベクトルは正規直交ベクトルであることが示されました。

ちなみに、逆行列を得るための一般的な手続きは、
①  行列を対角形に変形

②  対角要素を逆数に

③  逆変換
という流れで行われます。

時間発展演算子の導出

量子力学における時間発展演算子\(\displaystyle \hat{U}_{(t)}=e^{-i\frac{\hat{H}}{\hbar}t}\)の導出です。
これは、ハミルトニアンは時間依存しない場合に限ります。時間依存するときは時間順序積という概念が出てきます。


シュレディンガー方程式
\(
\begin{align}
i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}\psi_{(x,t)}
\end{align}
\)
は既知のものとして進めます。

形式的な導出でよく説明されるのは以下の手順です。

シュレディンガー方程式の両辺を\(i\hbar\)で割って、更に波動関数\(\psi_{(x,t)}\)で割り両辺を時刻\(0\sim t\)まで積分します。
すなわち、
\(
\begin{align}
\int_{\psi_{(0)}}^{\psi_{(t)}}\frac{d\psi}{\psi}&=\int_0^t\frac{\hat{H}}{i\hbar}dt^{\prime} \\
\ln\frac{\psi_{(t)}}{\psi_{(0)}}&=\frac{-i}{\hbar}\hat{H}t \\
\end{align}
\)

なので
\(
\displaystyle \psi_{(x,t)}=e^{-i\frac{\hat{H}}{\hbar}t}\psi_{(x,0)}
\)

証明終了,,

(;’;゚;ж;゚;`;)ブホォッ
(; ^ω^)は…

何で波動関数で割ってくれちゃってるんですか?ハミルトニアンは演算子ですよ。形式的にしてもひどくない?

ということで、ちゃんと求めましょう。結果は↑のものと同じになります。不思議ですね。


ハミルトニアンが時間依存しないという条件の下求めます。
まず、時間依存するシュレディンガー方程式より、
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}\psi_{(x,t)}
\)

となります。無限小時間\(\Delta t\)に対して、シュレディンガー方程式は、
\(
\displaystyle i\hbar\lim_{\Delta t\rightarrow 0}\frac{\psi_{(x,t+\Delta t)}-\psi_{(x,t)}}{\Delta t}=\hat{H}\psi_{(x,t)}
\)

です。以降\(\lim\)は省略します。

式を変形して、
\(
\displaystyle \psi_{(x,t+\Delta t)}=\left(1-i\frac{\hat{H}}{\hbar}\Delta t\right)\psi_{(x,t)}
\)

また、もう\(\Delta t\)だけ時間を進めると、
\(
\displaystyle \psi_{(x,t+2\Delta t)}=\left(1-i\frac{\hat{H}}{\hbar}\Delta t\right)^2\psi_{(x,t)}
\)

という結果が得られます。

ここで、\(\Delta t\)をN回作用させることを考えます。
すなわち、\(t^{\prime}=N\Delta t\)とおき、\(\Delta t=\frac{t^{\prime}}{N}\)として考えれば、任意の時間\(t^{\prime}\)に対して、
\(
\displaystyle \psi_{(x,t+t^{\prime})}=\left(1-i\frac{\hat{H}}{\hbar}t^{\prime}\frac{1}{N}\right)^N\psi_{(x,t)}
\)

という式が成立します。

無限回、微小区間\(\Delta t\)を動かす操作を考えればいいので、\(N\rightarrow \infty\)として考えれば、
\(
\displaystyle \psi_{(x,t+t^{\prime})}=\lim_{N\rightarrow \infty}\left(1-i\frac{\hat{H}}{\hbar}t^{\prime}\frac{1}{N}\right)^N\psi_{(x,t)}
\)

公式
\(
\displaystyle \lim_{x\rightarrow \pm\infty}\left(1+\frac{a}{x}\right)^x=e^a
\)

より、
\(
\displaystyle \psi_{(x,t+t^{\prime})}=e^{-i\frac{\hat{H}}{\hbar}t^{\prime}}\psi_{(x,t)}
\)

\(t=0\)と置き、\(t^{\prime}=t\)と文字を置き換えれば、
\(
\displaystyle \psi_{(x,t)}=e^{-i\frac{\hat{H}}{\hbar}t}\psi_{(x,0)}
\)

となり、ちゃんと時間発展演算子が導けました。めでたしめでたし。


ハミルトニアンが時間依存する場合、
\(
\displaystyle i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}_{(t)}\psi_{(x,t)}
\)

となります。無限小時間\(\Delta t\)に対して、シュレディンガー方程式は、
\(
\displaystyle i\hbar\lim_{\Delta t\rightarrow 0}\frac{\psi_{(x,t+\Delta t)}-\psi_{(x,t)}}{\Delta t}=\hat{H}_{(t)}\psi_{(x,t)}
\)

次の時刻では
\(
\displaystyle \psi_{(x,t+2\Delta t)}=
\left(1-i\frac{\hat{H}_{(t+\Delta t)}}{\hbar}\Delta t\right)\left(1-i\frac{\hat{H}_{(t)}}{\hbar}\Delta t\right)\psi_{(x,t)}
\)

となります。2乗にはなってくれません。
と、どんどん作用させていくわけです。時間依存しない時と明らかに違ったものになりますが、詳しくは次の機会のお話で。

ルンゲ=クッタ法の説明と刻み幅制御

ルンゲ=クッタ法(Runge-Kutta method、RK法)とは?
僕の知る限りの知識で紹介します。

特に良く使われる陽的ルンゲ=クッタ法は、
・実装が簡単
・良いアルゴリズムではない
という手法です。

良いアルゴリズムである陰的ルンゲ=クッタ法は、
陰的ルンゲ=クッタ法
をご覧ください。

もくじ

  1. ルンゲ=クッタ法の系統
  2. 陽的ルンゲ=クッタ法の段数と次数について
  3. 誤差について
  4. Butcher tableによるルンゲ=クッタ法の記述
  5. 埋め込まれた陽的ルンゲ=クッタ法
  6. ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)
  7. ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御のプログラム(実数/複素数 で 1/2 階微分方程式を解くプログラム)
  8. 不連続な点を含む場合
  9. 刻み幅制御のベンチマーク(振り子)
  10. (追記)ルンゲ=クッタ=ドルマンド=プリンス法
  11. 陽的ルンゲ=クッタ法の導出
  12. 参考文献

理論はいいから4次ルンゲ=クッタ法の数値計算ではどうやるかだけ教えてくれ!という方は4次ルンゲ=クッタ法へどうぞ。

ルンゲ=クッタ法の系統


ルンゲ=クッタ法は微分方程式の数値計算解を得るための手法のことを指します。

通常の参考書で微分方程式を解くために良く紹介されているのは、オイラー法と中点法、4次ルンゲ=クッタ法でしょうか。
オイラー法も中点法も4次ルンゲ=クッタ法も、”陽的ルンゲクッタ法” と呼ばれる枠組みの1つです。

オイラー法は正確には “陽的1段1次ルンゲ=クッタ法” と呼ばれ、
中点法は “陽的2段2次ルンゲ=クッタ法”、
4次ルンゲクッタ法(RK4)は、”陽的4段4次ルンゲ=クッタ法” と呼ばれています。

“段”と”次”とはなんなんでしょう?それは、
計算の大変さ(段)と、計算の正確さ(次)
です。”“の値が小さければ小さいほど計算時間が少なくて済みますし、”“の値が高ければ高いほど計算が正確です。

オイラー法は1という計算コストで正確さ1が得られますし、RK4は4という計算コストで正確さ4が得られます。

4次ルンゲ=クッタ法が使われる理由

理由は実装が簡単でそれなりの精度を持つから。です。

陽的ルンゲ=クッタ法において、pという計算コスト(p段)で、pより大きな正確さq(q次)を得ることはできません
Derivation of Runge–Kutta methodsによれば、
\(q\)次の正確さ(q次のオーダー)を得たい場合、最低限必要な段数\(p_{\mbox{min}}(q)\)は

という関係にあります。

ここで注目するべきは4次の時までは計算コストに比例して計算精度が上がっていきます。
しかし、4次以上では、計算コストの増加と計算の正確さが見合わなくなっていきます。故に
計算効率が良いのは4次だろう、
と予想できます。
また、重要な理由として、4次ルンゲ=クッタ法に現れる係数が5次以降と比べて圧倒的にシンプルであることが挙げられます。4次では\(0,1/2,1/6\)程度の係数だけが使われ、プログラム作成時の入力ミスがほぼ生まれません。しかし、5次では\(28561/56430, -7200/2197\)といった係数が数多く出てきます。
これらの理由から,4次ルンゲ=クッタ法(RK4)が数値計算科学の世界でよく使われるのです。

陽的ルンゲ=クッタ法に限って言えばプログラムの実装が非常に簡単であることが挙げられます。陰的ルンゲ=クッタ法と呼ばれるアルゴリズムもあり、これは陽的ルンゲ=クッタよりも優れていますが、計算量が多くなり、若干複雑なアルゴリズムになります。陰的ルンゲ=クッタ法を詳しく知りたい方は陰的ルンゲ=クッタ法をご覧ください。

Q. オイラー法もものすごく細かい分点を取れば、その計算効率はRK4と同じなんじゃないの?
A. 刻み幅の乗数で効いてくるのでそうではありません。高次を使っても計算が信頼できるのであれば、大きなステップで進めるほうが早いです。例え、目標精度への計算時間が同じだとしても、計算機の有効桁数によって否定されてしまいます。
RK4で典型的にとられる時間ステップの間隔は、おおよそ\(10^{-2}\sim 10^{-4}\)程度であり、RK4のエラーのオーダーは\(O(h^5)\)です。
そして、科学計算で使う際の有効桁数は倍精度型で16桁です。
1ステップ当りの誤差は\(h\)の減少に伴い、解が\(h^4\)の早さで収束していく、と言えます。
だから16桁の計算では\(h=10^{-1}\to h=10^{-4}\)に変化させる時、誤差は\(O(h^5)=c 10^{-5}\to c 10^{-20}\)と変化します。
おおよそ\(c\approx 1\)と見積もれば、(有効桁数16桁を多少超えてしまいますが、)有効桁数いっぱいまで正しい値が出るであろうと期待できます。

これに対し、オイラー法で同じことをするには\(h\)を\(10^{-16}\)にしなくてはいけません。
\(t\)の値が\(10^{-16}\)変わった時に、桁落ちの問題を回避できるほど関数\(x\)の値に変化が生まれるか?
が問題になり、まぁそんな急激な変化は生まれないでしょう、と予想できます。これでは桁落ちの問題を回避するほどの変化は到底望めません。

よって計算の効率と有効桁数の限界から、RK4なのです。

また、あまりに高次の方法を使うとルンゲ現象に代表される不安定性といったことが起こるかもしれません。
高次は高精度という意味とイコールではないことに注意しましょう。この高次≠高精度については等間隔の分点における数値積分の時に書いたので気になる方はどうぞ。。

誤差について


4次ルンゲ=クッタ法の、1ステップ当りの誤差は\(h^5\)に比例,すなわち\(O(h^5)\)です。
しかし、通常は細かいルンゲ=クッタ法を何度も繰り返して計算します。
区間\([a,b]\)を刻み幅\(h\)の4次ルンゲ=クッタ法で\(N\)回のステップを繰り返し計算する場合、誤差は
\(
\displaystyle N\cdot O(h^5)=\frac{b-a}{h}\cdot O(h^5)=O(h^4)\)

となります。よって、\(N\)回繰り返すような計算では、オーダーが1つ落ちる事に注意しましょう。

[adsense1]

陽的ルンゲ=クッタ法の段数と次数について


さて、ここまで“段”は計算コスト、で“次”は計算の正確さ、という曖昧な表現でしたが、その表現をちゃんと知りましょう。
段と次を知るためにはルンゲ=クッタ法の計算方法を知る必要があります。
具体例を載せます。
\(
\displaystyle \frac{dx}{dt}=f(t,x)
\)

の、\(t_{n}\rightarrow t_{n}+h\ \ (=t_{n+1})\)における初期値問題に対する、
陽的1段1次ルンゲ=クッタ法(オイラー法)の計算スキームは、
\(
x_{n+1}=x_{n}+h\cdot f(t_{n},x_{n})
\)
です。

陽的4段4次ルンゲ=クッタ法(RK4)の計算スキームは、
\(
\begin{align}
k_1&=f(t_n, x_n) \\
k_2&=f(t_n+h/2, x_n+h k_1/2) \\
k_3&=f(t_n+h/2, x_n+h k_2/2) \\
k_4&=f(t_n+h, x_n+h k_3) \\
x_{n+1}&=x_{n}+{(k_1+2k_2+2k_3+k_4)}h/6
\end{align}
\)
として与えられます。

一般的に、陽的s段のルンゲ=クッタ法とは
\(
\begin{align}
g_i&=x_n+h\sum_{j=1}a_{i,j}k_j\ \ \ (j\lt i, \ i=1,2,…,s) \\
k_i&=f(t_n+c_ih, g_i) \\
x_{n+1}&=x_n+h\sum_{i=1}^s b_ik_i
\end{align}
\)
として書けます。
ここで行列形式で与えられる係数\(a_{i,j}, b_{i},c_{i}\)によって、そのs段ルンゲ=クッタ法が持つ次数が決められます。段数はここから由来します。

点\((t_n, x(t_n))\)周りで関数をテーラー展開し、その関数が点\((t_n+h\ \ (=t_{n+1}), x(t_{n+1}))\)で作る点を近似解とするのがルンゲ=クッタ法です。
故に、\(x(t_{n+1})\)は、
\(
\begin{align}
x(t_{n+1})=x(t_n)+\left.\frac{h}{1!}\frac{dx}{dt}\right|_{t=t_n}+\left.\frac{h^2}{2!}\frac{d^2x}{dt^2}\right|_{t=t_n}+\left.\frac{h^3}{3!}\frac{d^3x}{dt^3}\right|_{t=t_n}+\left.\frac{h^4}{4!}\frac{d^4x}{dt^4}\right|_{t=t_n}+…
\end{align}
\)
と書けます。
ここで、テイラー展開としてどの程度一致させて\(x(t_n+h)\)を決定するか?を表すのが次数に当たります。

言葉で書くなら、

局所的な精度は,高次項が解のテイラー(Taylor)展開とどの程度マッチしているかによって測定される.オイラー法は一次精度なので,誤差は一次高く\(h^2\)というベキ乗から始まる. ~NDSolveの”ExplicitRungeKutta”メソッドより

ともあります。

Butcher tableによるルンゲ=クッタ法の記述


行列形式で与えられるルンゲ=クッタ法での係数\(a_{i,j}, b_{i},c_{i}\)は何なのか?
具体的に記述してみましょう。
オイラー法(1段1次)はもっとも単純で、係数は
\(
\begin{align}
a_{1,1}&=0 \\
b_{1}&=1 \\
c_{1}&=0
\end{align}
\)
です。これを一般的な表記法の式に当てはめれば、
\(
\begin{align}
g_1&=x_n+h a_{1,1}k_1 \\
k_1&=f(t_n+c_1h, g_1) \\
x_{n+1}&=x_n+h b_1k_1
\end{align}
\)
となります。

中点法は、
\(
\begin{align}
a_{1,1}&=0 \\
a_{1,2}&=0 \\
a_{2,1}&=1/2 \\
a_{2,2}&=0 \\
b_{1}&=0 \\
b_{2}&=1 \\
c_{1}&=0 \\
c_{2}&=1/2
\end{align}
\)
という組で与えられます。

この係数行列の組はまとめてButcher tableと呼ばれる表記をするのが便利です。

これは、\(a_{i,j}, b_{i},c_{i}\)を
Buchertable
としてまとめて書く表記法です。

再び、オイラー法はButcher tableで書くと
euler_butcher
とまとめて書くことができます。
中点法は
midpoint_butcher
RK4は
rk4_butcher
です。

高次のルンゲ=クッタ法(10,12,14次)


4次、5次…とずっとあるわけです。
こんなページがありました[3]。
High-Order Explicit Runge-Kutta Methods
この上のページには

  • 17段10次(8次が埋め込まれてる)
  • 25段12次(10次が埋め込まれてる)
  • 35段14次(12次が埋め込まれてる)

といったButcher tableにおける係数の値が書かれています。埋め込まれてる、の意味は次の節で説明します。
ただし、上のページのbutcher tableは
highorderbutcher2
となっているので注意が必要です。

埋め込まれた陽的ルンゲ=クッタ法


埋め込まれた“という表現が出てきたのでその説明を行いましょう。
日本語では『埋め込まれた陽的ルンゲ=クッタ法』、英語では『embedded explicit runge-kutta method』と呼ばれるものがあります。
これは、p段q次陽的ルンゲ=クッタ法を作ったら、別の次数の陽的ルンゲ=クッタ法も、係数行列\(a_{i,j}, c_{i}\)を使って作れるじゃありませんか!
というものです。

Butcher tableは、この場合extended Butcher tableと呼ばれ、こういう形式で書かれます。
embedded_rk
この埋め込まれたルンゲ=クッタ法のいいところは、

  1. 計算誤差の評価ができる
  2. 刻み幅を自動的に制御できる、適応刻み幅制御。(応用として。)

という点です。ルンゲ=クッタ法によって得られた解が真の解とどのくらい違っているのか?が評価できるんです。

例えば、4次のルンゲ=クッタ法を使って得られた解\(x^{(4)}(t)\)と5次のルンゲ=クッタ法を使って得られた解\(x^{(5)}(t)\)があったとします。
もしも、\(x^{(4)}(t)\)と解\(x^{(5)}(t)\)の解の差を調べ、その差が無かったらその数値計算解は真の解に限りなく近い、と判断することができ、差が大きかったらその解は真の解から離れていて、数値計算の精度が足らない、と判断することができます。どちらも1つだけの解では出来なかったことです。

精度が足らない場合、刻み幅を小さくすれば精度が上がります。また逆に、精度が十分に足りている場合、刻み幅を大きくし、計算時間を減らすことができます。
これが適応刻み幅制御なのです。

違った次数のルンゲ=クッタ法を、まるで別々に計算してもいいのですが、そうすると計算時間が単純に考えておおよそ2倍になります。
解を評価するために2倍の計算時間が必要というのは良くない計算効率です。
そこで考えられたのが埋め込まれたルンゲ=クッタ法なのです。

具体例を挙げましょう。
一番簡単な埋め込まれたルンゲ=クッタ法は、ホイン法と呼ばれています。
butcher_heun
1行目は2次のオーダーを持ち、2行目は1次のオーダーを持ちます。

また、4次と5次を持つ埋め込まれたルンゲ=クッタ法は、ルンゲ=クッタ=フェールベルグ(Runge-Kutta-Fehlberg)法と呼ばれています。
その埋め込まれたルンゲ=クッタ法は、
butcher_rkf45
と書かれます。1行目は5次のオーダー、2行目は4次のオーダーを持ちます。

ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)


さて、次数の違う2つのルンゲ=クッタ法を用いて、適応刻み幅制御を行いたいと考えます。
刻み幅を制御するにあたって、適当に精度良かったから2倍にしてもまだ大丈夫だろ、とか差が大きすぎるから刻み幅半分にしよう、ということをやってはいけません。
適当にやったら計算時間が余計にかかり、精度が良くない変な結果が得られます。

[5~9]によれば、ルンゲ=クッタ=フェールベルグ法において区間\(i\)での最適な刻み幅\(h’\)は区間\(i\)の誤差評価の結果を使って、
\(
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(5)}_{i}-x^{(4)}_{i}|}\right)^{1/4} h
\)

と予想できます。ここで\(\varepsilon\)はエラーを制御する精度の目安で、おおよそ計算時に要求する相対誤差です。もちろん、この\(h’\)は区間\(i\)の最適な幅ですが、関数に劇的な変化は無いだろうとして、次の区間の計算の一番初めに用いる刻み幅を推定するのです。
なので、\(i+1\)番目の計算区間では、計算するときはこの\(h’\)の値を使えばいいんです。
(ちなみに、m次ルンゲ=クッタ法の場合では
\(
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(m+1)}_{i}-x^{(m)}_{i}|}\right)^{1/m} h
\)

と推測されます。)

詳しい理由は分かりませんが、5次オーダーではなく、4次です。5次のオーダーは誤差評価のためだけに用いられているようです。
ルンゲ=クッタ=フェールベルグ法の計算スキームは[7]に詳しく書かれています。
日本語訳して、その計算スキームを書けば下のようになります。


ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(プログラム)


講義のレポート等の宿題で使うのは僕の意向と異なるので使用はご控えください。
研究目的、趣味、確かめの場合はミスがあるかもしれないことを念頭に置いたうえならば使用と改変をご自由にしてください。
このプログラムを使用して生じた責任は取りません。

fortran90によるプログラムです。ほぼ上の説明をそのままプログラミングしたものです。

  • 実数、1階微分方程式の場合
  • 微分方程式
    \(
    \displaystyle \frac{d}{dx}y=y\cos x,~~y(x=0)=1
    \)

    を倍精度実数、刻み幅制御で\(x=10\)まで解く事を考えます。
    解析解は
    \(
    \displaystyle y(x)=exp(sin(x))
    \)

    です。コードは以下の通りです。


  • 実数、2階微分方程式の場合
  • 微分方程式
    \(
    \displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~y(x=0)=1,~ y'(x=0)=0
    \)

    を倍精度実数、刻み幅制御で\(x=20\)まで解く事を考えます。コードは以下の通りです。

    ※ここで使われているサブルーチンdrkf45は実数一階微分のプログラム内にあるルーチンと一字一句同一です。


  • 複素数、1階微分方程式の場合
  • 微分方程式
    \(
    \displaystyle \frac{d}{dx}y=y\cos x,~~ y(x=0)=1+i\frac{1}{2}
    \)

    を倍精度実数、刻み幅制御で\(x=10\)まで解く事を考えます。コードは以下の通りです。


  • 複素数、2階微分方程式
  • 微分方程式
    \(
    \displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~ y(x=0)=1+i\frac{1}{2},~y'(x=0)=0+i0
    \)

    を倍精度実数、刻み幅制御で\(x=20\)まで解く事を考えます。コードは以下の通りです。

    ※ここで使われているサブルーチンcrkf45は複素数一階微分のプログラム内にあるルーチンと一字一句同一です。


等間隔の出力の場合は、以下の通りで実行できます。
サブルーチンはdrkf45は変わっていません。

[adsense2]

不連続な点を含む場合


不連続な点を含む場合、境界条件を指定しないと解くことはできません。

さて、ここで微分方程式
\(
\begin{eqnarray}
\frac{dy}{dx}=
\left\{
\begin{aligned}
0\;\;(x\le 0)\\
1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

を初期条件\(y(-0.095)=0\)の下で考えます(意図的に境界条件は考えません)。
これを4次ルンゲ=クッタと適応刻み幅ルンゲ=クッタで解いてみましょう。
プログラム上ではそのまま解くことが出来ます。
実際に解かせてみますと、
rk4_rkf451_disc
となり、2つの結果(赤:4次ルンゲ=クッタ、緑:刻み幅制御ルンゲ=クッタ)は異なってしまいます。これは1階微分の不連続性のため発生します。
不連続点\(x=0\)で関数\(y(x)\)に境界条件を指定しない限り、どちらも正しい解なのです。

さて、なぜこんなことが発生するのでしょうか?以下のように問題を表すことにします。

不連続な点を含む1階の微分方程式を考えます。
ここで不連続、という意味は関数\(y(x)\)の一階微分が、点\(x’\)で
\(
\displaystyle \left. \frac{dy}{dx}\right|^{x=x’+0}_{x=x’-0}=a,\;\;\;(a\ne 0)
\)

であるような点を指しているとします。

上記の例題を考えてみましょう。
上記の例題では\(a=1\)です。微分方程式を解析的に解いてみますと、
\(
\begin{eqnarray}
y(x)=
\left\{
\begin{aligned}
C_0\;\;(x\le 0)\\
x+C_1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

となります。ここで\(C_0, C_1\)は定数です。
\(C_0, C_1\)は\(y(x)\)が解きたい問題の境界条件によって決まります。

例えば、\(y(x)\)は全領域に対して繋がっている、という条件を課しましょう。この場合、不連続点\(x’\)で
\(
\displaystyle \left. y(x)\right|^{x=x’+0}_{x=x’-0}=0
\)

という境界条件を満たさなければなりません。この条件を課すと、\(C_1=C_0\)となり、初めて関数\(y(x)\)を一意に決めることが出来ます。

1階微分方程式を解く場合、適応刻み幅制御では関数\(y(x)\)は計算領域内で繋がっている事が課されています。しかし、4次ルンゲ=クッタではその条件は課されません。\(C_1\)の値は初期条件に依存し、一意に関数が決まりません。
どちらが悪いという話ではありません。

通常は適応刻み幅でも、4次ルンゲ=クッタでも\(y(x)\)にどこか連続ではない変な点がある場合、その点で区間を別々に分けて解きます。その後、境界条件に従って値を調節して全体の関数を構成します。

ベンチマーク用


微分方程式の解法がどれくらい正しそうかのベンチマーク問題として振り子(角度が大きい時)を考えましょう(振り子の詳しい解説はこちら)。
以下の\(\omega=1\)としたときの運動方程式
\(
\displaystyle \frac{d^2\theta}{dt^2}=-\sin\theta
\)


\(
t=0,\; \theta=0; \frac{d\theta}{dt}=1.9\cdots (k=0.95)
\)
の初期条件の下解いた場合、
1周期\(T\)は
\(
T=4K(0.95)=10.360044923498004876778\cdots
\)
となります。
この値はwolfram alphaから求めました。
4EllipticK[0.950.95] wolfram alpha

刻み幅制御を行い、45000周期目の値を考えます。45000周期目は時刻
\(
T_{45000}=466202.0215574102\cdots
\)
です。刻み幅制御による精度を\(10^{-12}\)に設定し、数値計算を行わせます。

すると実行結果として”fort.10″に

0.4662020113E+006  -0.2103356901E-001   0.1899883579E+001   0.4224109363E-002
0.4662020155E+006  -0.1300808922E-001   0.1899955473E+001   0.4223843994E-002
0.4662020198E+006  -0.4982881533E-002   0.1899993468E+001   0.4223658015E-002
0.4662020240E+006   0.3042061693E-002   0.1899997567E+001   0.4223520921E-002

というデータが出力されます。
1列目が時刻\(t\)、2列目が\(\theta(t)\),3列目が\(\frac{d\theta(t)}{dt}\),4列目が刻み幅\(h\)です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は77,852,488回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
4桁合っていればいい状況で8桁もあっているのは、被積分関数が変な振る舞いをしないからでしょう。

また、60000周期で確認してみると(\(T_{60000}=621602.695409880292\cdots\))

0.6216026888E+006  -0.1531918479E-001   0.1899938246E+001   0.4223959920E-002
0.6216026930E+006  -0.7293808996E-002   0.1899986003E+001   0.4223717084E-002
0.6216026973E+006   0.7312355417E-003   0.1899999862E+001   0.4223582630E-002
0.6216027015E+006   0.8756011575E-002   0.1899979827E+001   0.4223462029E-002

です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は103,803,513回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
・・・まだまだ大丈夫そうですね。

少し特殊な初期条件(\(k=1\))でやってみましょう。
この\(k=1\)は、振り子の質点がちょうど真上に来て静止する非常に不安定な状態です。
何秒間静止していられるか試してみましょう。刻み幅の制御等は上記条件と同じです。
横軸に時間\(t\)、縦軸に\(\theta(t)\)を取った時のグラフです。
ellk1
すぐに破綻しました。正しい値は\(t=20\)位までですかね。これは、機械的な誤差があることによって不安定な平衡点からほんの少しだけ外れてしまったからです。だからカオスとかを考えるときとかは注意しなければなりません。

ルンゲ=クッタ=ドルマンド=プリンス法


フェールベルグ法は初期の頃に考えられた埋め込まれた方法です。
研究が進み、実用上では別の係数の組み合わせが良いことが分かってきました。
その一つが1980年に発見されたドルマンド=プリンス(Dormand-Prince)法です。

これは、7段4(5)次の方法です。
フェールベルグ法は6段4(5)次の方法ですので、次数は同じです。

良く調べていませんが、この違いは、4次の結果を基準にして求めたブッチャー係数(フェールベルグ法)か、5次の結果を基準に求めたブッチャー係数(ドルマンド=プリンス法)か?の違いのようです。

単純に考えて、同じ次数なのにドルマンド=プリンス法の方が段数が増えていて効率が悪いです。
しかし、本来は7段なのですが、7段目に呼び出した結果を取って置けば、次のステップの1段目に同じ値が使えるように設計されているので、プログラム上は6段と(ほぼ)同じ関数の呼び出し回数になります。

プログラムはこんな感じになるでしょう。

適当な刻み幅で出力

等間隔(サブルーチンは上のものと同じなので省略)

program main
  implicit none
  integer::i,N,info,Nx
  double precision::h,tol
  double precision::xa,xb,tx
  double precision,allocatable::y(:),x(:),work0(:)
  double precision,external::grk

N=1
  allocate(y(1:N),work0(1:N))

Nx=101
  allocate(x(1:Nx))
  xa=0d0
  xb=10d0
  do i=1,Nx
     x(i)=(i-1)*(xb-xa)/dble(Nx-1)+xa
  enddo

!initial conditions
  y(1)=1d0  ! x (0)

tol=1d-8
  write(10,'(2e25.10e3)')x(1),y(1)
  do i=2,Nx
     info=-1
     h=x(i)-x(i-1)
     tx=x(i-1)
     do while(info.le.0)
        call dDP45(grk,tx,h,N,y,x(i),info,tol,work0)
     enddo
     write(10,'(2e25.10e3)')x(i),y(1)
  enddo

stop
end program main

function grk(N,x,y,s)
  implicit none
  integer,intent(in)::N,s
  double precision,intent(in)::x
  double precision,intent(in)::y(1:N)
  double precision::grk

grk=0d0
  if(s.eq.1)then
     grk=y(1)<em>cos(x)
  else
     write(6,</em>)"***Error grk"; stop
  endif

return
end function grk

4倍精度ルーチン


4倍精度のサブルーチンです。
計算速度は倍精度の30~50倍かかるので、必要なとき以外使わないようにしましょう。

陽的ルンゲ=クッタ法の導出


ルンゲ=クッタ法の導出は煩雑です。単に複雑なだけです。
導出過程について詳しく述べられているページは、早川様が記述された以下のpdfを参照すると良いと思います。
Runge-Kutta法についてのノート(早川尚男)
計算過程を含め記述されているので分かりやすいです。

参考文献


[1]Derivation of Runge–Kutta methods
[2]NDSolveの”ExplicitRungeKutta”メソッド
[3]High-Order Explicit Runge-Kutta Methods
[4]List of Runge–Kutta methods
[5]Runge-Kutta-Fehlberg Method (RKF45)
[6]Runge-Kutta-Fehlberg method
[7]Lecture:13Runge-Kutta-Fehlberg Method
[8]GPU acceleration of Runge Kutta-Fehlberg and its comparison with Dormand-Prince method
[9]William H. Pressら著『ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ– 』(1993)


↑この本が一番有益だと思います。Fortran版もありますが、英語しかありません。ちなみに、英語で,若干古いバージョンでいいのならば
公式のホームページ
Numerical Recipes in C (1992)
Numerical Recipes in Fortran 77 and Fortran 90 (1992, 1996)
で無料で公開されています。

運動量演算子の導出

量子力学の運動量演算子\(\hat{p}\)が
\(
\displaystyle \hat{p}=-i\hbar \frac{d}{d x}
\)

であらわされることを導出します。
ここでの”導出”は運動量演算子として妥当なものを導出する、と言った方がいいかもしれません。

平行移動演算子\(\hat{T}\)を出発点とします。
この平行移動演算子はある波動関数\(\psi_{(x)}\)を平行移動させる演算子で、これを作用させると
\(
\hat{T}_{(\Delta x)}\psi_{(x)}=\psi_{(x+\Delta x)}
\)

となる演算子です。

よく説明される、シュレディンガー方程式の平面波解からの出発はしません

この\(\hat{T}\)の持つであろう性質を考え、その特性から\(\hat{T}\)の具体的な形を推定し、運動量演算子を導きます。

平行移動演算子\(\hat{T}\)が満たすべき性質


当たり前と思われる性質から、平行移動演算子に関して以下の4つが言えます。

  1. 波動関数\(\psi_{(x)}\)が1に規格化済みであれば、平行移動した波動関数\(\psi_{(x+\Delta x)}\)もまた規格化されているはずである。
    すなわち、
    \(
    \begin{align}
    \langle\psi_{(x)}|\psi_{(x)}\rangle &= \langle\psi_{(x)}|\hat{T}^{\dagger}_{(\Delta x)}\hat{T}_{(\Delta x)}|\psi_{(x)}\rangle \\
    &= \langle\psi_{(x+\Delta x)}|\psi_{(x+\Delta x)}\rangle \\
    &\rightarrow \hat{T}^{\dagger}_{(\Delta x)}\hat{T}_{(\Delta x)}=\mathbf{1}
    \end{align}
    \)
    となり、3本目の式から\(\hat{T}\)はユニタリー演算子であることがわかります。
  2. \(\Delta x\)進めた後、\(\Delta x^{\prime}\)を進める操作は\(\Delta x + \Delta x^{\prime}\)進める操作に等しいはずである。
    すなわち、
    \(
    \hat{T}_{(\Delta x^{\prime})}\hat{T}_{(\Delta x)}=\hat{T}_{(\Delta x+\Delta x^{\prime})}
    \)
  3. 負の方向への平行移動\(\hat{T}_{(-\Delta x)}\)は正の方向への平行移動\(\hat{T}_{(\Delta x)}\)の逆のはずである。
    \(
    \hat{T}_{(-\Delta x)}=\hat{T}^{-1}_{(\Delta x)}
    \)
  4. \(\Delta x\)が0のとき、恒等操作なはずである。
    \(
    \begin{align}
    \hat{T}_{(0)}\psi_{(x)}&=\psi_{(x)} \\
    &\rightarrow \hat{T}_{(0)}=\mathbf{1}
    \end{align}
    \)

平行移動演算子の具体的な形


1,より、\(\hat{T}_{(\Delta x)}\)はユニタリー演算子である事が分かりました。ユニタリー演算子を作る1つの方法として、エルミート演算子\(\hat{A}_{(\Delta x)}\)を用いて、
\(
\hat{T}_{(\Delta x)}=e^{i\hat{A}_{(\Delta x)}}
\)

と置けば、\(\hat{T}_{(\Delta x)}\)はユニタリー演算子となります(エルミート演算子とユニタリー演算子の性質)。
この仮定の下、
\(
\hat{T}_{(\Delta x)}\psi_{(x)}=\psi_{(x+\Delta x)}
\)

に代入して、
\(
e^{i\hat{A}_{(\Delta x)}}\psi_{(x)}=\psi_{(x+\Delta x)}
\)

無限小移動を考えて\(\Delta_{(x)}\rightarrow 0\)と約束すれば、
\(
(1+i\hat{A}_{(\Delta x)})\psi_{(x)}=\psi_{(x+\Delta x)}
\)

変形して、
\(
\psi_{(x+\Delta x)}-\psi_{(x)}=i\hat{A}_{(\Delta x)}\psi_{(x)}
\)

辺々を\(\Delta x\)で割れば、
\(
\begin{align}
\displaystyle \lim_{\Delta x\rightarrow 0}\frac{\psi_{(x+\Delta x)}-\psi_{(x)}}{\Delta x} = \lim_{\Delta x\rightarrow 0} i\frac{\hat{A}_{(\Delta x)}}{\Delta x}\psi_{(x)}
\end{align}…(i)
\)

となります。左辺は位置\(x\)における波動関数の傾きです。
ここで右辺を見ます。\(\hat{A}_{(\Delta x)}\)は\(\Delta x\)のなんらかの関数であるはずですが、\(\hat{A}_{(\Delta x)}\)が分母の\(\Delta x\)と打ち消す形でなければ0、もしくは発散する事になります。
この場合、波動関数の傾きが到る点で0になるということは、全空間で定数ということになり非物理的な解です。また、発散してしまうのであれば波動関数が至る所で発散してしまい、これまた日物理的です。

具体的に、もしも\(\hat{A}_{(\Delta x)}=(\Delta x)^2\)だったら、
\(
\displaystyle \lim_{\Delta x\rightarrow 0} \frac{(\Delta x)^2}{\Delta x}=\lim_{\Delta x\rightarrow 0} \Delta x =0
\)

であるし、もしも\(\hat{A}_{(\Delta x)}=(\Delta x)^{1/2}\)だったら
\(
\displaystyle \lim_{\Delta x\rightarrow 0} \frac{(\Delta x)^{1/2}}{\Delta x}=\lim_{\Delta x\rightarrow 0} \frac{1}{\sqrt{\Delta x}} =\infty
\)

となります。
すなわち、言いたいことは、\(\hat{A}_{(\Delta x)}\)は\(\Delta x\)に関して1次であるべきで、
\(
\hat{A}_{(\Delta x)}=\hat{B}\cdot \Delta x
\)…(ii)
と書けるはずです。ここで\(\hat{B}\)は\(\Delta x\)に依らない演算子です(下に補足説明あり)。

よって、式(i)は\(\Delta x\rightarrow 0\)のとき、
\(
\begin{align}
\displaystyle \lim_{\Delta x\rightarrow 0}\frac{\psi_{(x+\Delta x)}-\psi_{(x)}}{\Delta x} &= \lim_{\Delta x\rightarrow 0} i\frac{\hat{A}_{(\Delta x)}}{\Delta x}\psi_{(x)} \\
\frac{d \psi_{(x)}}{dx}=i\hat{B}\psi_{(x)}
\end{align}…(iii)
\)

と書けるため、
\(
\psi_{(x)}=e^{i\hat{B}\cdot x}\psi_{(0)}
\)

となります。(この式で右辺の\(\psi_{(0)}\)が0である必要はなく、定数であればいいんです。)

平行移動演算子と運動量演算子との関係


次元について考えましょう!
今、\(\hat{T}_{(\Delta x)}=e^{i\hat{B}\cdot \Delta x}\)であり、\(\Delta x\)の次元は\(\mathrm{[m]}\)であるため、\(\hat{B}\)の次元は\(\mathrm{[m^{-1}]}\)でなければなりません。よって、\(\hat{B}\)は波数の次元を持つ演算子であると考えられるため、\(\hat{B}\)を\(c\hat{k}\)と記述することにします。ここで、\(c\)は無次元の定数です。
式(iii)より、
\(
\displaystyle \frac{d \psi_{(x)}}{dx}=ic\hat{k}\psi_{(x)}
\)

辺々に\(-i\hbar\)を掛けると、
\(
\displaystyle -i\hbar \frac{d}{dx} \psi_{(x)}=c\hbar \hat{k}\psi_{(x)}
\)
…(iv)

!!!

\(p=\hbar k\)は運動量を表すものでした。よって,\(\hat{p}=\hbar \hat{k}\)と書けば、式(iv)より
\(
\begin{align}
\hat{p}\psi_{(x)}&=\frac{1}{c}\left[-i\hbar\frac{d}{d x}\right] \psi_{(x)} \\
\end{align}
\)
となります。よって運動量演算子は
\(
\begin{align}
\hat{p}\propto -i\hbar\frac{d}{d x} \\
\end{align}
\)
となります。この考えからは定数倍を含めてしまうので、比例する、という事しか言えません。

なぜ決まらないのか考えてみますと、時間の依存性に関する考え方が式のどこにも含んでいないからです。
時間の単位が決まらないと速度の定義が出来ないわけで、これ以上進めることが出来ないのです。

各成分について言えるので、多次元の場合では
\(
\hat{p}_x\propto -i\hbar\frac{\partial}{\partial x}
\)
としても同じことです。
申し訳ないですが、数学的に厳密か?は保証しません。

補足説明


果たして運動量演算子として考えた\(e^{i\hat{B}\cdot \Delta x}\)は\(\hat{T}\)の特性1.~4.を満たすでしょうか?

  1. \(
    \left(e^{i\hat{B}\cdot \Delta x}\right)^{\dagger} \left(e^{i\hat{B}\cdot \Delta x}\right)=e^{-i\hat{B}^{\dagger}\cdot \Delta x}e^{i\hat{B}\cdot \Delta x}=\mathbf{1}
    \)
    ここで\(\hat{B}\)はエルミート演算子を考えているので、\(\hat{B}^{\dagger}=\hat{B}\)です。
    \(
    e^{-i\hat{B}^{\dagger}\cdot \Delta x}e^{i\hat{B}\cdot \Delta x}=\mathbf{1}
    \)
    ok!
  2. \(
    \hat{T}_{(\Delta x^{\prime})}\hat{T}_{(\Delta x)}=\hat{T}_{(\Delta x+\Delta x^{\prime})}
    \)
    か?
    \(
    (\mbox{左辺})=e^{i\hat{B}\cdot \Delta x^{\prime}} e^{i\hat{B}\cdot \Delta x}=e^{i\hat{B}\cdot (\Delta x^{\prime}+\Delta x)}
    \)
    ok!
  3. \(
    \hat{T}_{(-\Delta x)}=\hat{T}^{-1}_{(\Delta x)}
    \)
    か?
    両辺の左から\(\hat{T}_(\Delta x)\)を掛けて、
    \(
    \begin{align}
    &\hat{T}_{(\Delta x)}\hat{T}_{(-\Delta x)}=\mathbf{1} \\
    &= e^{i\hat{B}\cdot \Delta x}e^{i\hat{B}\cdot (-\Delta x)}=e^{i\hat{B}\cdot 0}=\mathbf{1}
    \end{align}
    \)

    ok!
  4. \(e^{i\hat{B}\cdot 0}=1\)
    ok!

よって、4つすべてを満たすので、\(\hat{T}_{(\Delta x)}=e^{i\hat{B}\cdot \Delta x}\)は平行移動演算子として適当であると考えることができるのです。

数値積分(等間隔)

数値計算での積分方法、特に等間隔の分点の場合であるニュートン・コーツ積分(とロンバーグ積分)に関する理論とのプログラムを載せます。

  1. ニュートンコーツ型の公式
  2. ルンゲ現象
  3. ロンバーグ積分
  4. サブルーチン”integral”のコードと使い方
  5. サブルーチン”romberg”のコードと使い方

等間隔の分点で積分を行う場合、よく使われる方法は台形則、もしくはシンプソン積分です。

台形則

台形則は、分点間を台形近似して面積を求める方法であり、下図のようなイメージで積分を近似する方法です。
外形_台形_c
積分\(\displaystyle \int_a^b f(x) dx\)を、
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N-1} \frac{f(x_{i+1})+f(x_i)}{2}h
\)

として近似します。この時、誤差のオーダーは\(O(h^3)\)となります。
言い換えると、
台形積分とは、隣り合う分点間を1次関数で近似して求める積分
言えます。

シンプソン積分

シンプソン積分は、隣り合う分点間を2次関数で近似して求める積分です。イメージは
外形_シンプソン_c
な感じです。数式では
\(
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N/2} \frac{f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})}{3}h
\)

となります。

高次≠高精度

これ以降もあります。分点間を3次、4次、…としていけば高次の積分をすることが可能になります。
この分点間の次数をどんどん上げていって積分をする方法をまとめてニュートンコーツの公式、と呼びます。
例えばの台形則はニュートンコーツの公式の1次に相当し、シンプソン則はニュートンコーツの公式の2次に相当しています。

じゃあ100次の公式作れば精度が凄い公式ができるんじゃないか?となりますが、これは違います。
高次であることと精度が高いことは違うのです。だから使うな、というわけではなくて知っておいてほしいことなんです。
ではなぜか?これはwikipediaのニュートンコーツの公式のページからとってきたものですが、この式に理由があります。
ニュートンコーツ型式_c

問題となるのが誤差項です。例えば台形近似を見ますと、誤差項は\(h^3 f^{(2)}(\xi)\)とあります。hは刻み幅なので、この項は刻み幅を小さくすれば小さくするほど精度がよくなることを言っています。
が、問題は\(f^{(2)}(\xi)\)の方です。この項は本当の関数が分からないと評価できません。
多項式だったら微分するごとに次数が1つづつ減少するのでこの項は高次になれば小さくなります。じゃあもしも微分するごとに値が増加していったら…?この場合、hを小さくしてもそれを上回るくらい\(f^{(2)}(\xi)\)が大きくなることがあります。これをルンゲ現象と呼びます。

ルンゲ現象

例えば関数
\(
\displaystyle f(x)=\frac{1}{1+25x^2}
\)

を考えます。wikipediaのページを参考にすれば、その導関数と値は、
ルンゲ現象_c
となります。高次になればなるほど値が増加する関数なのです。
実際には\(h^n f^{(m)}(\xi)\)の積の値が小さいか大きいか?なので、一概には高次は悪い!とは言えません。高次は危ないくらいでしょう。

これを可能な限り回避するために、低次の積分を組み合わせて使う方法がとられます。
こうすることで高次の値を増加を抑えるのです。

そのほかの手法もいろいろあります。分点の位置に縛られなければガウス求積法が最高でしょう。素晴らしい方法です。ぜひ調べてみてください。

[adsense1]

ロンバーグ積分

さて、高次を追い求めるのもいいですが、一線を駕す…とまではいきませんが、ロンバーグ(Romberg)積分というものがあります。
このロンバーグ積分とは、台形則による計算を基本とし、無限小区間での積分結果を補外によって求めよう、というものです。
ここでいう補外の概念は以下のような感じです。

刻み幅1で台形則の結果から積分値S1が分かる

刻み幅0.5で台形則による結果から積分値S2がわかる

刻み幅0.25で台形則による結果から積分値S3がわかる

刻み幅0.175で台形則による結果から積分値S4がわかる

とある程度計算していきます。そうすると刻み幅hを変化させるとそれに応じて積分値がS1→S2→S3→S4…と推移していくわけです。
この推移の仕方を計算すれば、刻み幅無限小の場合の積分値Sが分かる、こういう仕組みになっています。収束加速法と言ったり、リチャードソン補外、なんて言ったりします。
この補外による方法はかなり優れています。先ほどの高次がどうこう、といった問題が発生しないのです。台形則しか使ってないんですから。有限桁計算に適した方法である、とかwikipediaに載ってたりします。

ロンバーグ積分の詳しい計算方法や理論は
ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ


を見ると詳しく書かれています。

日本語では第6章 数値積分と数値微分にロンバーグ積分の記述があります。
英語ではRomberg integralで調べると良いでしょう。参考までに、英語では
Romberg integration exampleや、Romberg Integrationなにかがいいと思います。

2016/03/01
上のRomberg Integrationを参考にして書いたプログラムはこちらです。これは、解析関数fをRomberg積分するプログラムです。
解析関数ではなく、サイズ\(2^{n}+1\)の配列に格納されている場合はこのページの下の方にあるプログラムを参照してください。

call romberg(a,b,s,pre)

で呼び出し、区間a~bの積分\(s=\int_a^b f(x) dx\)を精度preで求めるプログラムです。

サブルーチン”integral”のコードと使い方

下のサブルーチン”integral”は1,2,3次元の短冊近似、台形則、シンプソン積分、シンプソン3/8積分、ブール則(ニュートン・コーツ型積分)をカバーします。
ただし、分点の個数に注意してください。分点の個数を気にせず使える近似方法は、短冊近似、台形近似のみです。

ロンバーグ積分では引数の関係上、ニュートン・コーツ型積分との間をうまく処理できなかったため、別のサブルーチンとして書きます。もしもロンバーグ積分を使いたい場合はもう少し先に進んでください。

まずはニュートン・コーツ型積分です。
下のサブルーチンを使ってください。
使い方は、

 call integral(配列の大きさ,配列,刻み幅,積分結果,積分手法)

という感じです。
実際のプログラムでは、以下のように使用してください。配列yは複素数解列(ただし実軸上の積分)でもokです。

1次元

サブルーチンintegralの使い方
積分則 呼び方
短冊近似(長方形近似)
 call integral(size(w,1),w,h,s,"box")
配列yの大きさに指定はない。
台形近似
 call integral(size(y,1),y,h,s,"trapezoid")
配列yの大きさに指定はない。
シンプソン則
 call integral(size(y,1),y,h,s,"simpson")
配列yの大きさが3,5,7,…個、\(2n+1\)個でないといけない
シンプソン3/8則
 call integral(size(y,1),y,h,s,"simpson38")
配列yの大きさが4,7,10,…個、\(3n+1\)個でないといけない
ブール則
 call integral(size(y,1),y,h,s,"boole")
配列yの大きさが5,9,13,…個、\(4n+1\)個でないといけない

2次元配列の場合は以下のように指定してください。

 call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")

3次元の場合はこう指定してください。

 call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")

積分のためのモジュール

module integral_mod
  !developer --> sikino
  !date --> 2015/04/07
  implicit none
  interface integral
     module procedure &
          dintegral, &
          dintegral2d, &
          dintegral3d, &
          cintegral, &
          cintegral2d, &
          cintegral3d
  end interface integral
contains
  subroutine dintegral(N,y,h,s,method)
    integer,intent(in)::N
    double precision,intent(in)::h,y(1:N)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::y0,y1,y2,y3
       
    s=0.d0; y0=0.d0; y1=0.d0; y2=0.d0; y3=0.d0
   
    if(trim(method).eq."box")then
       s=h*sum(y(1:N-1))
    elseif(trim(method).eq."trapezoid")then
       y1=y(1)+y(N)
       do i=2,N-1
          y2=y2+y(i)
       enddo
       s=(y1+2.d0*y2)*h*0.5d0
    elseif(trim(method).eq."simpson")then
       if(mod(N,2).ne.1)then
          write(6,*)"=====cannot calculation with simpson"
          write(6,*)"=====program stop"
          stop
       endif

       y1=y(1)+y(N)
       do i=2,N-1,2
          y2=y2+y(i)
       enddo
       do i=3,N-2,2
          y3=y3+y(i)
       enddo

       s=(y1+4.d0*y2+2.d0*y3)*h/3.d0

    elseif(trim(method).eq."simpson38")then

       if(mod(N,3).ne.1)then
          write(6,*)"=====cannot calculation with simpson38"
          write(6,*)"=====program stop"
          stop
       endif

       y0=y(1)+y(N)
       do i=2,N-2,3
          y1=y1+y(i)
       enddo
       do i=3,N-1,3
          y2=y2+y(i)
       enddo
       do i=4,N-3,3
          y3=y3+y(i)
       enddo
       s=(y0+3.d0*(y1+y2)+2.d0*y3)*3.d0*h/8.d0        

    elseif(trim(method).eq."boole")then
       if(mod(N,4).ne.1)then
          write(6,*)"=====cannot calculation with boole"
          write(6,*)"=====program stop"
          stop
       endif

       y0=y(1)+y(N)
       do i=5,N-4,4
          y1=y1+y(i)
       enddo
       do i=2,N-1,2
          y2=y2+y(i)
       enddo
       do i=3,N-2,4
          y3=y3+y(i)
       enddo

       s=(14.d0*y0+28.d0*y1+64.d0*y2+24.d0*y3)*h/45.d0
    else
       write(6,*)"=====cannot calculation in integral"
       write(6,*)"=====program stop"
       stop
    end if

    return
  end subroutine dintegral
 
  subroutine dintegral2d(Nx,Ny,z,hx,hy,s,method)
    implicit none
    integer,intent(in)::Nx,Ny
    double precision,intent(in)::hx,hy,z(1:Nx,1:Ny)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::ty(1:Ny),r(1:Nx)

    s=0.d0
    ty(1:Ny)=0.d0
    r(1:Nx)=0.d0
    do i=1,Nx
       ty(1:Ny)=z(i,1:Ny)
       call integral(Ny,ty,hy,s,method)
       r(i)=s
    enddo
    call integral(Nx,r,hx,s,method)
       
    return
  end subroutine dintegral2d
 
  subroutine dintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
    implicit none
    integer,intent(in)::Nx,Ny,Nz
    double precision,intent(in)::hx,hy,hz,w(1:Nx,1:Ny,1:Nz)
    character(*),intent(in)::method
    double precision,intent(out)::s
    integer::i
    double precision::tyz(1:Ny,1:Nz),r(1:Nx)
       
    s=0.d0
    tyz(1:Ny,1:Nz)=0.d0
    r(1:Nx)=0.d0
    do i=1,Nx
       tyz(1:Ny,1:Nz)=w(i,1:Ny,1:Nz)
       call integral(Ny,Nz,tyz,hy,hz,s,method)
       r(i)=s
    enddo
    call integral(Nx,r,hx,s,method)    
       
    return
  end subroutine dintegral3d
 
  subroutine cintegral(N,y,h,s,method)
    integer,intent(in)::N
    complex(kind(0d0)),intent(in)::y(1:N)
    double precision,intent(in)::h
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(N,dble(y),h,res,trim(method))
    call integral(N,dimag(y),h,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral

  subroutine cintegral2d(Nx,Ny,z,hx,hy,s,method)
    integer,intent(in)::Nx,Ny
    complex(kind(0d0)),intent(in)::z(1:Nx,1:Ny)
    double precision,intent(in)::hx,hy
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(Nx,Ny,dble(z),hx,hy,res,trim(method))
    call integral(Nx,Ny,dimag(z),hx,hy,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral2d

  subroutine cintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
    integer,intent(in)::Nx,Ny,Nz
    complex(kind(0d0)),intent(in)::w(1:Nx,1:Ny,1:Nz)
    double precision,intent(in)::hx,hy,hz
    character(*),intent(in)::method
    complex(kind(0d0)),intent(out)::s

    double precision::res,ims

    s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

    call integral(Nx,Ny,Nz,dble(w),hx,hy,hz,res,trim(method))
    call integral(Nx,Ny,Nz,dimag(w),hx,hy,hz,ims,trim(method))
   
    s=dcmplx(res,ims)
   
    return
  end subroutine cintegral3d
end module integral_mod

サブルーチン”integral”を用いた例題

実際の使い方。参考にどうぞ。
下のコードは2次元の積分
\(
\displaystyle \int\int xe^{-x^2}e^{-y}dxdy=\frac{e^{-x^2}}{2}e^{-y}+\mbox{const}
\)

を数値積分します。積分範囲を\(x=0\sim 3, \; y=0 \sim 5\)にした場合、その解析解は
\(
\begin{align}
\int_0^3 dx\int_0^5 dy xe^{-x^2}e^{-y}&=\frac{1}{2}(1-e^{-5}-e^{-9}+e^{-14})\\
&=0.4965697373627734784608751894356320535936864348993604\cdots
\end{align}
\)

となります。実行すると

./a.out
simpson   :    0.496569737366759E+000
romberg,3 :    0.496569737362773E+000

となります。コードは、

program main
  use integral_mod
  implicit none
  integer::i,j,n,Nx,Ny
  double precision::hx,hy,ans
  double precision::xmin,xmax,ymin,ymax
  double precision,allocatable::x(:),y(:)
  double precision::s
  double precision,allocatable::w(:,:)
  double precision::f
  external::f
 
  n=10

  xmin=0d0
  xmax=3d0

  ymin=0d0
  ymax=5d0

  Nx=2**(n)
  Ny=2**(n)
  allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
  x=0d0; y=0d0; w=dcmplx(0d0,0d0)

  hx=(xmax-xmin)/dble(Nx)
  hy=(ymax-ymin)/dble(Ny)
 
  do i=0,Nx
     x(i)=xmin+hx*dble(i)
  enddo
  do i=0,Ny
     y(i)=ymin+hy*dble(i)
  enddo

  do i=0,Nx
     do j=0,Ny
        w(i,j)=f(x(i),y(j))
     end do
  end do

  call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")
  write(6,'(A,e25.15e3)')"simpson   : ",s
 
  call integral(size(w,1),size(w,2),w,hx,hy,s,"romberg",3)
  write(6,'(A,e25.15e3)')"romberg,3 : ",s
 
  return
end program main

function f(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::f
 
  f=x*exp(-x*x)*exp(-y)
 
  return
end function f

3次元の計算の例題です。
下のコードは
\(
\displaystyle \int\int\int x^4y^4z^4 dxdydz=\left(\frac{1}{5}\{0.7^5-(-1)^5\}\right)^3
\)

を数値積分します。積分範囲を\(x,y,z=-1\sim 0.7\)にした場合、その値は
\(
\displaystyle \int_{-1}^{0.7}\int_{-1}^{0.7}\int_{-1}^{0.7} x^4y^4z^4 dxdydz = 0.127496010896795\cdots
\)

となります。実行すると、

$ gfortran integralmod.f90 main.f90
$ ./a.out
analysis  :   0.127496010896795E-01
simpson   :   0.10507613E-006  0.10507613E-006
romberg,3 :  -0.17347235E-017 -0.17347235E-017

と得られます。中身はこうです。分点の数はx,y,z軸各々で違っていて構いません。

program main
  use integral_mod
  implicit none
  integer::i,j,k,n,Nx,Ny,Nz
  double precision::hx,hy,hz,ans
  double precision,allocatable::x(:),y(:),z(:)
  complex(kind(0d0))::s
  complex(kind(0d0)),allocatable::w(:,:,:)
 
  n=6
  Nx=2**(n)
  Ny=2**(n-1)
  Nz=2**(n+1)
 
  allocate(x(0:Nx),y(0:Ny),z(0:Nz),w(0:Nx,0:Ny,0:Nz))
  x=0d0; y=0d0; z=0d0; w=dcmplx(0d0,0d0)
 
  hx=1.7d0/dble(size(w,1)-1)
  hy=1.7d0/dble(size(w,2)-1)
  hz=1.7d0/dble(size(w,3)-1)
 
  do i=0,Nx
     x(i)=hx*dble(i)-1d0
  enddo
  do i=0,Ny
     y(i)=hy*dble(i)-1d0
  enddo
  do i=0,Nz
     z(i)=hz*dble(i)-1d0
  enddo

  do i=0,Nx
     do j=0,Ny
        do k=0,Nz
           w(i,j,k)=dcmplx((x(i)*y(j)*z(k))**4,(x(i)*y(j)*z(k))**4)
        end do
     end do
  end do

  ans=(0.7d0**5.d0+1.d0)/5.d0
  ans=ans**3.d0
  write(6,'(A,e23.15e2)')"analysis  : ",ans
 
  call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")
  write(6,'(A,2e17.8e3)')"simpson   : ",dble(s)-ans,dimag(s)-ans
 
  call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"romberg",3)
  write(6,'(A,2e17.8e3)')"romberg,3 : ",dble(s)-ans,dimag(s)-ans
 
  return
end program main

ロンバーグ積分のプログラム


ロンバーグ積分
 call romberg(jx,x,y,s,pre)
配列yの大きさは2,3,5,10,…,\(2^{jx}+1\)個でないといけない ロンバーグ積分では収束精度”pre”を指定すること。もしも収束精度に達しなかった場合、警告と共に、与えられた値での収束限界の積分結果を返す。積分精度を高めたければ刻み幅を小さくする操作(例えば分点数を増やす等)をしてください。”pre”は、積分結果の絶対値が1より大きくなる場合は相対誤差を、小さくなる場合は絶対誤差を取ります。この理由は以前の結果の補正を加え、収束をさせるためであり、非常に小さい誤差の場合いつまでたっても機械誤差のため収束しないからです。
module romberg_mod
  !developer --> sikino
  !date --> 2016/03/01
  !         2016/03/03
  !
  ! 1D case :
  !  romberg(jx,x,y,s,pre)
  !          |  | | | +- (in)  precision (e.g. 1d-8)
  !          |  | | +--- (out) integration result
  !          |  | +----- (in)  y(1:2**jx+1) f(x)
  !          |  +------- (in)  x(1:2**jx+1) x
  !          +---------- (in)  related to array size
  !
  ! 2D case :
  !  romberg(jx,jy,x,y,z,s,pre)
  !          |  |  | | | | +-- (in)  precision (e.g. 1d-8)
  !          |  |  | | | +---- (out) integration result
  !          |  |  | | +------ (in)  z(1:2**jx+1,1:2**jy+1) f(x,y)
  !          |  |  | +-------- (in)  y(1:2**jy+1) y
  !          |  |  +---------- (in)  x(1:2**jx+1) x
  !          |  +------------- (in)  related to array size y
  !          +---------------- (in)  related to array size x
  !
  ! 3D case :
  !  romberg(jx,jy,jz,x,y,z,w,s,pre)
  !          |  |  |  | | | | | +-- (in)  precision (e.g. 1d-8)
  !          |  |  |  | | | | +---- (out) integration result
  !          |  |  |  | | | +------ (in)  w(1:2**jx+1,1:2**jy+1,1:2**jz+1) f(x,y,z)
  !          |  |  |  | | +-------- (in)  z(1:2**jz+1) z
  !          |  |  |  | +---------- (in)  y(1:2**jy+1) y
  !          |  |  |  +------------ (in)  x(1:2**jx+1) x
  !          |  |  +--------------- (in)  related to array size z
  !          |  +------------------ (in)  related to array size y
  !          +--------------------- (in)  related to array size x
  !
  implicit none
  interface romberg
     module procedure &
          dromberg, &
          dromberg2d, &
          dromberg3d
  end interface romberg
contains
  subroutine dromberg(jx,x,y,s,pre)
    implicit none
    integer,intent(in)::jx
    double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
    double precision,intent(out)::s

    integer,parameter::jm=6 !--> precision: O(h^(2*jm))
    integer,parameter::nm=2**jm
    double precision,allocatable::tx(:),ty(:)
    double precision::ts
    integer::k

    s=0d0
    if(jx.ge.jm)then
       allocate(tx(1:nm+1),ty(1:nm+1)); tx=0d0; ty=0d0
       do k=0,2**(jx-jm)-1
          ts=0d0
          tx(1:nm+1)=x(k*nm+1:(k+1)*nm+1)
          ty(1:nm+1)=y(k*nm+1:(k+1)*nm+1)

          call romberg_sub(jm,tx,ty,ts,pre)
          s=s+ts
       enddo
       deallocate(tx,ty)    
    else
       call romberg_sub(jx,x,y,s,pre)
    endif

    return
  end subroutine dromberg

  subroutine romberg_sub(jx,x,y,s,pre)
    ! reference "http://www.math.usm.edu/lambers/mat460/fall09/lecture29.pdf"
    implicit none
    integer,intent(in)::jx
    double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
    double precision,intent(out)::s

    integer::i,j,k,n,dn
    double precision::h,ps,tmp
    double precision::T(1:jx+1,1:jx+1)

    n=2**jx+1

    h=x(n)-x(1)
    dn=(n-1)/2

    T(1,1)=0.5d0*h*(y(1)+y(n))
    s=T(1,1)
    ps=s
    h=0.5d0*h
    do j=2,jx+1

       ! trapezoidal rule
       tmp=0d0
       do i=1,2**(j-2)
          tmp=tmp+y(1+(2*i-1)*(dn))
       enddo
       T(j,1)=0.5d0*T(j-1,1)+h*tmp

       do k=2,j
          ! Richardson extrapolation
          T(j,k)=T(j,k-1)+(T(j,k-1)-T(j-1,k-1))/(dble(4**(k-1))-1d0)
       enddo
       s=T(j,j)

       ! precision check
       if(abs(s).ge.1d0)then
          if(abs((ps-s)/s).le.pre)exit
       else
          if(abs((ps-s)).le.pre)exit
       endif
       ps=s
       h=0.5d0*h
       dn=dn/2
    enddo

    if(j-1.eq.jx)then
       write(6,'(A)')" -+-+- didn't converge at romberg integral -+-+- "
       write(6,'(A)')"       Please change stepsize h of array x  "
    endif

    return
  end subroutine romberg_sub

  subroutine dromberg2d(jx,jy,x,y,z,s,pre)
    implicit none
    integer,intent(in)::jx,jy
    double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1)
    double precision,intent(in)::z(1:2**jx+1,1:2**jy+1)
    double precision,intent(out)::s

    integer::i,nx,ny
    double precision::ty(1:2**jy+1),r(1:2**jx+1)
    nx=2**jx+1
    ny=2**jy+1

    s=0.d0
    ty(1:ny)=0.d0
    r(1:nx)=0.d0
    do i=1,nx
       ty(1:ny)=z(i,1:ny)
       call romberg(jy,y,ty,s,pre)
       r(i)=s
    enddo
    call romberg(jx,x,r,s,pre)
       
    return
  end subroutine dromberg2d
 
  subroutine dromberg3d(jx,jy,jz,x,y,z,w,s,pre)
    implicit none
    integer,intent(in)::jx,jy,jz
    double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1),z(1:2**jz+1)
    double precision,intent(in)::w(1:2**jx+1,1:2**jy+1,1:2**jz+1)
    double precision,intent(out)::s
    integer::i,nx,ny,nz
    double precision::tyz(1:2**jy+1,1:2**jz+1),r(1:2**jx+1)
       
    nx=2**jx+1; ny=2**jy+1; nz=2**jz+1
    s=0.d0
    tyz(1:ny,1:nz)=0.d0
    r(1:nx)=0.d0
    do i=1,nx
       tyz(1:ny,1:nz)=w(i,1:ny,1:nz)
       call romberg(jy,jz,y,z,tyz,s,pre)
       r(i)=s
    enddo
    call romberg(jx,x,r,s,pre)  
       
    return
  end subroutine dromberg3d
end module romberg_mod

ロンバーグ積分(romberg_mod)の例


必要なのは、上で紹介したモジュール “romberg_mod” と下のメインプログラムです。
以下のプログラムは1次元の定積分
\(
\int_1^10 \frac{1}{x^2} dx = 0.9
\)
を分点数\(2^8+1(jx=8)\)個でロンバーグ積分するものです。
精度は\(O(h^{2\cdot min{jm,jx}})\)です。
ここで\(jm\)はモジュールromberg_modの中にパラメータとして宣言されています。

コンパイルは

ifort romberg_mod.f90 main.f90

とでもすればいいでしょう。

program main
  use romberg_mod
  implicit none
  integer::i,jx,n
  double precision::h,a,b,s
  double precision,allocatable::x(:),y(:)
  double precision::f
  external::f
 
  a=1d0
  b=10d0
  jx=8
  n=2**jx+1
  allocate(x(1:n),y(1:n))
  h=(b-a)/dble(n-1)
  do i=1,n
     x(i)=a+h*dble(i-1)
     y(i)=f(x(i))
  enddo
 
  call romberg(jx,x,y,s,1d-8)
  write(6,*)s,"romberg"
 
  stop
end program main

function f(x)
  implicit none
  double precision::f,x
 
  f=1d0/(x*x)
 
  return
end function f

実行例(要求精度は\(10^{-8}\))

>./a.out
  0.900000000062642  
>

2次元ではこう。

program main
  use romberg_mod
  implicit none
  integer::i,j,Nx,Ny,jx,jy
  double precision::hx,hy
  double precision::xmin,xmax,ymin,ymax
  double precision,allocatable::x(:),y(:)
  double precision::s
  double precision,allocatable::w(:,:)
  double precision::f
  external::f
 
  jx=10
  jy=8
  xmin=0d0
  xmax=3d0

  ymin=0d0
  ymax=5d0

  Nx=2**(jx)
  Ny=2**(jy)
  allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
  x=0d0; y=0d0; w=dcmplx(0d0,0d0)

  hx=(xmax-xmin)/dble(Nx)
  hy=(ymax-ymin)/dble(Ny)
 
  do i=0,Nx
     x(i)=xmin+hx*dble(i)
  enddo
  do i=0,Ny
     y(i)=ymin+hy*dble(i)
  enddo

  do i=0,Nx
     do j=0,Ny
        w(i,j)=f(x(i),y(j))
     end do
  end do

  call romberg(jx,jy,x,y,w,s,1d-8)
  write(6,*)s
 
  stop
end program main

function f(x,y)
  implicit none
  double precision,intent(in)::x,y
  double precision::f
 
  f=x*exp(-x*x)*exp(-y)
 
  return
end function f

実行例

>./a.out
  0.496569737374163
>

[adsense2]

HTMLで行間をちょっとだけ開ける?

HTMLで行間を半行開けたい!
そう思ったのが始まりでした。が!
無理らしいです。CSSでやらないといけないみたいです。
まず、通常の改行<br />によって開けると、
———————–
AAA
AAA

AAA

AAA
AAA
———————–
のようにかなりの幅が開いてしまいます。これを
———————–
AAA
AAA
AAA
AAA
AAA
———————–
のように開けます。

ここでやっているのは半行開けるのではなく、行に2倍の広さを持たせて、その中央に文字を書く、ということを行っているのです。

すなわち、
<br /> AAA <br />
によって改行をあけるのではなく、
<span style=”line-height:200%;”>AAA</span>
と書くことです。

イメージとしては

な感じです。