はじめに
この記事は，物理学アドベントカレンダー2024 Advent Calendar 2024 24日目の記事です。この記事では，シュレーディンガー（Schrödinger）方程式の最も簡単な場合である、無限に深い井戸型ポテンシャル中のシュレーディンガー方程式を、変分法で数値的に解く場合を解説してみたいと思います（数値計算のコードにはC++20を用いることにします）。そしてこの場合、解析解が三角関数として求められることから、試行波動関数にはヤコビ（Jacobi）の楕円関数を用いてみようと思います。なお、この記事の半分以上（C++20のコードも含めて）はChatGPT o1によって書かれています。
シュレーディンガー方程式の変分法による解法
まず、シュレーディンガー方程式の変分法による解法を簡単に説明します。
変分法では、試行波動関数$\psi_{\mathrm{trial}}$を用いて期待値として定義されるエネルギー汎関数
E[\psi_{\mathrm{trial}}]
= \frac{\langle \psi_{\mathrm{trial}} \mid \hat{H} \mid \psi_{\mathrm{trial}} \rangle}
{\langle \psi_{\mathrm{trial}} \mid \psi_{\mathrm{trial}} \rangle}
を考えます。変分法の基本定理により、いかなる試行波動関数$\psi_{\mathrm{trial}}$ についても、このエネルギー期待値は系の基底エネルギーの値以上となることが知られています。
E[\psi_{\mathrm{trial}}] \ge E_{\mathrm{ground}}
したがって、試行波動関数を適切にパラメータ化し、そのパラメータを最適化 (variation)することで、実際の基底状態エネルギーを上から近似的に求めることが可能となります。試行波動関数は完全系であればいかなる関数でも構いませんが、物理的に適切な関数を選べばエネルギー期待値はより実際の基底エネルギーに近くなります。具体的には、例えば試行波動関数$\psi_{\mathrm{trial}}$を適当な関数$\chi(x)$の$N$個の線形結合として表し、
\left\vert \psi_{trial}\right\rangle =%
%TCIMACRO{\dsum \limits_{n=1}^{N}}%
%BeginExpansion
{\displaystyle\sum\limits_{n=1}^{N}}
%EndExpansion
C_{n}\left\vert \chi_{n}\right\rangle =%
%TCIMACRO{\dsum \limits_{n=1}^{N}}%
%BeginExpansion
{\displaystyle\sum\limits_{n=1}^{N}}
%EndExpansion
C_{n}\chi_{n}\left( x\right)
とすると、エネルギー汎関数は、
E=\dfrac{%
%TCIMACRO{\dsum \limits_{p,q=1}^{N}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{N}}
%EndExpansion
C_{p}^{\ast}C_{q}H_{pq}}{%
%TCIMACRO{\dsum \limits_{p,q=1}^{N}}%
%BeginExpansion
{\displaystyle\sum\limits_{p,q=1}^{N}}
%EndExpansion
C_{p}^{\ast}C_{q}S_{pq}}
となります。ここで、
H_{pq}=\left\langle \chi_{p}\right\vert H\left\vert \chi_{q}\right\rangle
S_{pq}=\left\langle \chi_{p}\right\vert \left. \chi_{q}\right\rangle
です。ここでは前者をハミルトニアン積分（あるいは行列と見てハミルトニアン行列）、後者を重なり積分（あるいは行列と見て重なり行列）と呼ぶことにします。ここで、パラメータ$C_n$の最適化を行うとすると、エネルギー汎関数の$C_n$に関する導関数がゼロになるという条件から、
%
%TCIMACRO{\dsum \limits_{q=1}^{N}}%
%BeginExpansion
{\displaystyle\sum\limits_{q=1}^{N}}
%EndExpansion
\left( H_{pq}-ES_{pq}\right) C_{q}=0,\left( p=1,\ldots,N\right)
となります。これは結局、一般化固有値問題
\mathbf{HC}=E\mathbf{SC}
を解くことに帰着します。
無限に深い井戸型ポテンシャル中のシュレーディンガー方程式のヤコビの楕円関数による展開
以下では、1次元の無限に深い井戸型ポテンシャル
V(x) =
\begin{cases}
0 & (0 \le x \le 1),\\
\infty & (x < 0 \ \text{または}\ x > 1)
\end{cases}
での定常シュレーディンガー方程式
-\dfrac{\hbar^{2}}{2M}\dfrac{d^{2}\psi\left( x\right) }{dx^{2}}+V\left(
x\right) =E\psi\left( x\right)
を、ヤコビの楕円関数 $ \mathrm{sn}(u,m) $ を用いた基底関数展開で扱う場合の、ハミルトニアン行列
H_{ij} \;=\; \langle \chi_i \,\bigl|\bigr.\, \hat{H} \,\bigl|\bigr.\, \chi_j \rangle
および重なり行列
S_{ij} \;=\; \langle \chi_i \,\bigl|\bigr.\, \chi_j \rangle
の具体的な形を示します。
1. 基底関数χ(r)の構築
無限に深い井戸型ポテンシャル中の波動関数$\psi(x)$では、境界条件
\psi(0) = 0,
\quad
\psi(1) = 0
を満たす必要があります。ヤコビの楕円関数 $ \mathrm{sn}(u,m) $ は
\mathrm{sn}(0,m) = 0
が成り立つので、$x = 0$で消える条件を自然に満たします。さらに$x = 1$でもゼロとなるようにするには、
\mathrm{sn}\bigl(\omega_n \cdot 1,\,m\bigr) = 0
が必要です。
ヤコビの楕円関数 $ \mathrm{sn}(u,m) $ には周期的性質があり、第一種完全楕円積分$ K(m) $を用いると
\mathrm{sn}\bigl(2K(m),m\bigr) = 0,
\quad
\mathrm{sn}\bigl(4K(m),m\bigr) = 0,
\quad
\mathrm{sn}\bigl(6K(m),m\bigr) = 0,
\ldots
と、偶数倍($ 2nK(m) $)で零点を持つことが知られています。したがって
\omega_n = 2\,n\,K(m)
\quad
(n = 1, 2, 3,\dots)
とすれば、
\mathrm{sn}\bigl(\omega_n \cdot 0, m\bigr)
= \mathrm{sn}(0,m) = 0,
\quad
\mathrm{sn}\bigl(\omega_n \cdot 1, m\bigr)
= \mathrm{sn}\bigl(2n\,K(m),m\bigr) = 0,
となり、境界条件を満たします。
よって、たとえば基底関数を
\chi_n(x) \;=\; N_n \, \mathrm{sn}\bigl(\omega_n\,x,\,m\bigr),
\quad
\omega_n = 2\,n\,K(m),
のように定めることが考えられます。ここで$ N_n $は規格化定数です。
2. 重なり行列の構築
基底関数$ \chi_n(x) $を用いて、重なり行列
S_{ij}
\;=\;
\langle \chi_i \mid \chi_j \rangle
\;=\;
\int_{0}^{1} \chi_i(x)\,\chi_j(x)\,dx
を計算すると、
S_{ij}
\;=\;
\int_{0}^{1}
N_i\,N_j \,
\mathrm{sn}(\omega_i x, m)\,\mathrm{sn}(\omega_j x, m)\,
dx.
したがって、
S_{ij}
\;=\;
N_i \, N_j
\int_{0}^{1}
\mathrm{sn}(\omega_i x, m)\,\mathrm{sn}(\omega_j x, m)\,
dx.
となります。一般には、$ \mathrm{sn}(u,m),\mathrm{sn}(v,m) $の積分は楕円積分の形をとり、解析的に単純な閉形式を得るのは難しいので、この記事では数値積分を用いて求めます。
さらに、同じ$ i = j = n $のときの規格化条件
1
=
\langle \chi_n \mid \chi_n \rangle
=
N_n^2
\int_{0}^{1}
\mathrm{sn}^2(\omega_n x, m)\,dx
から$ N_n $が求められます。
3. ハミルトニアン行列の構築
ポテンシャルが無限に深い井戸の内部でのみ$V(x) = 0$となる場合、ハミルトニアンは運動エネルギーの演算子のみ
\hat{H}
=
-\frac{\hbar^2}{2M}\,\frac{d^2}{dx^2}
で与えられます。よって
H_{ij}
=
\int_{0}^{1}
\chi_i(x)
\biggl[
-\frac{\hbar^2}{2M}\,\frac{d^2}{dx^2}
\biggr]
\chi_j(x)\,
dx.
これを部分積分により整理します。一般に
\int_{0}^{1}
f(x)\,\frac{d^2}{dx^2}g(x)\,dx
=
\Bigl[f(x)\,\frac{d}{dx}g(x)\Bigr]_{0}^{1}
-
\int_{0}^{1}
\frac{d}{dx}f(x)\,\frac{d}{dx}g(x)\,dx.
です。無限に深い井戸型ポテンシャルの境界では$ \chi_i(0)=\chi_i(1)=0 $と$ \chi_j(0)=\chi_j(1)=0 $であるため、境界項は消えます。すると、
\int_{0}^{1}
\chi_i(x)\,\frac{d^2}{dx^2}\chi_j(x)\,dx
=
-\int_{0}^{1}
\frac{d}{dx}\chi_i(x)\,\frac{d}{dx}\chi_j(x)\,dx.
となり、したがって
H_{ij}
=
-\frac{\hbar^2}{2M}
\int_{0}^{1}
\chi_i(x)\,\frac{d^2}{dx^2}\chi_j(x)\,dx
=
\frac{\hbar^2}{2M}
\int_{0}^{1}
\frac{d}{dx}\chi_i(x)\,\frac{d}{dx}\chi_j(x)\,dx.
です。ここで、基底関数$ \chi_n(x) = N_n\mathrm{sn}(\omega_n x, m) $を微分すると、ヤコビの楕円関数の性質
\frac{d}{du}\,\mathrm{sn}(u,m)
=
\mathrm{cn}(u,m)\,\mathrm{dn}(u,m)
を用いて、
\frac{d}{dx}\,\mathrm{sn}(\omega_n x, m)
=
\omega_n \,\mathrm{cn}(\omega_n x, m)\,\mathrm{dn}(\omega_n x, m).
よって
\frac{d}{dx}\chi_n(x)
=
N_n \,\omega_n \,\mathrm{cn}(\omega_n x, m)\,\mathrm{dn}(\omega_n x, m).
すると、
\frac{d}{dx}\chi_i(x)\,\frac{d}{dx}\chi_j(x)
=
N_i\,N_j \,\omega_i\,\omega_j \,
\mathrm{cn}(\omega_i x, m)\,\mathrm{dn}(\omega_i x, m)\,
\mathrm{cn}(\omega_j x, m)\,\mathrm{dn}(\omega_j x, m).
となります。ゆえに結局、
H_{ij}
=
\frac{\hbar^2}{2M}
\int_{0}^{1}
\frac{d}{dx}\chi_i(x)\,\frac{d}{dx}\chi_j(x)\,dx
=
\frac{\hbar^2}{2M}\,
N_i\,N_j \,\omega_i\,\omega_j
\int_{0}^{1}
\mathrm{cn}(\omega_i x, m)\,\mathrm{dn}(\omega_i x, m)\,
\mathrm{cn}(\omega_j x, m)\,\mathrm{dn}(\omega_j x, m)\,
dx.
を得ます。ここでも、$ \mathrm{cn}(u,m)\mathrm{dn}(u,m)\mathrm{cn}(v,m)\mathrm{dn}(v,m) $といった形の積分が現れ、数値積分が必要となります。
4. この節のまとめ
以上を要約すると、無限に深い井戸型ポテンシャル内で$\omega_n = 2nK(m)$と取り、
\chi_n(x) = N_n\,\mathrm{sn}(\omega_n x,m)
を基底関数系として用いるとき、ハミルトニアン行列及び重なり行列は次のように書けます。
-
重なり行列:
S_{ij} = \langle \chi_i|\chi_j\rangle = N_i \, N_j \int_{0}^{1} \mathrm{sn}(\omega_i x, m)\,\mathrm{sn}(\omega_j x, m)\, dx. -
ハミルトニアン行列:
H_{ij} = \langle \chi_i|\hat{H}|\chi_j\rangle = \frac{\hbar^2}{2M}\, N_i\,N_j \,\omega_i\,\omega_j \int_{0}^{1} \mathrm{cn}(\omega_i x, m)\,\mathrm{dn}(\omega_i x, m)\, \mathrm{cn}(\omega_j x, m)\,\mathrm{dn}(\omega_j x, m)\, dx.
これらを有限次元の基底数$N$だけ取り出して、
\mathbf{HC}=E\mathbf{SC}
を解けば、今回の目的が達せられたことになります。
今回の数値計算で用いたC++20のソースコード
今回は、モジュラス$m$を$0.5$に固定し、$N$を可変として以下のC++20のソースコードで数値計算を行いました。ただし、ヤコビの楕円関数の計算とGauss-Legendre積分にBoost C++ Librariesを、一般化固有値問題を解くためにEigenライブラリを用いました。また、コンパイルにはgcc (g++) 13などのC++20に対応したコンパイラが必要な点に注意してください。
#include <cmath> // for std::sin, std::sqrt
#include <cstdint> // for std::int32_t
#include <format> // for std::format
#include <fstream> // for std::ofstream
#include <functional> // for std::function
#include <iostream> // for std::cout, std::endl
#include <numbers> // for std::numbers::pi
#include <vector> // for std::vector
// ---- Boost 関連 (Jacobi elliptic, ellint_1, Gauss–Legendre) ----
#include <boost/math/quadrature/gauss.hpp>
#include <boost/math/special_functions/ellint_1.hpp>
#include <boost/math/special_functions/jacobi_elliptic.hpp>
// ---- Eigen 関連 (generalized eigenvalue problem) ----
#include <Eigen/Dense>
//============================================================
// 物理定数やパラメータ (ここでは簡単のため \hbar=1, 質量=1)
//============================================================
static auto constexpr GAUSS_ORDER = 64;
static auto constexpr HBAR = 1.0;
static auto constexpr MASS = 1.0;
static auto constexpr OUTPUTNUM = 1000;
//============================================================
// Gauss–Legendre 積分 (Boost) を使った数値積分
//============================================================
double numeric_integral(std::function<double(double)> const & func, double x_min, double x_max)
{
boost::math::quadrature::gauss<double, GAUSS_ORDER> integrator;
return integrator.integrate(func, x_min, x_max);
}
//============================================================
// Boost: Jacobi elliptic sn, cn, dn をまとめて返す構造体
//============================================================
struct JacobiElliptic
{
double sn;
double cn;
double dn;
};
//============================================================
// jacobi_elliptic(u, m):
// Boost の jacobi_sn(u, k), jacobi_cn(u, k), jacobi_dn(u, k)
// ただし k^2 = m => k = sqrt(m).
//============================================================
JacobiElliptic jacobi_elliptic(double u, double m)
{
auto const k = std::sqrt(m);
JacobiElliptic res;
res.sn = boost::math::jacobi_sn(k, u);
res.cn = boost::math::jacobi_cn(k, u);
res.dn = boost::math::jacobi_dn(k, u);
return res;
}
//============================================================
// \omega_n = 2*(n+1)*K(m)
// (n は 0-based, n+1 が本来の "モード番号")
//============================================================
double omega_n(std::int32_t n, double m)
{
auto const k = std::sqrt(m);
auto const Kval = boost::math::ellint_1(k); // 完全楕円積分 K(m)
// n=0 -> 2*1*K(m) = 2K(m)
// n=1 -> 2*2*K(m) = 4K(m), など
return 2.0 * static_cast<double>(n + 1) * Kval;
}
//============================================================
// phi_n_raw(x, n, m):
// sn( \omega_n * x, m ) [正規化定数は掛けない]
//============================================================
double phi_n_raw(double x, std::int32_t n, double m)
{
auto const w = omega_n(n, m);
auto const je = jacobi_elliptic(w * x, m);
return je.sn; // sn(...)
}
//============================================================
// 正規化係数 N_n を計算
// N_n = 1 / sqrt( \int_0^1 sn^2(\omega_n x, m) dx )
//============================================================
double normalization_factor(std::int32_t n, double m, std::int32_t gauss_order = 32)
{
auto const integrand = [&](auto xx) {
auto const val = phi_n_raw(xx, n, m);
return val * val;
};
auto const I = numeric_integral(integrand, 0.0, 1.0);
return 1.0 / std::sqrt(I);
}
//============================================================
// (正規化込み) \phi_n(x) = N_n * sn(\omega_n x, m)
//============================================================
double phi_n(double x, std::int32_t n, double m, double Nn)
{
return Nn * phi_n_raw(x, n, m);
}
//============================================================
// S_{ij} = \int_0^1 \phi_i(x) \phi_j(x) dx
//============================================================
double overlap_Sij(std::int32_t i, std::int32_t j, double m, std::vector<double> const & Norm)
{
auto const integrand = [&](auto xx) {
auto const pi_val = phi_n_raw(xx, i, m) * Norm[i];
auto const pj_val = phi_n_raw(xx, j, m) * Norm[j];
return pi_val * pj_val;
};
return numeric_integral(integrand, 0.0, 1.0);
}
//============================================================
// H_{ij} = (hbar^2 / 2m) \int_0^1 \phi_i'(x)* \phi_j'(x) dx
// ( 無限井戸内 V=0 => 運動エネルギー項のみ )
//============================================================
double hamiltonian_Hij(std::int32_t i, std::int32_t j, double m, std::vector<double> const & Norm)
{
auto const phi_n_prime = [&](auto x, auto nIndex) {
// d/dx [ N_n * sn(\omega_n x, m ) ]
// = N_n * [ \omega_n * cn(...) * dn(...) ]
auto const w = omega_n(nIndex, m);
auto const je = jacobi_elliptic(w * x, m);
auto const deriv = w * je.cn * je.dn;
return Norm[nIndex] * deriv;
};
auto integrand = [&](auto const xx) {
auto const dpi = phi_n_prime(xx, i);
auto const dpj = phi_n_prime(xx, j);
return dpi * dpj;
};
auto const val = numeric_integral(integrand, 0.0, 1.0);
auto const factor = (HBAR * HBAR) / (2.0 * MASS);
return factor * val;
}
//============================================================
// 行列 S, H を構築し、一般化固有値問題を解いて
// 固有値と固有ベクトルを得る。
//============================================================
struct DiagResult
{
std::vector<double> norms; // 各基底の正規化定数 [n=0..N-1]
Eigen::VectorXd eigenvalues;
Eigen::MatrixXd eigenvectors; // 各列が固有ベクトル
};
DiagResult solve_eigen_problem(std::int32_t N, double m)
{
DiagResult res;
res.norms.resize(N);
// 1) 各基底の正規化定数
for (auto n = 0; n < N; ++n) {
res.norms[n] = normalization_factor(n, m);
}
// 2) 行列 S, H
Eigen::MatrixXd S = Eigen::MatrixXd::Zero(N, N);
Eigen::MatrixXd H = Eigen::MatrixXd::Zero(N, N);
for (auto i = 0; i < N; i++) {
for (auto j = 0; j < N; j++) {
auto const sij = overlap_Sij(i, j, m, res.norms);
auto const hij = hamiltonian_Hij(i, j, m, res.norms);
S(i, j) = sij;
H(i, j) = hij;
}
}
// 3) 一般化固有値問題 H c = E S c を解く
Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> solver(H, S);
res.eigenvalues = solver.eigenvalues(); // 昇順に並ぶ
res.eigenvectors = solver.eigenvectors(); // (N x N)行列。各列が固有ベクトル
return res;
}
//============================================================
// ヤコビ楕円基底で得た k番目(0-based) の固有状態の波動関数 \psi_k(x)
// psi_k(x) = sum_{n=0..N-1} c_{n,k} * [ N_n * sn(\omega_n x, m) ]
// c_{n,k} = evecs(n, k) (行: n, 列: k)
//============================================================
double wavefunc_k(double x, std::int32_t k, const Eigen::MatrixXd & evecs, std::vector<double> const & norms, double m,
std::int32_t N)
{
auto sum = 0.0;
for (auto n = 0; n < N; n++) {
auto const c_nk = evecs(n, k);
// phi_n(x) = norms[n] * sn(omega_n*x, m)
auto const phi_val = phi_n_raw(x, n, m) * norms[n];
sum += c_nk * phi_val;
}
return sum;
}
//============================================================
// 無限井戸の厳密解 (節の数=ell) に対応する関数
// node = 0 -> \ell=1
// node = 1 -> \ell=2
// ...
// \psi_exact_node(x) = sqrt(2)*sin( (node+1)*pi*x )
//============================================================
double wavefunc_exact(std::int32_t node, double x)
{
// node=0 => l=1
// node=1 => l=2
// ...
auto const l = node + 1;
return std::sqrt(2.0) * std::sin(l * std::numbers::pi * x);
}
//============================================================
// メイン
//============================================================
std::int32_t main()
{
// ---------------------------
// 1) パラメータ設定
// ---------------------------
auto constexpr m = 0.5; // モジュラス m (固定)
auto constexpr N = 15; // 基底数
auto constexpr num_states_to_compare = 5; // 基底状態(節0)～第4励起(節4)
std::cout << std::format("=== Jacobi-sn basis (m={:.3f}, N={}) ===\n", m, N);
// ---------------------------
// 2) 一般化固有値問題を解く
// ---------------------------
auto const diag = solve_eigen_problem(N, m);
// 固有値, 固有ベクトルを取得
auto & Evals = diag.eigenvalues;
auto & Evecs = diag.eigenvectors;
// Evals[k] が k番目(0-based)の固有値 (昇順)
// Evecs.col(k) が k番目の固有ベクトル
// ---------------------------
// 3) 固有値を表示 (0～4まで)
// ---------------------------
std::cout << "\nEigenvalues (lowest 5):\n";
for (auto k = 0; k < num_states_to_compare; k++) {
std::cout << std::format(" E[{}] = {:.6f} E(exact)[{}] = {:.6f}\n", k, Evals[k], k,
std::numbers::pi * std::numbers::pi * 0.5 * static_cast<double>((k + 1) * (k + 1)));
}
// ---------------------------
// 4) 波動関数をファイルに出力して比較
// x=0.0～1.0(0.01刻み) で
// (a) 数値解 psi_k(x) (k=0..4)
// (b) 厳密解
// node=0..4 => psi_exact(x) = sqrt{2} sin((node+1)*pi*x)
// ---------------------------
std::ofstream ofs("wavefunctions_compare.txt");
if (!ofs) {
std::cerr << "Cannot open wavefunctions_compare.txt\n";
return 1;
}
ofs << "# x psi_num[0..4] psi_exact[0..4]\n";
for (auto i = 0; i <= OUTPUTNUM; i++) {
auto const x = static_cast<double>(i) / static_cast<double>(OUTPUTNUM);
ofs << std::format("{:.3f}", x);
// 数値解
for (std::int32_t k = 0; k < num_states_to_compare; ++k) {
auto const val_num = wavefunc_k(x, k, Evecs, diag.norms, m, N);
ofs << std::format(" {:.6f}", val_num);
}
// 厳密解
for (std::int32_t node = 0; node < num_states_to_compare; ++node) {
auto const val_exact = wavefunc_exact(node, x);
ofs << std::format(" {:.6f}", val_exact);
}
ofs << std::endl;
}
std::cout << std::format("\nWavefunctions (k=0..4) and exact solutions have been saved to {}\n",
"wavefunctions_compare.txt");
std::cout << "Done.\n";
return 0;
}
Juliaのソースコード
ついでにJuliaのソースコードもChatGPTに書いてもらいました。以下に示します。
#!/usr/bin/env julia
using Printf
using LinearAlgebra
using FastGaussQuadrature
using Elliptic # sn(u, m), cn(u, m), dn(u, m), elliptic_k(m) など
# ------------------------------------------------------------
# 定数の定義
# ------------------------------------------------------------
const GAUSS_ORDER = 64
const HBAR = 1.0
const MASS = 1.0
const OUTPUTNUM = 1000
# ------------------------------------------------------------
# Gauss–Legendre を用いた数値積分
# ------------------------------------------------------------
function numeric_integral(func::Function, x_min::Float64, x_max::Float64)
xg, wg = gausslegendre(GAUSS_ORDER)
half_range = 0.5 * (x_max - x_min)
mid = 0.5 * (x_max + x_min)
return half_range * sum(wg[i] * func(mid + half_range*xg[i]) for i in eachindex(xg))
end
# ------------------------------------------------------------
# Jacobi 3種函数をまとめて返す構造体
# ------------------------------------------------------------
struct JacobiElliptic
sn::Float64
cn::Float64
dn::Float64
end
function jacobi_elliptic(u::Float64, m::Float64)::JacobiElliptic
JacobiElliptic(Jacobi.sn(u, m), Jacobi.cn(u, m), Jacobi.dn(u, m))
end
# ------------------------------------------------------------
# ωₙ = 2 * n * K(m)
# ------------------------------------------------------------
function omega_n(n::Int, m::Float64)
Kval = Elliptic.K(m)
return 2.0 * n * Kval
end
# ------------------------------------------------------------
# φₙ(x) の生の形 (正規化なし)
# ------------------------------------------------------------
function phi_n_raw(x::Float64, n::Int, m::Float64)
w = omega_n(n, m)
je = jacobi_elliptic(w*x, m)
return je.sn
end
# ------------------------------------------------------------
# 正規化係数 Nₙ
# ------------------------------------------------------------
function normalization_factor(n::Int, m::Float64)
integrand(xx) = (phi_n_raw(xx, n, m))^2
I = numeric_integral(integrand, 0.0, 1.0)
return 1.0 / sqrt(I)
end
# ------------------------------------------------------------
# (正規化込み) φₙ(x) = Nₙ * sn(ωₙ x, m)
# ------------------------------------------------------------
function phi_n(x::Float64, n::Int, m::Float64, Nn::Float64)
return Nn * phi_n_raw(x, n, m)
end
# ------------------------------------------------------------
# Sᵢⱼ = ∫₀¹ φᵢ(x) φⱼ(x) dx
# ------------------------------------------------------------
function overlap_Sij(i::Int, j::Int, m::Float64, norms::Vector{Float64})
integrand(xx) = phi_n_raw(xx, i, m)*norms[i] * phi_n_raw(xx, j, m)*norms[j]
numeric_integral(integrand, 0.0, 1.0)
end
# ------------------------------------------------------------
# Hᵢⱼ = (ℏ² / 2m) ∫₀¹ φᵢ'(x) φⱼ'(x) dx
# ------------------------------------------------------------
function hamiltonian_Hij(i::Int, j::Int, m::Float64, norms::Vector{Float64})
function phi_n_prime(x::Float64, n::Int)
w = omega_n(n, m)
je = jacobi_elliptic(w*x, m)
return norms[n] * w * je.cn * je.dn
end
integrand(xx) = phi_n_prime(xx, i)*phi_n_prime(xx, j)
val = numeric_integral(integrand, 0.0, 1.0)
factor = (HBAR * HBAR) / (2.0 * MASS)
return factor * val
end
# ------------------------------------------------------------
# 一般化固有値問題を解く
# ------------------------------------------------------------
mutable struct DiagResult
norms :: Vector{Float64}
eigenvalues :: Vector{Float64}
eigenvectors:: Matrix{Float64}
end
function solve_eigen_problem(N::Int, m::Float64)
norms = [normalization_factor(n, m) for n in 1:N]
S = zeros(Float64, N, N)
H = zeros(Float64, N, N)
for i in 1:N
for j in 1:N
S[i,j] = overlap_Sij(i, j, m, norms)
H[i,j] = hamiltonian_Hij(i, j, m, norms)
end
end
es = eigen(H, S)
DiagResult(norms, es.values, es.vectors)
end
# ------------------------------------------------------------
# k番目 (1-based) の固有状態 ψₖ(x)
# ------------------------------------------------------------
function wavefunc_k(x::Float64, k::Int, evecs::Matrix{Float64}, norms::Vector{Float64}, m::Float64, N::Int)
s = 0.0
for n in 1:N
c_nk = evecs[n, k]
phi_val = phi_n_raw(x, n, m)*norms[n]
s += c_nk * phi_val
end
return s
end
# ------------------------------------------------------------
# 無限井戸の厳密解 (節数でなくモード番号 k=1,2,...)
# ------------------------------------------------------------
function wavefunc_exact(k::Int, x::Float64)
return sqrt(2.0)*sin(k * π * x)
end
# ------------------------------------------------------------
# main
# ------------------------------------------------------------
function main()
m = 0.5
N = 15
nstates = 5
@printf("=== Jacobi-sn basis (m=%.3f, N=%d) ===\n", m, N)
diag = solve_eigen_problem(N, m)
Evals = diag.eigenvalues
Evecs = diag.eigenvectors
@printf("\nEigenvalues (lowest 5):\n")
for k in 1:nstates
E_exact = (π^2 * 0.5)*(k^2)
@printf(" E[%d] = %.6f E(exact)[%d] = %.6f\n",
k, Evals[k], k, E_exact)
end
filename = "wavefunctions_compare.txt"
open(filename, "w") do io
@printf(io, "# x psi_num[1..5] psi_exact[1..5]\n")
for i in 0:OUTPUTNUM
x = i/OUTPUTNUM
psi_nums = [wavefunc_k(x, k, Evecs, diag.norms, m, N) for k in 1:nstates]
psi_exacts = [wavefunc_exact(k, x) for k in 1:nstates]
@printf(io, "%.3f", x)
for val in psi_nums
@printf(io, " %.6f", val)
end
for val in psi_exacts
@printf(io, " %.6f", val)
end
@printf(io, "\n")
end
end
@printf("\nWavefunctions (k=1..%d) and exact solutions have been saved to %s\n", nstates, filename)
@printf("Done.\n")
end
main()
結果と議論
以下に、固有値（固有エネルギー）の計算結果と厳密値の比較を示します。ここで、粒子の質量$M$は$1$とし、また$\hbar = 1$なる単位系を使いました。ここで、基底状態は波動関数の節の数が$0$に対応しており、第$n$励起状態は節の数が$n$に対応しています。
| 状態 | N = 5 | N = 20 | N = 50 | 厳密 |
|---|---|---|---|---|
| 基底状態 | 4.935957 | 4.934816 | 0.131572 | 4.934802 |
| 第一励起状態 | 20.010264 | 19.740007 | 2.737483 | 19.739209 |
| 第二励起状態 | 45.023034 | 44.423137 | 4.361770 | 44.413320 |
| 第三励起状態 | 80.041058 | 78.974932 | 7.508978 | 78.956835 |
| 第四励起状態 | 125.064217 | 123.433592 | 24.031326 | 123.370055 |
結果を見ると分かるとおり、$N = 5$でも十分な精度を持ち、$N = 20$では、厳密値に非常に近くなっていますが、$N = 50$では極めて過小評価するという結果になっています。$N = 50$の結果は変分原理にも反しており、なぜこうなってしまうのか私には分かりませんでした。
次に、第一励起状態の近似波動関数（$N = 5$、$N = 20$、$N = 50$）と厳密な波動関数の比較のプロットを示します。
このプロットを見ると、$N = 5$でも十分な精度を持っており、$N = 20$だと、プロット上は完全に厳密解に一致することが分かります。$N = 50$の場合は、よく分かりませんが何か特殊なことが起きているようです。
全体のまとめ
本記事では、無限に深い井戸型ポテンシャル中のシュレーディンガー方程式を、ヤコビの楕円関数で展開し、三角関数で表される厳密解と比較しました。使用するヤコビの楕円関数の数（基底の数）$N$が増えても、必ずしも精度が向上しない（それどころか変分原理に反してしまう）という興味深い結果になりました。あとChatGPT o1は本当に賢いです。
参考サイト
[1] ChatGPT
参考文献
[1] J.M.ティッセン 『計算物理学』シュプリンガー・フェアラーク東京（2003）