refineE_schrodinger1d.py 技術ドキュメント

プログラムの動作

refineE_schrodinger1d.py は、1次元シュレーディンガー方程式のエネルギー固有値 \(E\) を、与えられた初期推定値からBrent法を用いて高精度に精密化するPythonプログラムです。

主な機能は以下の通りです。

残差関数の定義: 外部モジュール schrodinger1d の solve_schrodinger 関数を利用して、特定のエネルギー \(E\) に対するシュレーディンガー方程式を解き、計算された波動関数 \(\psi(x)\) の終端 \(x=L\) における値 \(\psi(L)\) を残差として返します。
符号反転区間の自動探索: Brent法が要求する、残差関数 \(\psi(L)\) が符号を反転させるエネルギー区間 \([E_{low}, E_{high}]\) を、与えられた初期推定値 \(E_0\) から自動的に探索します。
Brent法による根の探索: scipy.optimize.brentq 関数を使用して、残差関数 \(\psi(L; E) = 0\) となるエネルギー \(E\) を高精度に計算します。これにより、物理的な境界条件を満たす固有値 \(E\) が特定されます。
結果の出力: 精密化された固有値 \(E\) を標準出力に表示するほか、計算の詳細（初期値、探索区間、最終固有値、反復回数など）を refineE.csv ファイルに追記します。

このプログラムは、1次元シュレーディンガー方程式の固有値問題を数値的に解く際の、高精度な固有値探索という課題を解決します。

原理

このプログラムは、1次元シュレーディンガー方程式の固有値問題を、根探索アルゴリズムであるBrent法に帰着させて解決します。

1次元シュレーディンガー方程式

一般に、1次元の定常状態シュレーディンガー方程式は次のように記述されます。

\[ -\frac{\hbar^2}{2m} \frac{d^2\psi(x)}{dx^2} + V(x)\psi(x) = E\psi(x) \]

ここで、\(\hbar\) はディラック定数、\(m\) は粒子の質量、\(V(x)\) はポテンシャルエネルギー、\(E\) はエネルギー固有値、\(\psi(x)\) は波動関数です。物理的に許容される波動関数は、特定の境界条件（例: \(x \to \infty\) で \(\psi(x) \to 0\)）を満たす必要があり、この条件を満たす \(E\) の値が固有値となります。

残差関数と根の探索

プログラムでは、外部の schrodinger1d モジュールが、与えられたエネルギー \(E\) とポテンシャル \(V(x)\)（プログラム中には明示されていないが、schrodinger1d 内部で定義されているか、デフォルトとして仮定されている）に対してシュレーディンガー方程式を数値的に解き、波動関数 \(\psi(x)\) を計算します。

固有値 \(E\) を探索するために、プログラムは residual(E, args) 関数を定義します。この関数は、solve_schrodinger を呼び出して波動関数 \(\psi(x)\) を計算し、その終点 \(x=L\) における値 \(\psi(L)\) を返します。真の固有値 \(E_{true}\) では、波動関数が物理的な境界条件（例: \(\psi(L)=0\) あるいは漸近的にゼロに近づく挙動）を満たすため、\(\psi(L)\) はゼロになるはずです。

したがって、固有値探索問題は、関数 \(f(E) = \psi(L; E)\) の根 \(E\) を見つける問題に帰着されます。

Brent法

scipy.optimize.brentq は、指定された区間 \([a, b]\) 内で関数 \(f(x)\) の根を探索するロバストなアルゴリズムです。このアルゴリズムは、二分法、逆二次補間、およびセカント法を組み合わせることで、高速かつ確実に根を見つけます。Brent法を使用するには、根を含む区間 \([a, b]\) が与えられ、かつその区間の両端で関数値が符号反転している必要があります（すなわち、\(f(a) \cdot f(b) < 0\)）。

refineE_schrodinger1d.py では、以下の手順でBrent法を適用します。

find_bracket 関数が、初期推定値 \(E_0\) から開始して、\(\psi(L)\) の符号が反転する区間 \([E_{low}, E_{high}]\) を自動的に探索します。main 関数では direction='upper' オプションを指定しているため、初期値 \(E_0\) を下限とし、\(E_0\) より大きい方向に区間を広げていきます（例: \([E_0, E_0 + \delta], [E_0, E_0 + 2\delta], \dots\)）。
compute_maxiter 関数が、目標精度 E_tol と探索区間の幅から、Brent法の最大反復回数を決定します。
brentq 関数が、この区間と残差関数 residual(E) を用いて、目標精度 E_tol で固有値 \(E\) を精密化します。

これにより、効率的かつ高精度にシュレーディンガー方程式の固有値を求めることができます。

必要な非標準ライブラリとインストール方法

このプログラムは、以下の非標準Pythonライブラリを使用します。

numpy: 数値計算（特に配列操作）に使用されます。
scipy: 科学技術計算ライブラリで、特に scipy.optimize.brentq による根の探索に使用されます。
schrodinger1d: 1次元シュレーディンガー方程式を解くためのカスタムモジュールです。

numpy と scipy は pip を使ってインストールできます。

pip install numpy scipy

schrodinger1d は本プログラムのソースコードには含まれていないカスタムモジュールであり、別途用意されている必要があります。このモジュールは、refineE_schrodinger1d.py と同じディレクトリに置かれるか、Pythonのモジュール検索パス（PYTHONPATH）が設定されている場所に存在する必要があります。

必要な入力ファイル

refineE_schrodinger1d.py は、コマンドライン引数を通じて全てのパラメータを受け取るため、直接的な入力ファイルは必要ありません。ただし、以下の前提条件があります。

schrodinger1d.py モジュールの存在: プログラムが schrodinger1d モジュールを import して使用するため、このモジュールが実行環境で利用可能である必要があります。通常は、refineE_schrodinger1d.py と同じディレクトリに schrodinger1d.py ファイルが存在することを想定しています。

生成される出力ファイル

プログラムの実行により、以下のファイルが生成されるか、更新されます。

refineE.csv: 精密化された固有値の計算結果を記録するCSV（Comma Separated Values）ファイルです。
- このファイルが存在しない場合は、ヘッダ行が最初に書き込まれます。
- プログラムが実行されるたびに、新しい計算結果がこのファイルの末尾に追記されます。
- 内容: E0, E_low, E_high, E_root, f_low, f_high, iterations, funcalls
  - E0: 初期推定値
  - E_low: 符号反転区間の下限
  - E_high: 符号反転区間の上限
  - E_root: 精密化された固有値
  - f_low: E_low における残差 \(\psi(L)\) の値
  - f_high: E_high における残差 \(\psi(L)\) の値
  - iterations: Brent法の反復回数
  - funcalls: Brent法による残差関数の呼び出し回数

コマンドラインでの使用例 (Usage)

refineE_schrodinger1d.py は、コマンドライン引数でパラメータを指定して実行します。

python refineE_schrodinger1d.py --E <initial_E> [options]

引数

--E <initial_E> (必須): 固有値の初期推定値。float型。
--E_tol <tolerance> (オプション): 固有値の目標精度。デフォルトは 1e-6。float型。
--maxiter <iterations> (オプション): Brent法の最大反復回数。指定しない場合、プログラムが自動計算します。int型。
--L <length> (オプション): 空間領域の長さ \(L\)。schrodinger1d モジュールに渡されます。デフォルトは 10.0。float型。
--nx <grid_points> (オプション): 空間グリッドの点数。schrodinger1d モジュールに渡されます。デフォルトは 2000。int型。
--psi_max <max_psi> (オプション): 波動関数の最大値（発散チェック用）。schrodinger1d モジュールに渡されます。デフォルトは 1e10。float型。
--bc <boundary_condition> (オプション): 境界条件。"asymptotic" または "zero"。schrodinger1d モジュール内部で処理される境界条件です。デフォルトは "asymptotic"。
--eps <epsilon> (オプション): schrodinger1d モジュール内部で使用される可能性のある許容誤差などの微小量。デフォルトは 1e-6。float型。

コマンドラインでの具体的な使用例

ここでは、初期推定値 \(E=1.0\) を与え、計算領域の長さ \(L=10.0\)、グリッド点数 \(nx=2000\)、目標精度 \(E_{tol}=10^{-8}\) で固有値を精密化する例を示します。

python refineE_schrodinger1d.py --E 1.0 --L 10.0 --nx 2000 --E_tol 1e-8

実行例と結果の説明

上記のコマンドを実行すると、以下のような出力が標準出力に表示され、refineE.csv ファイルが更新されます。

標準出力の例:

[INFO] Found bracket: [1.000000, 1.060000]
[INFO] maxiter set to 40
E= 1.00000000  psi(L)= -1.234567e-01
E= 1.06000000  psi(L)=  9.876543e-02
E= 1.00006000  psi(L)= -1.234567e-03
... (Brent法の探索過程が表示されます) ...
E= 0.99999999  psi(L)= -1.234567e-09
E= 1.00000000  psi(L)=  8.765432e-10
====================================
 Refined eigenvalue:
   E = 1.000000000000
====================================

refineE.csv does not exist. Add header labels.
Append the result to refineE.csv

出力の説明:

[INFO] Found bracket: [...]: 初期推定値 E=1.0 から、残差関数が符号反転する区間 [1.000000, 1.060000] が見つかったことを示します。E0=1.0 で delta = 0.05 * 1.0 + 0.01 = 0.06 と計算され、direction='upper' のため [E0, E0+delta] が最初の区間候補となります。
[INFO] maxiter set to 40: 目標精度 E_tol=1e-8 と探索区間の幅から、Brent法の最大反復回数が自動的に 40 に設定されたことを示します。
E=... psi(L)=...: Brent法の各ステップで試行されたエネルギー \(E\) と、それに対応する波動関数の終端値 \(\psi(L)\) が表示されます。\(\psi(L)\) の値が徐々にゼロに近づいていくのがわかります。
Refined eigenvalue: E = 1.000000000000: 精密化された最終的な固有値が表示されます。
refineE.csv does not exist. Add header labels.: refineE.csv がまだ存在しない場合、ヘッダ行が書き込まれることを示します。
Append the result to refineE.csv: 計算結果が refineE.csv ファイルに追記されたことを示します。

refineE.csv ファイルの内容例:

初回実行時:

E0, E_low, E_high, E_root, f_low, f_high, iterations, funcalls
1.0,1.0,1.06,1.000000000000,-0.1234567,0.0987654,10,12

（上記の値は、schrodinger1d モジュール内のポテンシャルや具体的な境界条件、実装によって異なります。）

この出力により、指定された条件の下でシュレーディンガー方程式の固有値が精密に決定され、その過程が記録されたことが確認できます。