二階線型常微分方程式の固有値問題の高精度数値解法と量子力学への応用

石川 英明


Return

1 序論

線型常微分方程式の固有値問題は境界値問題という数学と数理物理学の基礎の一つであり[1 - 4]、その重要な応用例は量子力学である[5 - 9]。量子力学の固有値問題で解が特殊或いは数学関数[1 - 3, 10 - 12]を使って解析的に得られるのは自由粒子、調和振動子、水素原子、等に限られている。摂動論、変分法、あるいはWentzel-Kramers-Brillouin (WKB) 法を用いて近似的に求まる場合があるが、一般に近似法の適用範囲は限られている。これらの制限を除くため、数値的に解くことが量子力学の初期から行われてきた。そうした先駆的な仕事は原子構造計算に実り多い知見をもたらしたが[13, 14]、固有値の計算誤差は1.0D-8 (= 1.0 × 10-8)より大きく、計算精度を更に改善する必要がある。このため任意のポテンシャルに対して現象を正しく理解するために必要な要求精度を満たす数値解を出すことが求められている。本報告では、最近発展した高精度かつ単純な数値計算手法とその応用を示す[15, 16]。
二点境界値問題の数値解法のひとつは、微分方程式の空間変数に関する微分を離散化して実空間での行列固有値問題に変換し、その固有値方程式を解く方法である[17 - 20]。この方法において最近急激な進歩があった[15]:一次元Schrodinger方程式を高精度の数値微分法に基づいて離散化し、行列固有値問題に変換して、行列固有値問題に対する標準的な手法を用いて解き、高精度の固有値と固有関数を得た。更に、中心差分積分法を用いてハミルトニアンの行列要素を計算し、対角要素が厳密な固有値と良く一致することを示した。この解法の利点は、高い量子数まで高精度であり、また固有値を見落とさず確実に見つけるので、計算の信頼性が高いことにある。一方、課題は計算時間が長いこと(特に行列の三重対角化に時間がかかる)である。
固有値問題のもう一つの解法は古典的なshooting法あるいはmatching法と呼ばれる方法である:二点境界値問題を、全区間の両端を出発点とする初期値問題に書き換え、固有値と固有関数の初期値を適切に与えて、両端から中央方向に向かって微分方程式を数値的に解き、それらを中間の点(接続点)で固有関数とその一階微分を滑らかに接続するように固有値を修正し、その修正値が許容誤差範囲内に収束するまで反復計算する。この場合の課題は、(i) 固有値と固有関数の初期値の求め方、(ii) 初期値問題の高精度かつ安定な数値解法、そして(iii)適切に選ばれた接続点における固有値の修正法にある。これらは以下のように解決された[16]。
第1の課題に対しては、離散化行列固有値法を使う。但し、初期値にはある程度の精度があれば良いので、2階微分を最も簡単な三点公式で離散化することで、計算時間がかかる三重対角化を省く。以上により信頼性の高い固有値と固有関数の初期値を得ることができる[15]。第2の課題に対しては、これまで多くの研究が、数値解析[21, 22]、数値計算ソフトウェア[23, 24]、微分方程式の数値解法[25 - 27]及びそこに引用されている文献、で報告されている。実際問題では、微分方程式を解くための物理量が離散点でのみ与えられている場合が多いので、離散点の中間での値を必要とする1段法(Runge-Kutta法が代表例)よりは、与えられた離散点での値のみを計算に用いる多段法の方が適している。多段法で精度と安定性に関するDahlquist [28, 29] とHenrici [25] の先駆的な研究を推し進めて、大きな段数の公式を系統的に導き、高精度を得てパフォーマンスを確かめた。この方法は低次の公式を使ったもの[30, 31]よりはるかに高精度である。第3の課題に対しては、よく知られた方法[32 - 35] で対処できる。数値微分と積分[15]及び初期値問題の解の精度向上と、接続点を適切に選ぶことで高精度の固有値と固有関数を得ることができた。
本報告の構成は以下の通りである。第2節では量子力学と数値計算の関係、第3節では高精度数値計算法として、補間、微分、積分、それらに基づいた離散化行列固有値法、shooting法では、固有値と固有関数の初期値、微分方程式の初期値問題、固有値の修正法を述べ、更に既存の方法との比較を述べる。第4節では量子力学への応用として一次元の場合を報告する。調和振動子、非調和振動子、Morseポテンシャル、変形 Poschl -Teller ポテンシャルの計算結果を示す。

2 量子力学と数値計算の関係

量子力学では、1電子ハミルトニアンは一次元において

で与えられる。ここで第1項は運動エネルギー演算子、mは質量、 はDirac常数、第2項はポテンシャルエネルギー演算子である。量子数nのエネルギー固有値Enと固有関数或いは波動関数YnはSchrodinger方程式

と問題に即した両端での境界条件を満たす。つまり、量子力学の固有値問題は二階線形常微分方程式の境界値問題である。Schrodinger方程式はエネルギー単位 と長さの単位 x = axを導入し、ln = En/E0V(x) = U(x)/E0とスケールすると、Yn(x) = yn(x) に対して無次元化された形に書ける

ここでxlnをそれぞれ改めてxEnと置き、以後エネルギーと長さはスケールされたものとする

波動関数を使って演算子A = A(x, d/dx, (d/dx)2)の行列要素が計算できる

行列要素には、例えば、規格直交積分 <n|n'> = dnn', 座標<n| x |n'>, 運動量 <n| -i (d/dx) |n'>, ポテンシャルエネルギー <n| V(x) |n'>, および運動エネルギー <n| - (d/dx)2 |n'> がある。規格直交積分は、解析的な波動関数を用いた場合には数値積分の精度の検証、補間を用いた場合には被積分関数の精度の検証になる。座標とポテンシャルエネルギーの行列要素は数値積分のもう一つの検証になる。運動量と運動エネルギーの行列要素は1階と2階の数値微分の検証になる。固有値とハミルトニアンの期待値との関係は行列要素の精度の更なる検証になる

これは一種のセルフコンシステント条件と見なせる。即ち、右辺は固有関数とその2階微分、およびポテンシャルを用いて計算した行列要素であり、それが固有値と一致することを示している。

3 数値計算法

3. 1 補間

数値計算の基本は補間である。補間は数値計算の至る所で使われているので、単純かつ高精度でなければならない。その要請を満たすのは古典的なLagrange補間である。x軸上にメッシュ点を置き、関数y = f(x)がこれらの離散点において数値表で与えられているとする。(n + 1)個の点xk, k = 0, 1, 2, ...,n, ここでxkは小さな値から大きな値に並べられている、でyk = f(xk)が与えられているとする。中間の点xでのf(x)をn次の多項式で近似して計算する[36 - 39]

ここで

であり、n次の多項式 はLagrange 補間係数と言う

剰余項 Rn

で与えられる。ここで f[x0, x1, ..., xn, x] は差分商であり、f(n+1)(x) はx0 <= x <= xn における(n+1) 階微分である。等間隔区間の場合、即ち、h = (xn - x0)/nの場合、打ち切り誤差はh, nf(n+1)(x) に依存する。離散的なデータ点の間隔を十分短くして、データ点が滑らかにつながるようにする。そして補間に使うデータ点を多く取って高次の多項式を当てはめると、その区間の中央近くでは高精度の補間値が得られる。後程示すように、倍精度演算で15桁以上正確である。これまでの常識「低次の補間多項式を当てはめよ」(例えば[24])は正しくない。補間値は精度良く求まる。また、高精度が得られることと、補間におけるRungeの現象[40, 41]とは矛盾しない。Rungeの現象(高次の補間多項式を当てはめて得られる多項式は区間の端の方で真の値から大きくはずれる)は端の方での現象であり、区間の中央付近では当てはまらない。Rungeの現象は教育的な見地から強調されている。しかしながら、実用的には、hを十分小さく取ってデータ点が滑らかにつながるようにすれば、区間の中央付近では高精度の補間値が得られる。Lagrange補間係数を使う方法は、同一点で補間すべき関数が多数ある場合に有用である。何故ならLagrange補間係数は一回計算すれば、同じものを他の関数にも使えるからである。こうした状況は実際問題でしばしば起こるので、有用である。

3. 2 数値微分

任意の点における数値微分は二つのステップで計算できる。第1ステップでは、表のデータ点での微分を計算する。データ点における微分は、Lagrange補間多項式(6)を微分してデータ点の座標を代入した式を用いて計算する。等間隔の区間幅hに対して、xk, k = 0, 1, ..., n,におけるy = f(x)のm階微分は次のように書ける

ここで係数mnAkiと打ち切り誤差mnEkは10次までが[42, 37]に採録されている。二階微分では(n + 1)点の中央での12次と14次の公式が[15]に採録されている。最も精度が高いのは中央の点での微分である。中央でない点における微分公式は例外的(例えば全区間の端の方にある点では中央の点における公式が適用できない場合)に用いられる。微分は次数nと幅hを適切に選ぶことにより高精度で計算できる。等間隔データ点の間隔を十分短く取り、データ点数を多く取って、高次の多項式を使う。一階微分は10次までの公式で高精度である。倍精度演算で14桁以上正確である。2階微分は式(11)でm = 2の公式を直接使うか、あるいは1階微分を2回使う、のどちらでも良い。第2ステップでは、データ点以外のx座標における微分を、データ点における微分からLagrange補間で計算する。Lagrange補間は高精度であるから、データ点の中間にある点においても微分は高精度に計算できる。これまでの常識「数値微分は避けるべき」は正しくない。「数値微分を積極的に使う」ことができる。
微分の計算において、数値微分の誤差は打ち切り誤差 (m!/hm) mnEkと丸め誤差から成る[37, 43, 24]。表の点における1階微分の打ち切り誤差はcn,1hnf(n+1)(x)の形をしており、2階微分の打ち切り誤差はcn,2hvf(n+2)(x)の形をしている、ここでcn,1cn,2は係数でnが大きくなると絶対値が小さくなる。1階微分と2階微分の丸め誤差はそれぞれ1/h と 1/h2に比例する。

3. 3 数値積分

数値積分法には種々の方法がある。データ点での値のみを用いた高精度数値積分のひとつに、中心差分積分公式(central-difference integration formula)[36] がある。離散かつ相異なる(n + 1)個の点xk, k = 0, 1, 2, ...,n, ここでxkは小さな値から大きな値に並べられている、における関数y = f(x)の値をyk = f(xk)とする。ここで中心点をxiとすると、k = i-(n/2), i-(n/2)+1,..., i-1, i, i+1,..., i+(n/2)、但しnは偶数とする。三点間、即ち二区間[xi-1, xi+1]、における積分は以下のように与えられる

ここでBkは定数である。この方法では区間内の三点の関数値に加えて、区間外の関数値も使って積分を計算する。n = 2に対する積分公式はSimpson公式と同じである。n = 2と4の公式は[14]に、n = 6と8は[15]に、n = 10は[16]に採録されている。この公式により、微小区間での高精度数値積分ができる。また、全区間での積分はこれらの公式を繰り返して計算し、累積を取ることで得られる。全区間での積分値も高精度である。倍精度で15桁以上正確である。

3. 4 離散化行列固有値法

以上の基本的な演算に基づいて離散化行列固有値法を扱う。全区間で総数(N + 1) の離散点xj, j = 0, 1, 2, ...,N,を置く。各離散点xjにおける二階微分に対して、(n + 1)個の離散点の中央xiに対する離散化公式を用いることにより、二階常微分方程式は実空間での行列固有値方程式に書き換えられる

ここで


および

ここで上付き添え字tはベクトルの転置を示す。(N + 1)個の固有値と固有関数があるが、物理的に意味があるのは少数個である。量子力学においてハミルトニアンは通常Hermite行列なので、固有値は実数である。固有値と固有関数の精度を上げるには、二階微分を高次の公式で精度良く離散化することと、全区間を適切に選ぶ必要がある。精度の劣化を避けるため全区間は十分広くとり、両端で固有関数の絶対値が十分小さくなるようにすると同時に、数値誤差による劣化を避けるため固有関数の絶対値が小さ過ぎないようにする。経験から、全区間の選び方は以下の通りである。問題で要求される最大量子数を持つ固有値の相対誤差が1.0D-15の周辺にするには、規格化された固有関数の裾の絶対値が1.0D-15 から 1.0D-10の範囲に収まるようにする。
固有値問題の解法には標準的な方法を用いる[23, 24]。三重対角化にはHouseholder変換を用いる。固有値の計算にはBisection法を用いて固有値を確実に得る。固有ベクトルの計算には逆反復法を用い、高精度の中心差分積分で規格化積分を計算する。
式(4)の右辺の行列要素は高精度の数値微分と積分を用いて計算できる。固有値と固有関数は共に二階微分の公式の次数nを上げることで正確な値に収束する。数値計算の経験から、エネルギーの対角行列要素は固有値よりも速く収束する。これら二つの計算量が一致することは、精度の相互チェックに役立つ。
以上により、高精度の固有値と固有関数、あるいは固有ベクトル、が得られる。また、ハミルトニアンの行列要素を高精度で計算できる。これらは後程述べるように倍精度で13−15桁正確であった。
本節の終わりに、相対誤差と絶対誤差につき述べる。固有値と行列要素のような積分値に対しては相対誤差あるいは正確な桁数を示す。規格化された波動関数とそれらの微分のような被積分関数に対しては絶対誤差を示す。何故なら実際的な数値積分では被積分関数の絶対誤差を制御する必要があるから。

3. 5 Shooting法

3. 5. 1 固有値と固有関数の初期値

離散化行列固有値法で得た解をshooting法のための固有値と固有関数の初期値に用いることができる。ここで、二階微分を最低次の三点公式で離散化する、即ち、行列は始めから三重対角である。このおかげで、三重対角化のためのHouseholder変換は不要である。固有値の相対誤差の絶対値は通常1.0D-5より大きいが 、初期値としては十分な精度である。固有ベクトルは固有関数の初期値に用いられる。更にこの固有関数を用いてハミルトニアンの対角行列要素を計算し、固有値の初期値とすることも可能である。

3. 5. 2 初期値問題への書き換え

両端での境界条件を課された境界値問題は初期値問題に書き換えられる。Schrodinger方程式 (3)

において、右辺をxy及びパラメータE = Enの連続かつ滑らかな関数f(x, y, E)、但しyの一階微分を含まない、に置き換えると、以下の一般化された微分方程式を得る

ここでyの一階微分がなくても一般性を失わない。何故なら二階線型常微分方程式では、一階微分を消去した形に書き直すことが常に可能だからである。初期値問題を扱う際、パラメータE を関数fの中に含めた以下の形の方程式を扱う

3. 5. 3 線形多段法

二階常微分方程式(20)の初期値問題を線形多段法(Linear Multistep Method, LMM)、つまり離散的な点のみでの関数値を計算する、を使って解く。微分方程式(20)を差分方程式の形で近似する一般的な定式は以下に示す通りである[25 - 27, 29]

ここでym, m = 0, 1, 2, ...は等間隔な離散点x = xmにおける差分方程式の解で、fm = f(xm, ym)である。kは整数で、係数ambm, m = 0, ..., k,はnに依らない定数である。この公式においてyn+kは既知であるk個の点における値を使って計算される(k-ステップ法)。hは離散点間のステップ幅であり、二階の常微分方程式に対しては2乗の形で差分方程式に入ってくる。線形多段の線形という言葉は、式(21)の右辺がfmの線形結合であることに由来し、fyに関して線形であることを要求しているのではない。
Dahlquist [28, 29]によって始められ、Henrici [25]によって発展させられた最良作用素の構成方法に基づいて、差分方程式(21)の陰公式を具体的に導出した。我々はこれまで知られていた公式に加えて多くの公式を導出したので、表記法としていちいち人名を付けずに、公式の段数と各段数の中で何番目の公式と呼ぶことにする。得られた公式は[16]にまとめて示してある。以上のLMMの公式を用いることにより、常微分方程式の数値解が求まる。
LMMには他に陽公式もある[44, 45]。これらも高精度を出す[16]。陰公式、陽公式いずれもLMMは単純でかつ高精度を与える有用な手法である。

3. 5. 4 LMM法の線形常微分方程式への適用

線形常微分方程式では、f(x, y)はyの一次関数であり、座標xに依存する係数F(x)とG(x)を用いて次のように表わされる

式(22)を式(21)に代入することで線形多段法による解の表式を得る

3. 5. 5 固有値の修正法

接続点での固有値の修正法を以下に示す。方程式(18)に対して、パラメータEDEずらしたE' = E+ DE に対する微分方程式の解をy' = y + Dy とする

ずれ DE を決める式は以下のように求まる[16]。(18)にy' をかけたものから(24)に y をかけたものを引き、全区間の両端から接続点xmまで座標につき積分し、D [ (dy/dx)/y ] の一次までの和が接続点の両側で等しくなるように置く。更に、分母の積分において積yy'のyDyを無視してy2と置くと次式を得る(Ridleyの公式)[32 - 35]

ここで、xm-は-¥から出発する微分方程式の解を使って計算したxmでの量を意味し、xm+は+¥から解いてきたxmでの値を意味する。式(25)から以下のことが分かる:接続点を固有関数のゼロ点、あるいは節点、の近傍に置いてはならない、何故ならゼロ割と精度の低下を避けるためである。更に、接続点を固有関数の極値の近傍に置いてはならない、何故なら分子にある一階微分の差(固有関数で規格化されている)から生じる桁落ちを避けるためである。以上のことから、接続点はゼロ点と極値となる点の中間に置かねばならない。こうなる点の候補のなかで、全区間の中央付近にある点が最も良い。何故なら、初期値問題を解いて行くと誤差が積もるので、左端から出発した累積誤差と右端から出発した累積誤差は全区間の中央付近でバランス良く小さくなるからである。古典的転換点、それは式(18)でV(x) - En = 0の最小或いは最大のxである[7, 8]、が接続点にしばしば使われてきたが、その理由は固有値の初期値が確実に分かっていないので、固有値が正しい範囲にあるかチェックするため、最小と最大の古典的転換点の間でノード数を数える必要があったためである。しかしながら我々の方法では、固有値と固有関数の初期値を離散化固有値法で確実に得ることができる。従ってノードカウントは不要で、接続点は全区間のどこに置いても良いことから、中央付近に選べるというメリットがある。
更に、(25)の数値積分は高精度の中心差分積分法を用いて計算する。これまで展開してきた公式を用いて、固有値の補正計算をDE の絶対値が許容誤差範囲内に収束するまで繰り返す。

3. 5. 6 他の最高水準の方法との比較

Schrodinger方程式の解法に対する我々の古典的な方法と最近の方法とを比較する。
最初のカテゴリは、コードSLEIGN [46]、SLEIGN2 [47, 48], SLEDGE [49], SL02F [50 - 52], とSLTSTPAK [53] である。最新のコードSLEIGN2 [48]では、Schrodinger方程式はStrum-Liouville (SL)方程式

の特別な形として扱われる。ここでp(x), q(x), r(x)は係数関数、lはパラメータである。2階のSL方程式はPrufer変換

を用いることにより、位相qと振幅rに対する二つの1階非線形微分方程式に書かれる

固有値はqとその微分 に関する非線形微分方程式に対してshooting法を適用して計算される。固有値の初期値はWKB近似に類似の式(前期量子論の量子化条件に対応する式)で計算される。qに対する初期値問題は、誤差評価の補外手順により決められたステップ幅調整付きRunge-Kutta-Fehlberg法[54]を用いて解かれる。固有関数は、振幅rに対する非線形微分方程式を解き、Prufer変換を適用して得られる。SLEIGN2 コードは以下のように実際的な目的に不満足である:SL方程式の係数関数は解析的に与えられなければならない。qrに対する非線形微分方程式はsin q と cos q の関数を含むため、元の線形微分方程式(Schrodinger方程式)よりも計算が複雑である。ステップ幅調整付きRunge-Kutta-Fehlberg法の初期値問題に対する計算時間は線形多段法よりも長い。従って端的に言えば、Schrodinger方程式を線形多段法で直接解くことはSLEIGN2で採用されたアルゴリズムより遥かに速い。SLEIGN2における固有値の収束は遅く、時々間違った値に収束する。固有関数に対する規格化積分の精度は十分でない。更に、行列要素を計算するツールは提供されていない。SLEDGEコード [49]では、SL方程式の係数関数は区分的に一定値で近似され、その結果の方程式は三角及び指数関数を用いて解析的に解かれる。固有値の誤差は正規型に対してO(h2)、ここでhは最大のメッシュサイズ、一般のSL型に対してO(h)であるから、高精度の固有値を得るにはhを十分小さくしなければならない。メッシュ点は区間の中点での二分割を繰り返して追加されるので、全計算が終了した後で得られる。SL02Fコード [50 - 52]では、スケールされたPrufer変換がSL方程式に適用され、係数関数が区分的に一定値で近似され、その結果の方程式は三角および双曲型関数を用いて解析的に解かれる。これは微分方程式のスティッフネスを避けるためである。SL02Fコードを用いた計算時間はSLEIGN2 コードを用いたそれより同程度あるいは1桁短い[51, 55 - 57]。SLEDGE [49], SL02F [50 - 52], SLEIGN [46], とSLEIGN2 [47] はSLTSTPAK [53]に全て含まれている。後で、我々の数値結果と第1のカテゴリを代表するSLEIGN2 [48](但し単精度演算を倍精度演算に置き換えて計算した)の結果を比較する。
第2のカテゴリはCPM(Constant reference potential Perturbation Method)である。CPMではSchrodinger方程式の初期値問題は所謂伝達行列アルゴリズムを使って解かれる。全区間は小区間に分割され、各小区間でポテンシャルは一定値と一定値からのずれに分解される。一定ポテンシャルとずれは、それぞれ基準および摂動として扱われる。プロパゲーター(伝達行列の行列要素)は最大次数Qで打ち切られた摂動級数で構成される。実際的な目的から、ポテンシャルは最大次数Nで打ち切られたずれLegendre多項式で近似される。次数{12, 10} [55- 57]、{14, 12}, {16, 14}, {18, 16} [58, 59]のCPMが報告されている。摂動級数における項数は高次摂動で加速的(破局的)に増える。公式の手計算による導出は手に負えない仕事であり、数式処理、例えばMATHEMATICA, MAPLE, 等、が不可欠である。その結果得られる公式は複雑で、実装が大変であり、計算コストが高い。更に、この方法ではポテンシャルが解析的に与えられる必要がある。何故なら、ずれLegendre多項式を用いた展開の係数はGauss-Legendre積分公式を使って計算されるからである。これは手法の適用性を制限する。その例はポテンシャルが数表の形でのみ与えられる場合で、原子構造計算のように物理の問題でしばしば起こる。ユーザは全区間の両端の位置を指定しなければならないが、その選び方の指針は与えられていない。端点を適切に選ばないと、改良のためのヒントを与えることなしに、計算がしばしば落ちたり、あるいは正しくない固有値を与える (注:我々の方法における両端の選び方は3.4節で示されており、それらはCPMでも使うことができる。) 。CPMにおいて行列要素を計算するツールは提供されていない。高次のCPM公式を使った行列要素の計算時間は、低次のCPM公式を使ったそれよりも長いだろう、何故なら、行列要素の計算では、固有値の計算と違って、多くの点で固有関数を計算する必要があるから。後程、我々の計算結果とコードSLCPM12 [56]を使って計算した結果を比較する。何故なら、我々はMATSLISE[58, 59] に必要なSymbolic Math Toolbox付きMATLABソフトウェアを持っていないから (注:数値計算の手法を大規模計算の部品として使う場合、特別なソフトウェアの使用は、汎用性とポータビリティの観点から望ましくない。) 。我々の方法はCPMよりも遥かに単純で実装が容易であり、精度と計算時間でCPMに遜色がない。
我々はプロセッサー266MHz 、メモリ 128MBのパソコンを用いて計算した。

4 数値計算結果

4. 1 調和振動子

新しい計算法の典型的な応用例として、調和振動子を取り上げる。何故なら解の性質がよく知られているからである[1 - 9]。ポテンシャルは座標の2乗に比例するV(x) = x2。ポテンシャルは|x|®¥で無限大となる。束縛状態は無限個ある。Schrodinger方程式は

厳密解はHermite多項式Hn(x)を用いて、次式で与えられる

数値計算精度の検証では、この厳密解からの誤差をチェックした。

4. 1. 1 補間

Figure 1(a)は波動関数に対する補間の絶対誤差である。ここで誤差は絶対値をとっている。量子数nが0から7までの波動関数に対する各区間の中央での補間の誤差で、3次と9次の多項式の場合を示す。9次では、絶対誤差は1.0D-15より小さく、また所々下に降りているのは倍精度で誤差がゼロ、つまり厳密に正確な結果である。3次の補間で絶対誤差は1.0D-7以下であった。誤差のxnに対する依存性は、式(9)にある誤差項の微分f(n + 1)(x)から(即ち、多項式で最大次数の項xν+n+1の存在から)理解できる[15]。Figure 1(b)は量子数(n = 0)と座標(x = 0.4921875)を固定して、h = 1/64, 1/32, 1/16, と1/8に対して絶対誤差を補間多項式の次数の関数としてプロットしたものである。次数を上げると誤差は減少し、かつ区間幅hが小さい程、誤差が速く減少する。また誤差はh = 1/64 と 1/32では9次で既に1.0D-15よりも小さい。nhの関数としての誤差の挙動は、式(9)の因子pn(x)/(n+1)! により、nの減少関数でかつ因子hnを持つことから理解できる。以上の結果から、多項式補間は実用的な問題に対して、単純かつ高精度であることが示された。


Figure 1. Absolute errors of interpolation for (a) wave functions with quantum numbers n = 0-7 at centers of the interval h using Lagrange interpolation of degree 3 and degree 9, and for (b) wave function with n = 0 as a function of the degree of Lagrange interpolation for h = 1/64, 1/32, 1/16, and 1/8 [15]. Absolute values are taken for the errors; Reprinted with permission from IOP Publishing.

4. 1. 2 1階微分

Figure 2(a)に波動関数の一階微分の絶対誤差を示す。ここで誤差は量子数nが0から7までの波動関数に対するメッシュ点における微分の誤差の絶対値であり、h = 1/64である。10次で絶対誤差は1 .0D-14より小さく、殆ど垂直な線で示されるように、所々で厳密な値(倍精度で誤差がゼロ)である。4次の公式で計算した1階微分は絶対誤差が1.0D-6より小さい。誤差のxnへの依存性は打ち切り誤差cn,1hnf(n+1)(x)のf(n+1)(x)から理解できる。Figure 2(b)では、量子数と座標値を固定して(それぞれx = 0.5 と n = 1)、h = 1/64, 1/32, 1/16, と 1/8に対して誤差を次数の関数としてプロットしたものである。hが小さい程、誤差が速く減少する。hが小さく(1/64)、次数が8より大きい場合に誤差が平らになるのは、1/hに比例する丸め誤差のためである。以上の結果は、1階微分の数値計算が実用的に高精度で求まることを示している。


Figure 2. Absolute errors in numerical first derivative of (a) wave functions for quantum numbers n = 0-7 at mesh-points with the interval h using first derivative formulas of degree 4 and degree 10, and of (b) wave function with n = 1 as a function of the degree of first derivative formulas for h = 1/64, 1/32, 1/16, and 1/8 [15]. Absolute values are taken for the errors; Reprinted with permission from IOP Publishing.

4. 1. 3 2階微分

Figure 3(a)に波動関数の二階微分の絶対誤差を示す。量子数nが0から7までの波動関数に対するメッシュ点における微分の誤差の絶対値で、h = 1/64である。10次で2階微分の絶対誤差は1.0D-12より小さく、4次では1.0D-6より小さい。誤差のxnへの依存性は打ち切り誤差cn,2 hn f(n+2)(x)から理解できる。Figure 3(b)は量子数と座標値を固定して(それぞれx= 0.5 と n = 1)、h = 1/64, 1/32, 1/16, と 1/8に対して、式(11)を用いて計算した誤差を次数の関数としてプロットしたものである。二つの大きなh (h = 1/8 と h = 1/16)で誤差は次数が大きくなると単調に減少するあるいは一定になる。しかしながら、二つの小さなh ( h = 1/64とh = 1/32)で誤差が次数8から12まで単調に下がらないが、それ以外で単調である。小さなhに対してはっきり出るこの非単調性は、1/h2に比例する2階微分の丸め誤差に由来する。以上の結果は、2階微分の数値計算が実用的に高精度で求まることを示している。


Figure 3. Absolute errors in numerical second derivative of (a) wave functions for quantum numbers n = 0-7 at mesh-points with the interval h using second derivative formulas of degree 4 and degree 10, and of (b) wave function with n = 1 as a function of the degree of second derivative formulas for h = 1/64, 1/32, 1/16, and 1/8 [15]. Absolute values are taken for the errors; Reprinted with permission from IOP Publishing.

4. 1. 4 積分

Table 1に示すように、積分および被積分関数の計算精度は波動関数の規格直交積分の精度を検証することではっきり見ることができる。第2欄は正確な被積分関数に対して中心差分積分公式を用いて計算した規格直交積分である。8次の中心差分積分公式を用い、h = 1/64で計算した値は15桁の精度を示している。第3欄はGauss-Hermite積分公式を用い、14分点で厳密に計算した被積分関数を使った数値積分の結果であり、15桁の精度を示している。第4欄はGauss-Hermite積分公式を用い、補間を用いて計算した被積分関数を使った14分点の数値積分の結果である。これは積分および補間が共に15桁の精度を持っていることを示している。更に、表には出していないが、量子数n = 32までの規格直交積分を計算して同じ精度が出ていることを確認している。
行列要素の積分では、他にも、座標の冪乗 <n| xk |n'>, k = 1, 2, 3, 4, と1階微分 <n| d/dx |n'> それぞれの積分を計算して15桁の精度を確認しており、また等式 <n| x |n'> = <n| d/dx |n'> が15桁で成り立っていることも確認した。更に、2階微分<n| (d/dx)2 |n'>でも15桁の精度を確認している[15]。

Table 1. Orthonormal integrals <n|n'> with quantum numbers n and n' of the harmonic oscillator [15]. CDIF: Central-difference integration formula with the degree 8. GHE: Gauss-Hermite quadrature rule with exact integrand at 14 abscissas. GHI: Gauss-Hermite quadrature rule with interpolated integrand at 14 abscissas. Reprinted with permission from IOP Publishing.
<n|n'>CDIFGHEGHI
<0|0>1.000 000 000 000 0001.000 000 000 000 0001.000 000 000 000 000
<1|1>1.000 000 000 000 0000.999 999 999 999 9991.000 000 000 000 000
<2|2>1.000 000 000 000 0001.000 000 000 000 0001.000 000 000 000 000
<3|3>1.000 000 000 000 0001.000 000 000 000 0001.000 000 000 000 000
<4|4>0.999 999 999 999 9990.999 999 999 999 9991.000 000 000 000 000
<5|5>1.000 000 000 000 0001.000 000 000 000 0001.000 000 000 000 000
<6|6>1.000 000 000 000 0001.000 000 000 000 0001.000 000 000 000 000
<7|7>1.000 000 000 000 0001.000 000 000 000 0001.000 000 000 000 000
<8|8>0.999 999 999 999 9991.000 000 000 000 0001.000 000 000 000 000
<9|9>0.999 999 999 999 9990.999 999 999 999 9981.000 000 000 000 000
Others0.000 000 000 000 0000.000 000 000 000 0000.000 000 000 000 000

4. 1. 5 離散化行列固有値法

離散化行列固有値法の誤差をFigure 4に示す。Figure 4(a)は行列の固有値とハミルトニアン行列要素の相対誤差の絶対値で、横軸は二階微分公式の次数である。固有値の相対誤差は次数と共に単調に減少し、次数が12以上で5.0D-13未満に収束する。相対誤差はnが高くなると増大する。行列要素Enの相対誤差は、次数が低い場合、固有値よりも誤差が小さくて収束も速く、次数が6で収束している。いずれも誤差は5.0D-13未満に収まっている。Figure 4(b)は規格化された固有関数の絶対誤差の最大値を次数に対してプロットしたものである。誤差は絶対値を取っている。規格化された固有関数の絶対誤差の最大値は次数を上げると単調減少し、10次以上で5 .0D-13未満に収束した(Figure 4(b))。


Figure 4. (a) Relative errors in eigenvalues for quantum numbers n = 0-9 as a function of degree n of the second derivative formulas and relative errors in matrix elements En [15]. DM indicates the discretized matrix eigenvalue method and ME denotes the matrix elements. (b) Absolute errors in wave functions for n = 0-9 as a function of degree n of the second derivative formulas. Absolute values are taken for the errors; Reprinted with permission from IOP Publishing.

Table 2は調和振動子の結果である。全区間は(-10.0, 10.0)、h = 1/32である。第2欄の固有値は13から14桁正確で、行列要素Enは13から15桁正確で、固有値と同程度の正確さであった。規格化された固有関数の誤差は5.0D-13未満であった。離散化行列固有値法でここまで高精度の結果は世界で始めてであった。(ディフェクト補正を入れたNumerovの離散化行列固有値法ではnが高くなると精度が急速に低下した[60]。)ここで、離散化行列固有値法の固有値とハミルトニアンの対角行列要素の収束性を調べることは、固有値に対する良い相互チェックとなることを注意しておく。

Table 2. Eigenvalues of the discretized matrix equation and matrix elements En = <n| - (d/dx)2 + V(x) |n> / <n|n> with quantum number n of the harmonic oscillator by using central-difference integration formula [15]. The fifth column is the maximum of the absolute value of the absolute error for the wave functions (MAVAEWF); Reprinted with permission from IOP Publishing.
Quantum No.EigenvalueEnExactMAVAEWF
01.00000000000025D+001.00000000000000D+001.01.512D-13
13.00000000000027D+003.00000000000000D+003.02.476D-13
24.99999999999982D+004.99999999999999D+005.03.552D-13
37.00000000000025D+006.99999999999996D+007.03.413D-13
49.00000000000025D+008.99999999999987D+009.04.165D-13
51.10000000000002D+011.09999999999996D+0111.02.358D-13
61.30000000000002D+011.29999999999991D+0113.02.212D-13
71.49999999999993D+011.49999999999980D+0115.03.189D-13
81.69999999999998D+011.69999999999959D+0117.04.456D-13
91.90000000000004D+011.89999999999922D+0119.02.548D-13

4. 1. 6 線形多段法

Figure 5に正確な固有値をパラメータに用い、正確な固有関数を初期値に用いた、線形多段法による初期値問題の数値解の絶対誤差を示す。縦軸は誤差の絶対値である。量子数nは0から7までをとった。波動関数の絶対値が十分小さな値である座標値x = 10から出発してx = 0までステップ幅 h = 1/64で解いた。
公式2_1では、LMMの数値解全体の絶対誤差は1.0D-5 以下であった。n > 0 で絶対誤差がxの小さな所で下向きに尖って小さくなっているのは、波動関数のゼロ点(節点)近くを通るからである。また、nが高くなるとxの大きい所で絶対誤差が大きくなる。何故なら波動関数の解自身が多項式を含むため、有限の大きさを持つ範囲が高いnで広がるからである。この公式2_1(Numerov法あるいはCowell法)は広く使われているが、その精度が低いことは明らかである。そこで高段の公式における精度を順次述べる。4段の公式における絶対誤差は似通っている。誤差は量子数n = 0と座標x = 0 で最大あり、その大きさは公式4_1で2.4D-8である。4段の公式の絶対誤差は2段の公式のそれに比較して2桁小さい。これは公式の段数を上げると数値解の精度が高まることを示唆している。6段の公式では、絶対誤差が量子数n = 0と座標x = 0で最大となるのは4段の公式を使った場合と同様であるが、一方最大誤差は1D-10の程度であり、4段の公式から更に2桁誤差が小さい。6段の公式のなかで最も精度が高い公式6_2の最大誤差は8.7D-11である。8段の公式では、絶対誤差は同じ段数の公式の中では殆ど同じで、最大誤差は1D-12の程度であり、6段の公式から更に2桁精度が向上した。最も精度が高かったのは公式8_1で、最大誤差は9.1D-13であった。10段の公式10_4では、倍精度計算で絶対誤差がしばしばゼロになり、そこでは図で誤差が真下に降りている。また、低い量子数では最大誤差が1D-14の程度であり、8段の公式から更に2桁精度が向上した。公式10_4は最も精度が高く、最大誤差はn = 7で7D-14であった。


Figure 5. Absolute errors in the numerical solutions of the initial-value problems for the ordinary differential equations (ODEs) using the linear multistep method (LMM) for the harmonic oscillator [16]. The ordinate is the absolute value of the errors and the abscissa is coordinate x. The quantum state is taken from n = 0-7. Exact eigenvalue and exact eigenfunction are taken for the initial values. The ODE is solved from x = 10 to x = 0 with the step-size h = 1/64 using the implicit formulas (a) 2_1, (b) 4_1, (f) 6_2, (k) 8_1, (v) 10_4; Reprinted with permission from Elsevier.

LMMの陽公式[44, 45]を用いて同様の計算を行い、高次の公式で高精度解が得られることも示している[16]。陰公式と陽公式の多様な選択肢があり、かつその単純さと高精度の結果からLMMは有用なツールである。

4. 1. 7 Shooting法

shooting法を用いた固有値の収束性を示す。h = 1/32の場合を述べる。離散化固有値問題法で得た固有値を初期値に用いた場合、初期値の誤差は1.0D-4の程度であるが、収束のiteration回数は3〜4回であった。固有値の初期値として、ハミルトニアンの行列対角要素を用いた場合、収束は2〜3回のiterationで一定値に近づいた。
Figure 6(a)はshooting法を用いた収束後の固有値の相対誤差とハミルトニアン行列対角要素の相対誤差を、ステップ幅h = 1/32 と種々のnに対して、段数の関数としてプロットしたものである。固有値の相対誤差は低い段数に対しては段数の関数として単調に下がり、10段で一定値に落ち着く。固有値の相対誤差はnが高くなると大きくなる。一方、行列対角要素のそれは段数が変わってもほぼ一定であり、かつ1D-15近傍と小さいことが分かる。Figure 6(b)は固有関数の絶対誤差の最大値を、ステップ幅h = 1/32 と種々のnに対して、段数の関数としてプロットしたものである。固有関数の最大誤差は段数を上げると単調に減少し、かつ一定値に落ち着く。最大誤差はn が高くなると増加する。


Figure 6. Dependence of the errors in the eigenvalue and eigenfunction on the step-number in the shooting method [16]. The formula with 2-step is represented by that with 2_1, 4-step by 4_1, 6-step by 6_2, 8-step by 8_1, and 10-step by 10_4. The quantum state is shown for n = 0, 3, 6, and 9. (a) The ordinate is the relative error in the eigenvalue and the diagonal Hamiltonian matrix element, the abscissa is the step-number of the formulas with h = 1/32, (b) the ordinate is the largest absolute error in the eigenfunction and the abscissa is the step-number of the formulas with h = 1/32; Reprinted with permission from Elsevier.

Table 3に調和振動子の結果を示す。左から量子数、10_4公式を用いた場合の行列の固有値、ハミルトニアンの対角行列要素、正確な値、固有関数の絶対誤差の最大値、である。全区間(-10.0, 10.0)、区間幅はh = 1/32である。固有値、行列要素、共に、13から15桁正確であった。波動関数の絶対誤差は5.0D-13未満である。Shooting法で、離散化固有値問題法の結果に匹敵する高精度を得た。更に、SLEIGN2で計算した固有値の正確な桁数は12から14桁であり、我々の方法よりも1桁精度が低い。SLEIGN2の計算時間は固有値だけで我々の全計算時間(固有値、固有関数、行列要素)の183倍、即ち2桁、長い。SLCPM12で計算したそれは11から13桁であり、我々の方法よりも2桁精度が低い。SLCPM12の計算時間(固有値と固有関数)は我々の全計算時間(固有値、固有関数、行列要素)と同程度であった。

Table 3. Comparison of exact and computational results using the shooting method presented here and other methods for the harmonic oscillator with V(x) = x2 [16]. The formula 10_4 with h = 1/32 is used and n is the quantum number. Matrix elements of the wave function En = <n| - (d/dx)2 + V(x)|n> /<n|n>; LAVAEWF is the largest absolute value of the absolute (not relative) error in the wave functions. The underlining indicates erroneous digits; Reprinted with permission from Elsevier.
Quantum no.EigenvalueEnExactLAVAEWFOther methods
nSLEIGN2 (a)SLCPM12 (b)
09.99999999999950D-019.99999999999999D-011.01.055D-149.99999999999693D-011.00000000008220D+00
12.99999999999993D+003.00000000000000D+003.07.161D-153.00000000000071D+003.00000000005223D+00
25.00000000000005D+004.99999999999998D+005.01.749D-145.00000000000113D+004.99999999998599D+00
37.00000000000015D+006.99999999999996D+007.06.745D-146.99999999999918D+007.00000000005347D+00
49.00000000000004D+008.99999999999984D+009.05.640D-149.00000000000215D+008.99999999997320D+00
51.10000000000000D+011.09999999999996D+0111.08.732D-141.10000000000054D+011.09999999999966D+01
61.29999999999998D+011.29999999999992D+0113.09.556D-141.30000000000005D+011.30000000000202D+01
71.50000000000001D+011.49999999999979D+0115.05.054D-141.50000000000002D+011.49999999999268D+01
81.70000000000004D+011.69999999999955D+0117.02.259D-131.70000000000022D+011.70000000000204D+01
91.90000000000003D+011.89999999999919D+0119.02.769D-131.89999999999939D+011.89999999999413D+01
CPU-time (seconds)0.2604 (c)47.84 (d)0.2003 (e)
(a) The whole interval is (-¥, +¥) and the boundary conditions are limit points.
(b) The whole interval is the same one used in the present method.
(c) CPU-time for eigenvalues, eigenfunctions and matrix elements.
(d) CPU-time for eigenvalues.
(e) CPU-time for eigenvalues and eigenfunctions.

4. 2 非調和振動子 V(x) = mx2 + l x4

第2の例は非調和振動子V(x) = m x2 + l x4, mlはパラメータ、である[19, 61 - 66]。ポテンシャルはm>=0以上で単一極小、m < 0で二つの極小を持つ。ポテンシャルは|x|®¥で無限大となる。束縛状態は無限個ある。束縛状態の固有値は種々の目的で調べられており、高精度の固有値が他の方法で数値的に求められている[61 - 66]。しかしながら、その解法は特殊で、限られた範囲のmlにしか適用できない、他の問題には使えない、等の限界がある。離散化行列固有値法では、3つの代表的なパラメータの組(m, l) = (0.0, 1.0), (1.0, 1.0), 及び(-1.0, 1.0)に対して統一的に固有値と行列要素を計算し、13から15桁の精度を得た。他の方法がある場合、それらとも13から15桁の精度で一致している[15]。更に、Shooting法を用いて、上記の3つの代表的なパラメータの組(m, l)に対して、固有値と行列要素を計算した。Shooting法の結果は離散化固有値問題法を用いて得た高精度の結果と13から15桁一致している[16]。これら二つの独立した方法でそれぞれ統一的に得た数値が一致していることは、我々の計算の信頼性を示している。これまでの得られた最も精度が高い数値と13から15桁一致している。更に、SLEIGN2で計算した固有値の正確な桁数は12から14桁であり、SLCPM12で計算したそれは10から12桁であった。

4. 3 非調和振動子 V(x) = x2+lx2/(1+g x2)

第3の例は調和振動子にローレンツ型のポテンシャルを重ねあわせた非調和振動子V(x) = x2+lx2/(1+g x2), ここでlgはパラメータ、である[19, 67 - 70]。ポテンシャルは|x|®¥で無限大となる。束縛状態は無限個ある。このポテンシャルはレーザー発振のモデルに使われている。パラメータl, gnの特殊な組み合わせに対して厳密解が得られている。いくつかのパラメータで数値解が得られているが、その解法は特殊で他の問題には適用しにくい。離散化行列固有値法で、厳密解が得られる4つの典型的な場合を計算し、固有値が13桁、ハミルトニアンの対角行列要素が15桁正しいという結果を得た[15]。また、典型的な場合、l = g = 1.0、に対して10個の固有値と行列要素を計算した。Shooting法を用いて、上と同様のパラメータの組に対して計算した結果は離散化固有値問題法を用いて得た高精度の結果と13から15桁一致している[16]。パラメータlgの特殊な組み合わせに対して他の方法で得られた厳密解[19, 67 - 70]と13桁以上一致した。他の方法で結果が出ていない場合でも得られた結果は13桁以上正しいと思われる。更に、SLEIGN2で計算した固有値の正確な桁数は11から14桁であり、SLCPM12で計算したそれは10から13桁であった。

4. 4 Morse ポテンシャル

第4の例はMorseポテンシャルV(x) = V0 (e-2αx - 2 e-αx)である[7, 8]。このポテンシャルはx = 0に関して非対称で、x > 0に対して-V0とゼロの間の有限な範囲をとり、x ® -¥で無限大となる。束縛状態の固有値は有限個である。固有値は解析的に求まり、分子振動で非調和性を解析するのに良く用いられる。離散化行列固有値法を用いて、V0 = 1.0, 2.25, 6.25, 12.25に対して計算した。固有値とハミルトニアンの対角行列要素は13から15桁の精度であった[15]。Shooting法を用いて、上と同様のV0 に対して計算した結果は厳密解に対して13から15桁一致している[16]。SLEIGN2を用いた場合、一般に11から15桁正確であったが、一つ固有値が完全に誤っている場合があった。SLCPM12では10から13桁の精度であった。
Morseポテンシャルのもう一つの型はV(x) = D{1 - exp[-a(x - x0)]}2 [6]、D = we2/4wexe, a = (kwexe)1/2及びk = 1である[20, 71]。ここでx0 = 2.40873, we = 48.66888 およびwexe = 0.977888の場合[20, 71]、離散化行列固有値法で計算した結果、固有値とハミルトニアンの対角行列要素は13から15桁の精度があった[15]。Shooting法で計算した結果は、厳密解に対して固有値が13から15桁一致しており、離散化行列固有値法の結果に匹敵する結果である[16]。SLEIGN2を用いた場合、多くの場合で固有値が誤っていたが、いくつかの場合に11から13桁正確であった。SLCPM12では11から13桁の精度であった。

4. 5 変形Poschl-Tellerポテンシャル

第5の例はmodified Poschl-TellerポテンシャルV(x) = -V0 / cosh2(ax)、V0は定数、である[7, 8]。このポテンシャルはx = 0に関して対称で、-V0とゼロの間の有限な範囲をとる。束縛状態の固有値は有限個である。有限深さのポテンシャルでは、波動関数が広がるので、それに対応して全区間を広げる必要がある。離散化行列固有値法では、V0 = 1.0, 2.0, 6.0, 12.0に対して固有値とハミルトニアンの対角行列要素は厳密解と13桁以上一致した[15]。Shooting法では、厳密解に対して、固有値とハミルトニアンの対角行列要素が13-15桁一致しており、離散化行列固有値法に匹敵する結果である[16]。SLEIGN2を用いた場合、デフォールトの条件では計算できなかったが、全区間を我々の計算に合わせることで得た結果、12から15桁正確であった。SLCPM12では全区間が広すぎても狭すぎても定性的にも正しい解が得られない事態があった。全区間を調整して計算結果が得られた場合、10から13桁の精度であった。

5 結論

一次元の量子力学における固有値問題を高精度で数値的に解く方法を展開した。補間、微分、積分を高精度で計算する手段を提供した。離散化行列固有値法で、高精度の固有値とハミルトニアンの対角行列要素を計算した。常微分方程式の初期値問題を正確かつ安定に解くため、線形多段法の新しい公式を提出し、そのパフォーマンスを示した。また、shooting法により、固有値と固有関数を高精度に解いた。

参考文献

[ 1] V. I. Smirnov, A Course of Higher Mathematics, Vol. 3, Part 2, 6th ed., National Publisher for Technical-Theoretical Literatures, Moscow (1958), (in Russian).
日本語訳、スミルノフ、高等数学教程、第3巻、第2部、共立、東京 (1959), and English translation, Pergamon, Oxford (1964).
[ 2] E. T. Whittaker, G. N. Watson, A Course of Modern Analysis, 4th ed., Cambridge University Press, Cambridge (1927).
[ 3] R. Courant, D. Hilbert, Methods of Mathematical Physics, Vol. 1, Wiley, New York (1955).
[ 4] E. C. Titchmarsh, Eigenfunction Expansions Associated with Second-Order Differential Equations, Vol. 1, 2nd ed., Oxford University Press, Oxford (1962).
[ 5] E. Schrodinger, Collected Papers on Wave Mechanics, 3rd ed., Chelsea, New York (1982).
[ 6] L. Pauling, E. B. Wilson, Jr., Introduction to Quantum Mechanics, Dover, New York (1985).
[ 7] L. D. Landau, E. M. Lifshitz, Quantum Mechanics, 3rd ed., Elsevier, Amsterdam (1963).
[ 8] S. Flugge, Practical Quantum Mechanics, Springer, Berlin (1971).
[ 9] C. S. Johnson, Jr., L. G. Pederson, Problems and Solutions in Quantum Chemistry and Physics, Dover, New York (1986).
[10] M. Abramowitz, I. A. Stegun, ed., Handbook of Mathematical Functions, Dover, New York (1970).
[11] N. N. Lebedev, Special Functions & Their Applications, Dover, New York (1972).
[12] H. Hochstadt, The Functions of Mathematical Physics, Dover, New York (1986).
[13] D. R. Hartree, The Calculation of Atomic Structures, Wiley, New York (1957).
[14] C. Froese Fischer, The Hartree-Fock method for Atoms, Wiley, New York (1977).
[15] H. Ishikawa, J. Phys. A, 35, 4453 (2002).
http://stacks.iop.org/JPhysA/35/4453
[16] H. Ishikawa, J. Comput. Appl. Math. (2006), doi: 10.1016/j.cam.2006.10.035.
[17] L. Pauling, E. B. Wilson, Jr., Introduction to Quantum Mechanics, Dover, New York (1985), p. 202.
[18] G. D. Truhlar, J. Comput. Phys., 10, 123 (1972).
[19] V. Fack, G. Vanden Berghe, J. Phys. A, 18, 3355 (1985).
[20] G. C. Groenenboom, H. M. Buck, J. Chem. Phys., 92, 4374 (1990).
[21] D. R. Hartree, Numerical Analysis, 2nd ed., Oxford University Press, Oxford (1958).
[22] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, 3rd ed., Springer, New York (2002).
[23] M. Mori, Programming for Numerical Calculations with FORTRAN77, (in Japanese).
FORTRAN77数値計算プログラミング, 岩波, 東京 (1987).
[24] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes: The Art of Scientific Computing, 2nd ed., Cambridge University Press, Cambridge (1992).
[25] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations, Wiley, New York (1962).
日本語訳、ヘンリチ、計算機による常微分方程式の解法、離散変数法、サイエンス社、東京(1973)
[26] J. D. Lambert, Computational methods for Ordinary Differential Equations, Wiley, London (1973).
[27] E. Hairer, S. P. Norsett, G. Wanner, Solving Ordinary Differential Equations I, 2nd ed., Springer, Berlin (1993).
[28] G. Dahlquist, Math. Scand., 4, 33 (1956).
[29] G. Dahlquist, Trans. Roy. Inst. Technol., Stockholm, No. 130 (1959).
[30] L. Gr. Ixaru, M. Rizea, Comput. Phys. Commun., 38, 329 (1985).
[31] G. Vanden Berghe, F. Fack, H. E. De Meyer, J. Comput. Appl. Math., 28, 391 (1989).
[32] Ref. [13], Sec. 5.3.1.
[33] Ref. [21], Sec. 7.64.
[34] E. C. C. Ridley, Proc. Camb. Phil. Soc., 51, 702 (1955).
[35] A. S. Douglas, Proc. Camb. Phil. Soc., 52, 636 (1956).
[36] F. B. Hildebrand, Introduction to Numerical Analysis, 2nd ed., Dover, New York (1987).
[37] Z. Kopal, Numerical Analysis, 2nd ed., Chapman-Hall, London (1961).
[38] P. J. Davis, Interpolation and Approximation, Dover, New York (1975).
[39] L. M. Milne-Thomson, The Calculus of Finite Differences, Chelsea, New York (1981).
[40] G. E. Forsyth, M. A. Malcom, C. B. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, Englewood Cliffs, NJ (1977).
[41] M. Mori, Curves and surfaces, (in Japanese).
森正武, 曲線と曲面, 教育出版, 東京 (1974), pp. 30-31.
[42] W. G. Bickley, Math. Gaz, 25, 19 (1941).
[43] S. D. Conte, C. de Boor, Elementary Numerical Analysis, 3rd ed., McGraw-Hill, New York (1981).
[44] G. D. Quinlan, S. Tremaine, Astron. J., 100, 1694 (1990).
[45] G. D. Quinlan, preprint (astro-ph/9901136) 1999.
[46] P. B. Bailey, M. K. Gordon, L. F. Shampine, ACM Trans. Math. Software, 4, 193 (1978).
[47] P. B. Bailey, B. S. Garbow, H. G. Kaper, A. Zettl, ACM Trans. Math. Software, 17, 491 (1991).
[48] P. B. Bailey, W. N. Everitt, A. Zettl, ACM Trans. Math. Software, 27, 143 (2001).
See also SLEIGN2 homepage
http://www.math.niu.edu/~zettl/SL2/
[49] S. Pruess, C. T. Fulton, ACM Trans. Math. Software, 19, 360 (1993).
[50] J. D. Pryce, M. Marletta, Comput. Phys. Commun., 62, 42 (1991).
[51] M. Marletta, J. D. Pryce, J. Comput. Appl. Math., 39, 57 (1992).
[52] J. D. Pryce, Numerical Solution of Strum-Liouville Problems, Oxford University Press, Oxford (1993).
[53] J. D. Pryce, ACM Trans. Math. Software, 25, 21 (1999).
[54] L. F. Shampine, H. A. Watts, ACM Trans. Math. Software, 2, 172 (1976).
[55] L. Gr. Ixaru, H. De Meyer, G. Vanden Berghe, J. Comput. Appl. Math., 88, 289 (1997).
[56] L. Gr. Ixaru, H. De Meyer, G. Vanden Berghe, Comput. Phys. Commun., 118, 259 (1999).
[57] L. Gr. Ixaru, J. Comput. Appl. Math., 125, 347 (2000).
[58] V. Ledoux, M. Van Daele, G. Vanden Berghe, Comput. Phys. Commun, 162, 151 (2004).
[59] V. Ledoux, M. Van Daele, G. Vanden Berghe, ACM Trans. Math. Software, 31, 532 (2005).
See also
http://allserv.ugent.be/~vledoux/MATSLISE/
[60] B. Lindberg, J. Chem. Phys., 88, 3805 (1988).
[61] S. N. Biswas, K. Datta, R. P. Saxena, P. K. Srivastava, V. S. Varma, J. Math. Phys., 14, 1190 (1973).
[62] K. Banerjee, S. P. Bhatnagar, V. Choudhry, S. S. Kanwal, Proc. Royal Soc. London A, 360, 575 (1978).
[63] K. Banerjee, Proc. Royal Soc. London A, 364, 265 (1978).
[64] R. Basla, M. Plo, J. G. Esteve, A. F. Pacheco, Phys. Rev. D, 28, 1945 (1983).
[65] G. Schiffrer, D. Stanzial, Nuovo Cimento, 90B, 74 (1985).
[66] F. N. Fernandez, A. M. Meson, E. A. Castro, J. Phys. A, 18, 1389 (1985).
[67] G. P. Flessas, Phys. Lett., 83A, 121 (1981).
[68] V. S. Varma, J. Phys. A, 14, L489 (1981).
[69] V. Fack, H. de Meyer, G. Vanden Berghe, J. Math. Phys., 27, 1340 (1986).
[70] R. J. W. Hodgson, J. Phys. A, 21, 1563 (1988).
[71] M. Dagher, H. Kobeissi, J. Comput. Chem., 9, 647 (1988).


Return