「プログラミングと数値計算」カテゴリーアーカイブ

-1のn乗の計算速度

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

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

です。


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

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


\(
(-1)^n
\)


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


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


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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2次元の場合


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

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

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

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

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

となります。

基底関数の対応


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

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

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

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

プログラム


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

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

ifort -mkl integrate.f90 main.f90

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

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

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

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

具体的な計算例


三角形型のポテンシャル

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

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

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

1次元時間依存シュレーディンガー方程式を解かねばならぬ

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、1コア使用)で、\(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)系の固有状態を与える場合の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)\)について求めていけば、

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

数値計算プログラム


プログラムは以下のリンクにあります。
http://slpr.sakura.ne.jp/qp/supplement_data/highprecision_1Dtdsetise.tar.gz
そのままのプログラムでは、以下の設定となっています。
プログラムの利用に関してはhttp://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/連立線形微分方程式等にあるのでそちらを参照してください)。

gnuplot上で行う任意関数のフィッティング

gnuplotには、”fit”というコマンドがあります。
これを利用すると、データ列に対して、任意の自由度を含んだ関数\(f(x)\)でフィッティングを行うことが出来ます。
手法はMarquardt-Levenberg法に基づいた非線形最小二乗法によるフィッティングです。

gnuplot ver4.6の場合

データ列”data.d”を\(f(x)=ax^2+bx+c\)でフィッティングする\(a,b,c\)を求めたい場合、gnuplot上で

f(x)=ax**2+bx+c
fit f(x) "data.d" u 1:2 via a,b,c

※gnuplot ver5.0以降の場合、エラーバーで重みづけしたフィッティングが出来るようです。
http://gnuplot.sourceforge.net/demo_5.0/fit.html

説明


以下のデータ列(ファイル名”data.d”)に対してgnuplot上でフィッティングを行います。

データは、関数
\(
\displaystyle f(x)=\frac{1}{1+e^{x-15}}
\)

に乱数を与えて生成したデータです。

このデータ列に対して、
\(
\displaystyle f(x)=\frac{a}{1+e^{(x-b)/c}}
\)

でフィッティングを行い、未知の定数\(a,b,c\)を決めたいと思います。

gnuplot上で

f(x)=a/(1+exp((x-b)/c))
a=1e0
b=10e0
c=0.5e0
fit f(x) "data.d" u 1:2 via a,b,c

と打つと、\(a,b,c\)にフィッティングした後の定数が入ります。
ここで、フィッティング前に\(a,b,c\)の値を入れているのは初期値です。
初期値があまりにも違うとうまくフィッティングが失敗するので、出来るだけ近い値を入れましょう。

実際に動かすと

> fit f(x) "data.d" u 1:2 via a,b,c
...

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 0.995618         +/- 0.005679     (0.5704%)
b               = 15.3268          +/- 0.051        (0.3327%)
c               = 0.970465         +/- 0.04439      (4.574%)

...
>

という文章が出力されます。

フィッティングの結果、
\(a=0.995618,b=15.3268,c=0.970465\)と求められました。
乱数を与える前のデータは\(a=1,b=15,c=1\)ですので、元の関数に近いことが分かります。
fitコマンドを動かした後は既に\(a,b,c\)に値が代入されているので、そのまま

plot f(x)

とすれば、フィッティング結果をプロットすることが出来ます。

データと共に載せれば以下のようになります。

フィッティングする範囲を制限したければ、

fit [10:20] f(x) "data.d" u 1:2 via a,b,c

とすると、\(x\)の範囲を\([10,20]\)に制限することが出来ます。

2変数関数のフィッティング


2変数関数であろうとフィッティングは可能です。
元のデータ
(http://slpr.sakura.ne.jp/qp/supplement_data/data2.d)
は、関数
\(
g(x,y)=x\sin(xy+0.5)
\)

に適当な乱数を足して作られたデータです。グラフにすれば

となります。

このデータに対して、gnuplot上で関数
\(
g(x,y)=ax\sin(bxy+c)
\)

によるフィッティングを行います。コマンドは

g(x,y)=a*x*sin(b*x*y+c)
a=1.3e0
b=0.7e0
c=0.3e0
fit g(x,y) "data2.d" u 1:2:3 via a,b,c

です。
実行すると

>fit g(x,y) "data2.d" u 1:2:3 via a,b,c

...
Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 1.16313          +/- 0.004072     (0.3501%)
b               = 1.00073          +/- 0.001672     (0.1671%)
c               = 0.495134         +/- 0.003432     (0.6932%)
...

>

と得られますので、フィッティング結果と共にプロットすれば

となります。

初期値が答えに割と近いですが、これは失敗を防ぐためです。全ての初期値のパラメータを0で行った所失敗しました。

参考サイト


Gnuplotでの関数fitting
gnuplot demo script: fit.dem -ver5.0
gnuplot demo script: fit.dem -ver4.6
Fit (manual)
gnuplot で最小二乗フィッティングする -ゴルディアスの涙目

ゼロ点を探す(二分法、挟み撃ち法、Brent法、Newton法、Steffensen法)

ここでは、非線形方程式の解を数値的に求める方法である、
二分法、
挟み撃ち(False position)法、
Brent法、
Newton法、
Steffensen法
について紹介します。

  1. 問題設定
  2. まとめ
  3. 数値解法
  4. 複数のゼロ点を見つけたい場合

問題設定


区間\([a,b]\)で定義された実関数\(f(x)\)に
\(
f(x)=0
\)

を満たす\(x\)が1つある。\(x\)を求めよ。但し\(f(a),f(b)\ne 0\)。

まとめ


1変数の場合:Brent法(二分法と逆2次補間、挟み撃ち法の組み合わせ)
多変数の場合:Newton法
を使うのが良いでしょう。

数値解法


1変数の問題で、
教育上良く使われる方法は二分法
実用的な手法はBrent法
です。
またNewton法は多変数の場合にも複素数の場合でも容易に拡張できるので、その意味で実用的です。

区間\([a,b]\)にゼロ点が1つ存在する事を数学的に言えば、\(f(a)f(b)\lt 0\)ということです。

二分法


要点

・解が確実に見つかります。
・関数の振る舞いが非連続であっても、たちが悪い関数でも確実ですが、収束は遅いです。

計算方法

二分法は位置\([a,b]\)における関数の符号の情報だけを見てゼロ点を推測します。
なので、区間内にゼロ点がある事が分かっていれば、そこを境として符号が違うので確実に解を囲い込んでいくことが出来ます。
教育的に分かりやすい方法であり、根を探すにあたり失敗が無い方法ですが、
滑らかな関数であっても計算時間がかかるため、あまり実用的ではありません。

収束速度について

2分法では、1回の計算あたり範囲が\(1/2\)になるため、漸化式
\(
\displaystyle \varepsilon_n=\varepsilon_{n-1}\cdot 2^{-1}
\)

が成立するため、
\(
\varepsilon_n=\varepsilon_{0}\cdot 2^{-n}
\)

が導けます。ここで\(\varepsilon_{0}\)は初期のxの範囲(\(\varepsilon_0=|b-a|\))を表します。
収束精度と計算回数の関係を決めるためにはnについて解いて、
\(
\begin{align}
\varepsilon_n &=2^{-n} \varepsilon_{0} \\
2^n &=\frac{\varepsilon_0}{\varepsilon_n} \\
n &=\log_2\left(\frac{\varepsilon_0}{\varepsilon_n}\right) \\
\mbox{または}& \\
n &=\frac{\ln\left(\varepsilon_0/\varepsilon_n\right)}{\ln{2}}
\end{align}
\)
より、例えば初期の解の範囲(\(\varepsilon_0\))を1, 解を\(10^{-14}\)まで収束させようとすれば、大体\(n=47\)となり、47回、2分法を繰り返す必要があります。
逆に言えば、どんな状況でも47回繰り返せば相対誤差を14桁収束させることが出来るということです。

収束の判定は
要求精度\(\epsilon_{\text{tol}}\)、機械精度\(\epsilon_{\text{mac}}\)、解の存在範囲\(\delta\)、新たな解の推定値\(x\)、
とすると、
\(
\delta \lt 4\epsilon_{\text{mac}}|x|+\epsilon_{\text{tol}}
\)

もしくは
\(
f(x)=0
\)

で判定します。
この意味は、解の存在範囲を\(\pm 2\epsilon_{\text{mac}}\)の範囲まで特定するか、解の存在範囲が要求精度\(\epsilon_{\text{tol}}\)まで達するかのどちらかが満たされるかで判定します。


Fortranプログラムはこちら。

二分法でやると以下の通り20回必要になります。

$ gfortran bisection.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      5.2500000000000000
      2.8750000000000000      5.2500000000000000
      2.8750000000000000      4.0625000000000000
      2.8750000000000000      3.4687500000000000
      2.8750000000000000      3.1718750000000000
      2.8750000000000000      3.0234375000000000
      2.9492187500000000      3.0234375000000000
      2.9863281250000000      3.0234375000000000
      2.9863281250000000      3.0048828125000000
      2.9956054687500000      3.0048828125000000
      2.9956054687500000      3.0002441406250000
      2.9979248046875000      3.0002441406250000
      2.9990844726562500      3.0002441406250000
      2.9996643066406250      3.0002441406250000
      2.9999542236328125      3.0002441406250000
      2.9999542236328125      3.0000991821289062
      2.9999542236328125      3.0000267028808594
      2.9999904632568359      3.0000267028808594
      2.9999904632568359      3.0000085830688477
      2.9999995231628418      3.0000085830688477
      2.9999995231628418      3.0000040531158447
      2.9999995231628418      3.0000017881393433
      2.9999995231628418      3.0000006556510925
   3.0000000894069672    
$

6桁収束するために24回も評価回数が必要です。
計算時間を気にしないのならば良い方法ですが、やはり遅いです。

挟み撃ち法(False position法)



要点

・解が確実に見つかります。
・勢いよくゼロ点を横切っていない場合、収束が遅くなります。

計算方法

二分法との違いはゼロ点の推測に用いる情報量です。
二分法では位置\([a,b]\)における関数の符号だけを見てゼロ点を推測しますが、
挟み撃ち法では位置\([a,b]\)における関数のも見てゼロ点を推測しています。
この挟み撃ち法の推測は、点\((a,f(a)),(b,f(b))\)を結ぶ直線が通るゼロ点して値を推測するのです。

一つ、挟み撃ち法の説明を見ていた時に、解の囲い込み方法が2通りあるようです。
オリジナルのものは片方の点を絶対に変えない方法のようですが、二分法と同じように決めている物もあります。

ここでは、二分法と同じように決めていく後者の方法を載せておきます。
二分法の新しい点を決める所だけ変えればプログラムは完了です。

要求精度を6桁にすると、実行結果は

$ gfortran false.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      4.8581311942172727
      2.6310989812608572      4.8581311942172727
      2.6310989812608572      3.1997221817009311
      2.9967206978247356      3.1997221817009311
      2.9967206978247356      3.0000354381074144
      2.9999999998952847      3.0000354381074144
   3.0000000000000000    
$

となり、初めの1回は無視すると6回で収束に至っていることが分かります。早いですね。

Brent法



要点

・解が確実に見つかります。
・基本的には逆二次補間(inverse quadratic interpolation)で根を探索します。
・逆二次補間が原理的に出来ない場合(二つ以上の関数値が同じ場合など)では挟み撃ち法を行います。
・”探索が失敗”した時は二分法に切り替えて根を探します。
・1変数の場合、実用的な方法として推奨されています。

計算方法

逆二次補間は3点を結ぶ逆二次関数(上下に凸ではなく、左右に凸をもつ関数)としてラグランジュ補間をして、その根をゼロ点の位置として用いる方法です。
しかし、逆二次補間を行う際に必要な3点の内、2点以上が同じ関数値を持つと解が発散する、という問題があるため、その時は挟み撃ち法に切り替えて計算を行います。

二分法は、逆二次補間や挟み撃ち法を用いた為に解の収束が遅くなる、と判定されたときに使われます。

実際に動かしてみると逆二次補間よりも挟み撃ち法が使われる頻度が多いようです。

詳細は(Van Wijngaarden-Dekker-Brent法, William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery著, 丹慶勝市ら訳(1993) 『Numerical Recipes in C [日本語版]C言語による数値計算のレシピ』 技術評論社 261pp.)に書いてありますので、こちらを参照してください。

プログラム自体はパブリックドメインのNetlibにあるもの(zeroin.f)を基本としましょう。ただし、Fortran77で掛かれているので並列計算を考えた時に厄介です。文番号を消したものを置いておきます。

要求精度を6桁にすると、実行結果は

$ gfortran main.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      4.8581311942172727
      2.6310989812608572      4.8581311942172727
      2.6310989812608572      3.5172292241144740
      2.6310989812608572      3.0069570368988154
      2.9997489497326733      3.0069570368988154
      2.9997489497326733      3.0000000032534229
      2.9999995032534215      3.0000000032534229
   3.0000000032534229    
$

となり、初めの1回は無視すると7回で収束に至っていることが分かります。挟み撃ち法とだいたい同じです。

Netlibがpublicdomainであることの記述は以下の本にも見ることが出来ます。
https://books.google.de/books?id=2GPNBQAAQBAJ&pg=PA109&lpg=PA109&dq=netlib+public+domain&source=bl&ots=d6Uw8D5JMv&sig=Dw3wXTAQThFKLVrL7pm_qnCaq9s&hl=de&sa=X&ved=0ahUKEwjO05Xxt_XKAhXCwxQKHUBlC5sQ6AEIXjAI#v=onepage&q=netlib%20public%20domain&f=false

ちなみにですが、wikipediaにあるBrent法の疑似コードを実装したのですが、うまく働きませんでした。
Brent’s method -wikipedia en
ブレント法 -wikipedia ja

Newton法



要点

・解近傍のただ1点から収束させるので必要な情報が少なく済みます。
・解近傍のただ1点上における関数の微分を用いて収束させ、一回の計算で一致桁数が倍になります。
・本当に解の近くから始めないと、失敗します。
・多変数への拡張が容易ですが、関数の微分を数値微分で代用すると効率が落ちます。

計算方法

ニュートン法について詳しくは
ニュートン法(1、2次元、多次元)
に書いたので、こちらを参照してください。


プログラムは以下のもので、

実際に実行すると

$ gfortran newton.f90
$ ./a.out
      4.0000000000000000
      2.4339000410871483
      3.0980975267667068
      2.9994762813770324
      3.0000000000829825
   3.0000000000000000
$

という結果を得ます。

番外編:Steffensen法(Steffensen’s method)


Steffensen法は自己無撞着方程式を解くために使われる方法で、Newton法に似ています。
特に良く使われるのは、方程式がたまたま
\(
x=g(x)
\)

と書けている特定の方程式を満たす\(x\)を探す場合に使われます。

ですが、
\(
f(x)=x-g(x)
\)

を考えれば
\(
f(x)=0
\)

を考える問題と同じ、ということになります。
ここでは、方程式
\(
x=x^2/2
\)

を満たす\(x\)を探す問題を考えます。

この問題は
\(
x-x^2/2=0
\)

と同じです。

計算手法は
htt://park.itc.u-tokyo.ac.jp/kato-yusuke-lab/nagai/note_070420_self.pdf
でまとめられている通り、

\(
\begin{align}
a&=g(x_n), \\
b&=g(a),\\
x_{n+1}&=x_n-\frac{(a-x_n)^2}{b-2a+x_n}
\end{align}
\)
の漸化式で求められます。

プログラムではこちら。

複数のゼロ点を見つけたい場合



複数のゼロ点を見つけたい場合、一つの方法として、

解のある区間の特定→収束

を繰り返すことで得られる、と考えられるでしょう。

具体的には計算区間を荒く区切り、符号が変わる範囲を特定、その符号が変わる範囲に対してのみゼロ点を探すのです。

上の図の縦線は等間隔に荒く探した時の、関数が正になる時に青、負の時に赤色に表示しています。
収束した結果が緑の点です。

fortranプログラムはこんな感じになります。

ニュートン法(1、2次元、多次元)

ニュートン法に関するお話です。
”ニュートン法”と呼ばれる方法は
2種類(1.初期値近傍の極小値を求めるニュートン法、2. ゼロ点を求めるためのニュートン法)
があるようですが、ここでは2. 関数のゼロ点を求める方のニュートン法についてのお話です。

  1. まとめ
  2. 1次元のニュートン法
    1. 初期値の推定
    2. 導関数の近似
  3. 1次元ニュートン法のプログラム
  4. 2次元ニュートン法
  5. 2次元ニュートン法のプログラム
  6. 多次元ニュートン法
  7. 高次のニュートン法?
  8. 補)ニュートン法が使われる問題
  9. 参考文献

まとめ

1次元の場合

\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

ここで、\(x_0\)は解の近傍であること。

2次元の場合

\(
\begin{align}
\left(
\begin{array}{c}
x_{n+1} \\
y_{n+1}
\end{array}
\right)
=
\left(
\begin{array}{c}
x_n \\
y_n
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_n,y_n) \\
g(x_n,y_n)
\end{array}
\right)
\end{align}
\)

\(
\begin{align}
&\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1} \\
&~~~~=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

ここで、
\(
\begin{align}
f_x(x_n,y_n) &= \left.\frac{\partial f(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& f_y(x_n,y_n) = \left.\frac{\partial f(x,y)}{\partial y}\right|_{x=x_n,~y=y_n} \\
g_x(x_n,y_n) &= \left.\frac{\partial g(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& g_y(x_n,y_n) = \left.\frac{\partial g(x,y)}{\partial y}\right|_{x=x_n,~y=y_n}
\end{align}
\)

多次元の場合

\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

※\((~~)_n\)の添え字nはn回目の繰り返しにおけるベクトル\(\mathbf{x}\)を表しています。


スポンサーリンク


ここから本文

1次元のニュートン法


方程式
\(
f(x)=0
\)

を満たす\(x=a\)を求めることを考えます。

まず、この問題をグラフ
\(
y=f(x)
\)

のゼロ点(\(y=0\)となるような\(x=a\))を探す問題だ、という問いに置き換えます。

これを、解近傍で関数をテーラー展開で表し、その展開式のゼロ点を探すことで解を近似することで方程式を解きます。

関数\(f(x)\)が与えられたとき、任意の点\(x’\)周りにおけるテーラー展開は、
\(
f(x) = f(x’)+f'(x’)(x-x’) + O(\Delta^2),~~\Delta=|x-x’|
\)

と書くことが出来ます。もしも\(x\)と\(x’\)が近接していれば、\(\delta=|x-x’| \ll 1\)となるので、\(\Delta^2\)の項は無視できるほど小さくなると期待します。

さて、この式を導く際に用いた条件は、\(x\)と\(x’\)が近くにある、という条件のみです。なので、求めたい解\(a\)と\(a\)の近くの適当な点\(x_0\)を用いれば、
\(
\begin{align}
f(a) &= f(x_0)+f'(x_0)(a-x_0) +O(\Delta^2)~~~…(1)\\
f(x_0) &= f(a)+f'(a)(x_0-a) +O(\Delta^2) ~~~…(2)
\end{align}
\)

の2通りの式を考えることが出来ます。(2)に含まれる\(f'(a)\)は本当の\(a\)が分からないと評価するのが難しいので、(1)を選びます。

(1)の左辺は条件よりゼロです。残りの右辺の\(a\)について解けば、

\(
\begin{align}
0 = f(x_0)+f'(x_0)(a-x_0)+O(\Delta^2) \\
\to a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\end{align}
\)


ここで、\(f(x_0)=f'(x_0)=O(\Delta^0)\)を想定しているため、\(f'(x_0)\)の割り算を行っても誤差のオーダーは変化せず、\(\Delta\)に対して2乗のままとなります。

\(O(\Delta^2)\)を無視すると、得られる解は近似値となります。この近似値を\(x_1\)と置くと
\(\displaystyle
x_1 = x_0-\frac{f(x_0)}{f'(x_0)}
\)

得られる解\(x_1\)は\(x_0\)よりも解\(a\)に近いため、\(O(\Delta^2)\)はさらに小さくなります。なので、もう一度計算します。この時の近似値を\(x_2\)と置くと
\(\displaystyle
x_2 = x_1-\frac{f(x_1)}{f'(x_1)}
\)

ここで得られた\(x_2\)は更に\(a\)に近いため、・・・と繰り返していきます。

すなわち、解\(a\)に近い初期値\(x_0\)を与えた時、漸化式
\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

を繰り返すことで解\(a\)に近づいていく、ということです。
この式から分かる通り、解近傍で導関数がゼロになる場合はゼロ割発生するかもしれないので危険です。

もしもこのステップが無限回繰り返されれば、無限回目のステップの近似値\(x_\infty\)は
\(
x_\infty=a
\)

に収束します。

以上のように導関数を用いて方程式の解を求める方法は、ニュートン法(またはニュートン=ラフソン法)と呼ばれます。
この方法は、解に近いに所の値とその導関数が分かりさえすればいいので、この方法は1次元、多次元でも使用する事ができ、また複素関数でさえ使用することが出来ます。
さらに、誤差が\(O(\Delta^2)\)で減衰していくため、一度繰り返すごとに厳密解に一致する桁数が2倍になっていきます。

しかし問題があり、
1. 初期値を解の近傍に取らないと失敗する
2. 導関数が分からないと使えない
という大きな問題を抱えています。

また、解近傍でテーラー展開可能でなければならず、解近傍で導関数がゼロになる点があってはいけません。

  • 初期値の推定


    ニュートン法は解の近傍で、テーラー展開可能な領域から始めなければ用いることは出来ません。
    これは工夫をするしかありません。

    解近傍の初期値を推定する方法として、以下の2つの方法が良く取られます。
    1. 解の範囲を限定し、その近傍からニュートン法を行う(1次元のみ)。
    2. 適当な項を落とせば方程式が解けてしまう場合、その解を出発点とし、残りの項を摂動的に加えながら逐次解く。

    1. 解の範囲を限定する

      これは1次元の場合に使える方法です。
      単純な考え方で、例えば広い区間\([x_a,x_b]\)にゼロ点があることだけが分かっているとします。
      この区間を\(N+1\)分割し、
      \(
      x^{(0)}=x_a,x^{(1)},\cdots,x^{(k)},x^{(k+1)},\cdots, x^{(N-1)},x^{(N)}=x_b
      \)

      とします。もしも\(x^{(k)}\)と\(x^{(k+1)}\)の符号が違う場合、この区間でゼロを横切っているということになりますので、この点からニュートン法を始めれば良いのです。

      問題は、解の範囲を狭めるために関数を\(N+1\)回評価しなければならない、という点で計算コストがかかってしまう、という点です。

    2. 既知の解から摂動的に逐次解く

      これは若干特殊な方法で、あてずっぽうな上の考えとは違います。
      本来解きたい式を
      \(
      f(x)=f_0(x)+g(x)
      \)

      と分けます。ここで、
      \(
      f_0(x)=0
      \)

      の解は簡単に解けて、\(x=a_0\)が関数\(f_0(x)\)のゼロ点であるとします。この\(a_0\)を初期値とし、\(\lambda\)を十分小さい値にとり
      \(
      f_0(x)+\lambda_1 g(x)=0
      \)

      を解きます。ここで、\(\lambda_m\)は\(0\lt\lambda_1\lt\lambda_2\cdots \lt\lambda_{M-1} \lt \lambda_M=1\)を満たすとします。
      この方程式の解\(a_1\)が得られたら、この解を初期値とし、
      \(
      f_0(x)+\lambda_2 g(x)=0
      \)

      を解きに行き、これを\(\lambda_M=1\)まで続けます。
      こうすることにより、最後の答え\(a_M\)がもともと解きたかった式の答えとなっているのです。

      関数の一回の評価に時間がかかってしまう場合、この解き方の方が早く解けます。
      しかし、関数\(g(x)\)によってゼロ点の位置が増えたり、消えたりしてしまう時にこの方法は使うことが出来ず、連続的に解が動いてくれないと使うことが出来ない点に注意してください。

  • 導関数の近似


    導関数を求める簡単な方法は、差分を使うことです。すなわち、導関数\(f'(x)\)を差分に置き換えます。
    \(\displaystyle
    f'(x)\approx \frac{f(x+h)-f(x)}{h} + O(h)
    \)

    すると、漸化式は
    \(\displaystyle
    x_{n+1}=x_{n}-h\frac{f(x_{n})}{f(x_n+h)-f(x_n)},~~n=0,1,2\cdots,
    \)

    と得られます。
    この方法は導関数を必要とせず、実装が簡単であるという意味で有用な式です。

    導関数を求める際の前進差分の刻み幅\(h\)は、計算機の扱える丸め誤差\(\epsilon\)(単精度ならば約\(\epsilon=10^{-8}\)、倍精度ならば約\(\epsilon=10^{-16}\))の\(1/2\)乗の2倍、つまり
    \(
    \displaystyle h\sim 2\sqrt{\epsilon}
    \)

    とすると良いです。

    この刻み幅の導出は折りたたんでおきますが、以下のように導出することが出来ます。

参考文献[1][3][5]

1次元ニュートン法のプログラム


以下のプログラムは、1次元のニュートン法のプログラムで、
方程式
\(
x^2-4 = 0
\)

を初期値\(x_0=3\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに急激に解に近づいていっていることが分かります。
(デフォルトではコメントアウトしてあります。)

$ gfortran main.f90
$ ./a.out
    1       2.1666666616021075
    2       2.0064102565311974
    3       2.0000102400381690
    4       2.0000000000262466
    5       2.0000000000000000
   2.0000000000000000       ---result
$

上には出力していませんが、初期値の評価もしています。
なので、ニュートン法は6回呼び出されたことになります。
ニュートン法1回当たり関数は2回呼び出されますので、12回関数が呼び出されたことになります。

例えば比較対象として、二分法を考えますと、1回目に関数は3回,そのほかで1回ずつ呼び出されます。
二分法の解は1ステップ当たり区間の半分になりますので、16桁一致するまでには50回ほど繰り返す必要があります。そのため、52回位関数が呼び出されます。

ニュートン法の方が圧倒的に早いことが分かります。

複素関数の1次元ニュートン法のプログラム
$ gfortran main.f90
$ ./a.out
                    1 ( -1.2933333304320405     , 0.72000000272341036     )
                    2 (-0.78207788228079722     , 0.60930726266833468     )
                    3 (-0.43844290526861285     , 0.73503785880633421     )
                    4 (-0.50851135934742864     , 0.89043164725052448     )
                    5 (-0.50009896196621662     , 0.86666395460253720     )
                    6 (-0.49999991066297655     , 0.86602581130981315     )
                    7 (-0.49999999999985084     , 0.86602540378453374     )
                    8 (-0.50000000000000000     , 0.86602540378443860     )
 (-0.50000000000000000     , 0.86602540378443860     )  ----result
$

ソースコードはこちら。

2次元のニュートン法


続いて方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y)&=0,\\
g(x,y)&=0
\end{aligned}
\right.
\end{eqnarray}
\)

を満たす\(x=a,y=b\)を見つけることを考えます。

1次元の時と同様に、任意の点\(x’,y’\)周りにおけるテーラー展開は、
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y) &= f(x’,y’)+f_x(x’,y’)\cdot(x-x’)+f_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)\\
g(x,y) &= g(x’,y’)+g_y(x’,y’)\cdot(x-x’)+g_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

ここで、
\(\displaystyle \left. f_x(x’,y’)=\frac{\partial f(x,y)}{\partial x}\right|_{x=x’,y=y’}\)
\(\displaystyle \left. f_y(x’,y’)=\frac{\partial f(x,y)}{\partial y}\right|_{x=x’,y=y’}\)
\(\displaystyle \Delta x=|x-x’|=O(\varepsilon^1),~~\Delta y=|y-y’|=O(\varepsilon^1)\)
と書きました。

1次元の時と同様に初期値\((x_0,y_0)\)が解\((a,b)\)に近ければ、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
0 &= f(x_0,y_0)+f_x(x_0,y_0)\cdot(a-x_0)+f_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)\\
0 &= g(x_0,y_0)+g_x(x_0,y_0)\cdot(a-x_0)+g_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

を同時に満たす\((a,b)\)を探せばよい、ということになります。

行列形式で書けば
\(
\begin{align}
\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\left(
\begin{array}{c}
a-x_0 \\
b-y_0
\end{array}
\right)
=
-\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

と書けます。\(a,b\)について解けば
\(
\begin{align}
\left(
\begin{array}{c}
a \\
b
\end{array}
\right)
=
\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

となります。丁度一次元のニュートン法と似た形になりました。

右辺の
\(
\begin{align}
\mathbf{J}_0=\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\end{align}
\)

は\(x=x_0,y=y_0\)のヤコビアンです。また、
\(
\mathbf{a}=\left(
\begin{array}{c}
a \\
b
\end{array}
\right),~~
\begin{align}
\mathbf{x}_0=\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right),~~
\mathbf{f}_0=\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)
\end{align}
\)

と表記すれば、2次元のニュートン法は
\(
\mathbf{a}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0+O(\varepsilon^2)
\)

と書くことが出来ます。
近似として、\(O(\varepsilon^2)\)の項を無視すれば
\(
\mathbf{x_1}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0
\)

と書けます。繰り返せば、漸化式
\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

を解いていけば良い、ということになります。

ニュートン法を実行するにあたり、ヤコビアンの逆行列を解くことが一番の問題となります。

良く知られているように、2次元の逆行列は、計算コストもさほどかからずに解けてしまい、
\(
\begin{align}
\mathbf{J}_n^{-1}=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

として求める事が出来ます。

参考文献[2],[6]

スポンサーリンク

2次元ニュートン法のプログラム


以下のプログラムは、2次元のニュートン法のプログラムで、
方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f_1(x,y)&=x^2+y^2-1 = 0 \\
f_2(x,y)&=y-x^3 = 0
\end{aligned}
\right.
\end{eqnarray}
\)

を初期値\(x_0=2,~y_0=1\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに解に近づいていっていることが分かります。

$ gfortran main.f90
$ ./a.out
    1       1.3571428532359113       0.2857142813732352
    2       0.9844126813492445       0.4401111876910193
    3       0.8485699953676349       0.5590406123395349
    4       0.8265084709901402       0.5633730814110012
    5       0.8260315915076198       0.5636240770681390
    6       0.8260313576542442       0.5636241621612319
    7       0.8260313576541870       0.5636241621612585
  0.82603135765418700       0.56362416216125855       ---result
$

この問題の答えは\(f_1=0,~f_2=0\)の陰関数の交点です。
図示すれば、下図の紫色の経路を辿って解に収束していくことが分かります。

ひとつ、注目しておかなければならない点として、解に近いからと言って、その解に収束するとは限らないということに注意してください。
図中の水色の線は初期値こそ右上の解に近いですが、実際に計算してみるともう一つの解に収束していることが分かります。

適当な初期値を用意し、ニュートン法が収束する点を調べるという問題はニュートン写像と呼ばれる分野です。
力学系とかそういう言葉が出てきます。
私も過去に複素関数でやったことがありますので、リンクだけ載せておきます。
ニュートン写像に現れる綺麗な画像 -シキノート

2次元の複素関数のニュートン法

置いておきます。

多次元のニュートン法

さて、この形まで持ってくると多次元への拡張が簡単にできます。
多次元のニュートン法は
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

の形がそのまま使えます。
ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

です。

問題となるのは(ヤコビアンの逆行列)×(f)の計算です。
この計算を数値計算的に何とかするべく改良が重ねられてきたのが準ニュートン法(本稿では書きません)と呼ばれるアルゴリズムです。

逆行列を数値的に計算することは宜しくないですので、
\(
\begin{align}
\mathbf{r}&=\mathbf{J}_n^{-1}\mathbf{f}_n \\
\to \mathbf{J}_n \mathbf{r} = \mathbf{f}_n
\end{align}
\)

という方程式を解く問題に帰着させます。これが解ければ、ニュートン法の次のステップを
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{r}
\)

として求める事が出来るのです。

高次のニュートン法?


1次元のニュートン法は
\(\displaystyle
a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\)

の形で求めていきますが、二階微分が分かっていればもっと収束が早くなりそうです。
実際、
\(\displaystyle
a = x_0-\frac{f^{\prime}(x_0)}{f^{\prime\prime}(x_0)}\pm\sqrt{{f’}^2(x_0)-2f(x_0)f^{\prime\prime}(x_0)}+O(\Delta^3)
\)

と解けます。しかし、2次関数ですので
解が見付からないかもしれない
という可能性があります。1次のニュートン法であれば、導関数がゼロでない限り、必ず\(y=0\)と交わる点が見付かります。しかし、2次で解に近づけると交点が見つからないことがあるかもしれません。

なので、1次のニュートン法を繰り返して使う方が(2次と比べたら)安定な解法なのです。

参考ページ

[1]1 Newton 法:壱変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node6.html
[2]2 Newton 法:弐変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node7.html
[3]Newton法
http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-suchi/newton130805.pdf
[4]アルゴリズムによる誤差
http://www.aoni.waseda.jp/ykagawa/chap2html/node9.html
[5]Newton 法による方程式の近似解法
http://www.math.u-ryukyu.ac.jp/~suga/C/2004/7/node9.html
[6]多次元のニュートン・ラフソン法
http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%C2%BF%BC%A1%B8%B5%A4%CE%A5%CB%A5%E5%A1%BC%A5%C8%A5%F3%A1%A6%A5%E9%A5%D5%A5%BD%A5%F3%CB%A1
[7]山本 有作, 行列計算における高速アルゴリズム
http://www.cms-initiative.jp/ja/events/0620yamamoto.pdf
[8]中田 和秀, 大規模線形方程式を解くためのクリロフ部分空間法の前処理
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1288-6.pdf

微分方程式を解く際の計算量の工夫

運動方程式を解いて、ある特定の時刻の値が欲しいとします。
特定の時刻がいつなのか、あらかじめ分かっている場合、与えられた初期条件から特定の時刻まで計算し、データを貯めておいて後で使えば良いです。

特定の時刻が分からない時に計算量の問題が生じます。
例えば、時刻\(t=0\)の初期条件が与えられた時、時刻\(t_1=100,t_2=101, t_3=102\)の時刻が知りたいとします。ここで、時刻\(t_2\)の値は実は\(t_1\)の関数の値を評価してから初めて決まるとします。

この場合の計算方法としては、

  1. 毎回、時刻\(t_0\)から初めて\(t_1,t_2,t_3\)まで計算
  2. あらかじめ微分方程式を一度細かく計算しておき、欲しい時刻\(t=t’\)を補間で求める
  3. あらかじめ微分方程式を一度細かく計算しておき、欲しい時刻\(t=t’\)に近い点から計算

だと思います。最も簡単に実装できるのは1の方法で、欲しい時刻が数点であれば計算時間は問題にならないほど少なくなるでしょう。2の方法は精度を気にする場合には使用できません。ここでは3番目の方法を考えます。

計算方法を図示したものが以下のものです。

初期条件からいちいち出発して\(t=t’\)に到達するのではなく、あらかじめ計算していた参照点から出発して計算する、ということです。
ここで問題なのが、もし横軸が時刻だった場合、時間を逆向きに遡って良いか?という問題です。

これは物理の問題で言えば、解こうとしている問題は時間可逆か?ということで、言い換えれば系のエネルギーは保存しているか?ということです。もしも抵抗力など、エネルギーが散逸する場合は、いくら参照点が近いからといって逆向きに進んではいけないということです。

ここでは、微分方程式を解く際に戻ることはせず、必ず求めたい時刻よりも前の参照点からスタートする、というプログラムを作ります。

以下のように実装できます。

シンプルQUADPACK

QUADPACKのコード(http://www.netlib.org/quadpack/)は優秀ですが、2000行近くあって色々いじる際に面倒です(最速・高精度の数値積分)。

コードを覗くと、求積法の点数だけが違うプログラムがあるので長いことに気が付きます。
なのでコード量を減らしましょう。

使用する積分方法を20-41点ガウス求積法(key=4)に限定しました。key=4を選んだ理由は、最速・高精度の数値積分の結果から推測しました。

key=2の低次の方法(10-21点)であると滑らかな関数でも評価回数が多くなる反面、不連続点があっても少なく済みます。
key=6の高次の方法(30-61点)であると滑らかな関数で評価回数が少なくなる反面、不連続点があると多くなります。

QUADPACKにはkey=1,2,3,4,5,6のどれかを選べますが、基本的に滑らかな関数の積分を行うことが多いので若干高次のkey=4を選びました。

350行程度にまで減らし、ルーチンを一つ無くました。変更箇所が怖い方は、オリジナルのものを使用してください。

コードの変更をした結果、計算結果が違くなるなど問題が起こる場合、原因は私にあります。
著作者とは関係ありませんので、その点はよろしくお願いします。

——————————————

実軸上、複素関数でkey=4の場合

——————————————

key=4, 複素平面上の積分

——————————————
Key=1の場合
(7-15点ガウス求積法)
・不連続性や、端点特異性がある時に有利です。

exp(-ikx)/x の無限区間に渡る積分を数値的に計算したい。

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)


1/xのフーリエ変換である
\( \displaystyle
\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

数値的に計算します。

この積分値は一意には決まりません。
複素関数論によると、この積分値は\(x=0\)周りに存在する特異点を上周りに迂回して積分するか、下周りに迂回して積分するか、で答えが異なります。

特異点を上周りに回る場合、答えは
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
-2\pi i~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
であり、特異点を下周りに回る場合、
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x-i\varepsilon}dx =
\left\{
\begin{aligned}
0~~~(k\gt 0)\\
2\pi i~~~(k \lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
となります。
私が疑問に思うのは、積分すると本当にステップ状の関数が現れるのか?ということです。にわかには信じられません。

ここでは、上を回った時の計算にのみ注目します。
答えに出てくる係数を減らしたいため、被積分関数
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx
\)

を考えます。

数値計算方法


さて、問題を整理しましょう。今取り組む問題は、

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

です。読み取れる情報は

・始点は、実軸上のマイナス無限大
・終点は、実軸上のプラス無限大
・\(x=0\)に特異点が存在する

という情報です。これらは理由なく勝手に変えてはいけません。
被積分関数が収束するならば、xがどんなに大きくても無限遠まで振動していきますが、フレネル積分のように積分結果は収束することを期待しましょう。

さて、\(x=0\)に特異点があるために、実軸上のマイナス無限からプラス無限まで積分していこうとすると\(x=0\)を通過しようとする時に無限大になるため、オーバーフローします。
ということで、積分経路を特異点周りで迂回します。

始点と終点は決められていますが、積分経路の指定はありませんので何も問題ありません。
積分経路は以下の図のように取ります。

このようにとると、積分経路上に特異点はないですので計算することが出来ます。
始点と終点は十分大きい値、ここでは(±2000,0)という値にします。

横軸をkとして、積分を計算すると以下のようになります。

複素関数論から、特異点を複素平面上方を迂回した場合に得られる解は
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} -\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
1~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
です。数値計算によって得られた解は、\(k\lt 0\)で値0, \(k\gt 0\)で値1であるので、数値的にも正しそうです。
しかし、厳密にそれぞれの値になっているわけではないことに気が付きます。なぜならば、\(k=0\)周りでは値は収束していないように見えます。
この原因は始点と終点としてとった、大きい値(±2000,0)が十分ではない事を意味しています。kが小さいと、被積分関数の\(\exp(-ikx)\)による振動はゆっくりになるため、収束が遅くなっている、と考えられます。

ジョルダンの補題の利用


この振動を無くすためには被積分関数が漸近領域で収束するように積分経路を変えればいいのです。

始点と終点は変わりません。図は\(k\)が負の時の経路です。

経路の取り方は無数にあり、上の経路は一つ前の計算で選んだ経路を出来る限り変更しないように経路をとっているだけにすぎません。

ジョルダンの補題を利用すれば、被積分関数が収束する領域であれば積分の寄与はゼロなので、数値的に積分する場合の視点、終点を複素平面上へ移すことが出来ます。
この考えの下実際に計算すると積分結果はこうなります。

綺麗なエッジが作られました。これならば解析解と数値解が一致した、と言えるでしょう。複素平面上の数値計算は一筋縄ではいかないことが実感できました。

数値計算プログラム


プログラムは以下のものを用いています。

数値積分はquadpackのガウス=ルジャンドル求積法の自動積分コードを複素関数用に変更して用いています。
実際に使用しているquadpackコードは以下のもので、
http://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag.f90
これと、下のはmainプログラムです。

program main
  implicit none
  integer::ier,Nt,i,j,sn
  double precision::eps
  complex(kind(0d0))::a,b,s,s0
  complex(kind(0d0))::xa(1:5),xb(1:5)
  complex(kind(0d0)),external::g
     
  eps=1d-10
  ier=1

  Nt=1000
  do i=1,Nt
     t=i*(1d0-(-1d0))/dble(Nt)-1d0
     s=dcmplx(0d0,0d0)
     
     sn=1
     if(t.gt.0d0)sn=-1
     if(t.eq.0)sn=0
     
     xa(1)=dcmplx(-2000d0,sn*2000d0)
     xb(1)=dcmplx(-1d0,0d0)

     xa(2)=dcmplx(-1d0,0d0)
     xb(2)=dcmplx(-1d0,1d0)

     xa(3)=dcmplx(-1d0,1d0)
     xb(3)=dcmplx( 1d0,1d0)

     xa(4)=dcmplx( 1d0,1d0)
     xb(4)=dcmplx( 1d0,0d0)

     xa(5)=dcmplx( 1d0,0d0)
     xb(5)=dcmplx( 2000d0,sn*2000d0)  

     do j=1,5
        s0=dcmplx(0d0,0d0)
        call cqag_sk(g,xa(j),xb(j),eps,s0,3,ier)
        s=s+s0
     enddo
     write(10,*)t,dble(s),dimag(s)
  enddo  
     
  stop
end program main

function g(z)
  implicit none
  complex(kind(0d0))::g
  complex(kind(0d0)),intent(in)::z

  double precision::pi=dacos(-1d0)
 
  g=exp(-dcmplx(0d0,1d0)*z*t)/(z)
  g=g*(-1d0/(2d0*pi*dcmplx(0d0,1d0)))

  return
end function g

オリジナルのquadpackコードは
http://www.netlib.org/quadpack/(パブリックドメイン)です。

時間依存シュレーディンガー方程式の数値解法 クランク=ニコルソン法

解きたい問題は以下の形の時間依存シュレーディンガー方程式です。

これをクランクニコルソン法という、微分を差分で近似する数値的な時間発展方法で求めていきます。

解法


※途中、式番号が抜けたり(例えば式15)、前後していますが仕様です。

時間を初期時刻\(t=0\)から\(t=\Delta t\)まで時間発展したいと考えます。
微分を差分で離散化しますが、時刻\(t=\frac{\Delta t}{2}\)だけ進んだ時刻からの前進差分をまず考えます。
すると、

と導けます。時間,空間に関して\(O(\Delta t),O(\Delta x^2)\)です。

続いて、時刻\(t=\Delta t\)から\(\frac{\Delta t}{2}\)との後進差分をまず考えますと

です。時間,空間に関して\(O(\Delta t),O(\Delta x^2)\)です。

式(1)と式(2)を足して2で割りますと、時間に関して中央差分を取っていることになりますので、次数が上がります。代入しますと、

となります。時間,空間に関して\(O(\Delta t^2),O(\Delta x^2)\)です。

左辺に時刻\(t+\Delta t\)の波動関数をまとめますと、

となります。これで離散化が出来ました。これがクランクニコルソン法です。

式が複雑なので、少し簡単にまとめます。
位置\(x\)を等間隔に離散化して、以下の通り離散化と記述を変えます。

この表記に倣いますと、離散化した時間依存シュレーディンガー方程式は簡潔に書けて

と書けます。
ここで、係数\(a,b,c,r\)は式(5)の各係数をまとめたものですが、境界条件によって多少異なります。
ここでは、固定端条件と周期境界条件の場合を特に考えてみましょう。

固定端条件の場合


固定端条件は両端での関数の値が時間によって変化しないとする条件です。数式で書けば

という条件が波動関数に課されます。この時、行列表式で式(8)を記述すれば、

と書き表すことが出来ます。ここで、\(\psi_0^{(+)},\psi_N^{(+)}\)に掛かる係数は境界条件を満たすように決定しています。行列の演算で知られているように、第1行目を\(a_1\)倍して第2行目を引き、第N行目を\(c_{N-1}\)倍して第N-1行目を引けばオレンジで四角く囲った枠内のみを計算すればよいことが分かります。
各々の係数の明らかな形は以下の通り得られます。

周期境界条件の場合


周期境界条件は両端での関数の値が等しく、その微分も等しいという条件です。数式で書けば

という条件が波動関数に課されます。この時、行列表式で式(8)を記述すれば、

と書き表すことが出来ます。ここで、\(\psi_N^{(+)}\)に掛かる係数は境界条件を満たすように決定しています。行列の演算で知られているように、第N行目は既に解けているのでこの要素を考える必要はありません。なので、オレンジで四角く囲った枠内のみを計算すればよいことが分かります。
各々の係数の明らかな形は以下の通り得られます。


スポンサーリンク


三重対角行列の方程式の解法


さて、あとは式(9),(13)を数値的に解けば良いのです。
固有値問題ではありませんので、連立方程式を解けば良いのです。

固定端条件の式(9)を解くためのアルゴリズムはThomasのアルゴリズムと呼ばれる方法が一般的です。計算量のオーダーは\(O(N)\)です。

周期境界条件を含んだ三重対角行列の解法はLU分解を利用した方法[2]や、Sherman-morrison formulaという公式を利用した方法[3]があります。ここでは、擬コードが載せられている[2]を元に作成していきます。これも計算量のオーダーは\(O(N)\)ですが、固定端条件の場合よりも約二倍ほど計算量が多くなります。

メモとして書いておきますが、周期境界条件を含む5重対角行列の解法もありまして、[1]に見ることが出来ます。

フォン・ノイマンの安定性解析


これは、クランク=ニコルソン法のような線形変微分方程式を有限差分法で解く際に、数値的に安定なのか発散するのかを調べる手法です。
注意しなければならないのは、数値的に安定なので答えも正しい、という事ではない事です。
この解析により、有限差分を行う際に最低限満たさなければいけない空間の刻み幅と時間の刻み幅の条件が決定されます。

\(V(x)=0\)を考えます。
ざっくりと説明しますと、偏微分方程式における独立な解、すなわち固有の波数\(k\)で特徴づけられる解が発散しない条件から導かれます。
ある時刻の解が1ステップの時間発展をした時、複素数の増幅因子(amplification factor) \(\xi\)を伴って

\(
\begin{align}
\psi_j &= \xi^0 e^{ikj\Delta x}\\
&\downarrow \\
\psi^{(+)}_j &= \xi^1 e^{ikj\Delta x} \\
\end{align}
\)

で時間発展する、とします。時間発展によって、複素数の整数乗が掛かるのです。
もしも増幅因子の絶対値、すなわち\(|\xi|\)が1より小さければ、時間発展しても独立な解は収束するため、数値的に安定だと言えます。
しかし、\(|\xi|\)が1より大きければ、独立な解の時間発展が指数関数で発散していく事を意味するため、数値的に不安定になります。

この\(\xi\)を求めるためには、波数\(k\)をもつ固有モードの式を有限差分で書き表した式(ここでは式(5))に代入すれば求められます。kはいかなる値を持ち得ると考えます。

計算の結果、\(\xi\)は以下の通り求められます[3]。
\(
\begin{align}
\xi=\frac{1-2\beta \sin^2\left(\frac{k\Delta x}{2}\right)}{1+2\beta \sin^2\left(\frac{k\Delta x}{2}\right)} \\
\beta=\frac{\alpha}{i}\frac{\Delta t}{(\Delta x)^2}
\end{align}
\)

※検算はしていません。[3]を参照しました。

この式から分かることは、位置の間隔\(\Delta x\)が与えられたとき、いかなる\(\Delta t\)であっても\(|\xi|\lt 1\)を取るので、数値的に安定である、ということが分かります。

これはクランク=ニコルソン法(陰解法)の特徴となっており、数値的には安定な計算方法です。

詳細は[3]や[4]をご覧ください。

数値計算コード


数値計算コードは以下のように書けます。
以下のプログラムは下の例2を計算するプログラムです。
実用上は、固定端、周期境界のどちらかだけを使えばいいので、適宜消してください。

スポンサーリンク

実例


自由粒子

初期状態:ガウシアン型, 運動量を持つ波束が時間とともにどう発展していくか。
\(
\begin{align}
i\frac{\partial }{\partial t}\psi(x,t)&=-\frac{1}{2}\frac{\partial^2 }{\partial x^2}\psi(x,t) \\
\psi(x,t=0)&=\left(\frac{1}{\pi a^2}\right)^{1/4}\exp\left[-\frac{x^2}{2a^2}+ik_0 x\right]
\end{align}
\)

\(
\begin{align}
\psi(x,t)&=\frac{1}{\pi^{1/4}}\left(\frac{1}{a+it/a}\right)^{1/2}
\exp\left[-\frac{1}{2}\frac{(x-k_0 t)^2}{a^2+it}+i k_0 x -i\frac{k_0^2}{2} t\right]\\
|\psi(x,t)|^2&=\frac{1}{\pi^{1/2}}\left(\frac{1}{a^2+(t/a)^2}\right)^{1/2} \exp\left[-\frac{(x-k_0 t)^2}{a^2+(t/a)^2}\right]
\end{align}
\)


上:存在確率密度\(|\psi(x,t)|^2\)。黒→解析解、青→固定端条件の数値解、赤→周期境界条件の数値解
下:解析解との相対誤差\(|(|\psi(x,t)|^2-|f(x,t)|^2)/|f(x,t)|^2|\)。

gnuplotコードはこちら

調和振動子の固有状態の重ね合わせ

初期状態:調和振動子の基底状態+第一励起状態
\(
\begin{align}
i\frac{\partial }{\partial t}\psi(x,t)&=\left[-\frac{1}{2}\frac{\partial^2 }{\partial x^2}+\frac{1}{2}x^2\right]\psi(x,t) \\
\psi(x,t=0)&=\frac{1}{\sqrt{2}}\left[\frac{2^{1/2}}{\pi^{1/4}}xe^{-x^2/2}+\frac{1}{\pi^{1/4}}e^{-x^2/2}\right]
\end{align}
\)

\(
\displaystyle \psi(x,t)=\frac{1}{\sqrt{2}}\left[\frac{2^{1/2}}{\pi^{1/4}}xe^{-x^2/2}e^{-i\frac{3}{2}t}+\frac{1}{\pi^{1/4}}e^{-x^2/2}e^{-it/2}\right]
\)

gnuplotコードはこちら


上:存在確率密度\(|\psi(x,t)|^2\)。黒→解析解、青→固定端条件の数値解、赤→周期境界条件の数値解
下:解析解との相対誤差\(|(|\psi(x,t)|^2-|f(x,t)|^2)/|f(x,t)|^2|\)。

参考文献


[1]Tomohiro Sogabe, New algorithms for solving periodic tridiagonal and periodic
pentadiagonal linear systems, Appl. Math. Comput. 202 (2008) 850-856

[2]A. A. Karawia, A computational algorithm for solving periodic pentadiagonal linear systems, Appl. Math. Comput. 174 (2006) 613–618

[3]Van Wijngaarden-Dekker-Brent法, William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery著, 丹慶勝市ら訳(1993) 『Numerical Recipes in C [日本語版]C言語による数値計算のレシピ』 技術評論社

[4]フォン・ノイマンの安定性解析 -wikipedia