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次元の定常状態シュレーディンガー方程式は次のように記述されます。
ここで、\(\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' を使ってインストールできます。
'''bash 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' は、コマンドライン引数でパラメータを指定して実行します。
'''bash python refineE_schrodinger1d.py --E <initial_E> [options] '''
引数
'--E <initial_E>' (必須): 固有値の初期推定値。'float'型。
'--E_tol
' (オプション): 固有値の目標精度。デフォルトは '1e-6'。'float'型。 '--maxiter
' (オプション): Brent法の最大反復回数。指定しない場合、プログラムが自動計算します。'int'型。 '--L
' (オプション): 空間領域の長さ \(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
' (オプション): 'schrodinger1d' モジュール内部で使用される可能性のある許容誤差などの微小量。デフォルトは '1e-6'。'float'型。
コマンドラインでの具体的な使用例
ここでは、初期推定値 \(E=1.0\) を与え、計算領域の長さ \(L=10.0\)、グリッド点数 \(nx=2000\)、目標精度 \(E_{tol}=10^{-8}\) で固有値を精密化する例を示します。
'''bash 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' ファイルの内容例:
初回実行時:
'''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' モジュール内のポテンシャルや具体的な境界条件、実装によって異なります。）
この出力により、指定された条件の下でシュレーディンガー方程式の固有値が精密に決定され、その過程が記録されたことが確認できます。