無限に深い井戸型ポテンシャル中のシュレーディンガー方程式のヤコビの楕円関数による展開 with C++20

# はじめに
この記事は，[物理学アドベントカレンダー2024 Advent Calendar 2024](https://adventar.org/calendars/10569 "物理学アドベントカレンダー2024 Advent Calendar 2024") 24日目の記事です。この記事では，シュレーディンガー（Schrödinger）方程式の最も簡単な場合である、無限に深い井戸型ポテンシャル中のシュレーディンガー方程式を、変分法で数値的に解く場合を解説してみたいと思います（数値計算のコードにはC++20を用いることにします）。そしてこの場合、解析解が三角関数として求められることから、試行波動関数にはヤコビ（Jacobi）の楕円関数を用いてみようと思います。**なお、この記事の半分以上（C++20のコードも含めて）はChatGPT o1によって書かれています**。

<font color="red">**2024年12月30日
コメントでのご指摘を受けて記事を修正しました。**</font>

# シュレーディンガー方程式の変分法による解法
まず、シュレーディンガー方程式の変分法による解法を簡単に説明します。
変分法では、試行波動関数$\psi_{\mathrm{trial}}$を用いて期待値として定義されるエネルギー汎関数

```math
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}}$ についても、このエネルギー期待値は**系の基底エネルギーの値以上となる**ことが知られています。

```math
E[\psi_{\mathrm{trial}}] \ge E_{\mathrm{ground}}
```

したがって、試行波動関数を適切にパラメータ化し、そのパラメータを最適化 (variation)することで、実際の基底状態エネルギーを上から近似的に求めることが可能となります。試行波動関数は完全系であればいかなる関数でも構いませんが、物理的に適切な関数を選べばエネルギー期待値はより実際の基底エネルギーに近くなります。具体的には、例えば試行波動関数$\psi_{\mathrm{trial}}$を適当な関数$\chi(x)$の$N$個の線形結合として表し、

```math
\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)
```

とすると、エネルギー汎関数は、

```math
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}}
```

となります。ここで、

```math
H_{pq}=\left\langle \chi_{p}\right\vert H\left\vert \chi_{q}\right\rangle
```

```math
S_{pq}=\left\langle \chi_{p}\right\vert \left.  \chi_{q}\right\rangle
```

です。ここでは前者をハミルトニアン積分（あるいは行列と見てハミルトニアン行列）、後者を重なり積分（あるいは行列と見て重なり行列）と呼ぶことにします。ここで、パラメータ$C_n$の最適化を行うとすると、エネルギー汎関数の$C_n$に関する導関数がゼロになるという条件から、

```math
%
%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)
```

となります。これは結局、一般化固有値問題

```math
\mathbf{HC}=E\mathbf{SC}
```

を解くことに帰着します。

# 無限に深い井戸型ポテンシャル中のシュレーディンガー方程式のヤコビの楕円関数による展開

以下では、1次元の無限に深い井戸型ポテンシャル

```math
V(x) = 
\begin{cases}
0 & (0 \le x \le 1),\\
\infty & (x < 0 \ \text{または}\ x > 1)
\end{cases}
```

での定常シュレーディンガー方程式

```math
-\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) $ を用いた基底関数展開で扱う場合の、ハミルトニアン行列

```math
H_{ij} \;=\; \langle \chi_i \,\bigl|\bigr.\, \hat{H} \,\bigl|\bigr.\, \chi_j \rangle
```

および重なり行列

```math
S_{ij} \;=\; \langle \chi_i \,\bigl|\bigr.\, \chi_j \rangle
```

の具体的な形を示します。

## 1. 基底関数χ(r)の構築

無限に深い井戸型ポテンシャル中の波動関数$\psi(x)$では、境界条件

```math
\psi(0) = 0, 
\quad
\psi(1) = 0
```

を満たす必要があります。ヤコビの楕円関数 $ \mathrm{sn}(u,m) $ は

```math
\mathrm{sn}(0,m) = 0
```

が成り立つので、$x = 0$で消える条件を自然に満たします。さらに$x = 1$でもゼロとなるようにするには、

```math
\mathrm{sn}\bigl(\omega_n \cdot 1,\,m\bigr) = 0
```

が必要です。

ヤコビの楕円関数 $ \mathrm{sn}(u,m) $ には周期的性質があり、第一種完全楕円積分$ K(m) $を用いると

```math
\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) $)で零点を持つことが知られています。したがって

```math
\omega_n = 2\,n\,K(m)
\quad
(n = 1, 2, 3,\dots)
```

とすれば、

```math
\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,
```

となり、境界条件を満たします。

よって、たとえば基底関数を

```math
\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) $を用いて、重なり行列

```math
S_{ij}
\;=\;
\langle \chi_i \mid \chi_j \rangle
\;=\;
\int_{0}^{1} \chi_i(x)\,\chi_j(x)\,dx
```

を計算すると、

```math
S_{ij}
\;=\;
\int_{0}^{1} 
  N_i\,N_j \,
  \mathrm{sn}(\omega_i x, m)\,\mathrm{sn}(\omega_j x, m)\,
dx.
```

したがって、

```math
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 $のときの規格化条件

```math
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$となる場合、ハミルトニアンは運動エネルギーの演算子のみ

```math
\hat{H}
= 
-\frac{\hbar^2}{2M}\,\frac{d^2}{dx^2}
```

で与えられます。よって

```math
H_{ij}
= 
\int_{0}^{1}
  \chi_i(x)
  \biggl[
    -\frac{\hbar^2}{2M}\,\frac{d^2}{dx^2}
  \biggr]
  \chi_j(x)\,
dx.
```

これを部分積分により整理します。一般に

```math
\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 $であるため、境界項は消えます。すると、

```math
\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.
```

となり、したがって

```math
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) $を微分すると、ヤコビの楕円関数の性質

```math
\frac{d}{du}\,\mathrm{sn}(u,m)
=
\mathrm{cn}(u,m)\,\mathrm{dn}(u,m)
```

を用いて、

```math
\frac{d}{dx}\,\mathrm{sn}(\omega_n x, m)
=
\omega_n \,\mathrm{cn}(\omega_n x, m)\,\mathrm{dn}(\omega_n x, m).
```

よって

```math
\frac{d}{dx}\chi_n(x)
=
N_n \,\omega_n \,\mathrm{cn}(\omega_n x, m)\,\mathrm{dn}(\omega_n x, m).
```

すると、

```math
\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).
```

となります。ゆえに結局、

```math
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)$と取り、

```math
\chi_n(x) = N_n\,\mathrm{sn}(\omega_n x,m)
```

を基底関数系として用いるとき、ハミルトニアン行列及び重なり行列は次のように書けます。

- **重なり行列:**

  ```math
  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.
  ```

- **ハミルトニアン行列:**

  ```math
  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$だけ取り出して、

```math
\mathbf{HC}=E\mathbf{SC}
```

を解けば、今回の目的が達せられたことになります。

# 今回の数値計算で用いたC++20のソースコード
今回は、モジュラス$m$を$0.5$に固定し、$N$を可変として以下のC++20のソースコードで数値計算を行いました。ただし、ヤコビの楕円関数の計算とGauss-Legendre積分に[Boost C++ Libraries](https://www.boost.org/ "https://www.boost.org/")を、一般化固有値問題を解くために[Eigen](https://eigen.tuxfamily.org/index.php?title=Main_Page "https://eigen.tuxfamily.org/index.php?title=Main_Page")ライブラリを用いました。また、コンパイルにはgcc (g++) 13などのC++20に対応したコンパイラが必要な点に注意してください。なお、計算に時間がかかるのでOpenMPでスレッド並列化を行っています。

```cpp
#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 <omp.h>       // for #pragma omp parallel
#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 N           = 50;  // 基底数
static auto constexpr GAUSS_ORDER = 6 * N;
static auto constexpr HBAR        = 1.0;
static auto constexpr MASS        = 1.0;
static auto constexpr OUTPUTNUM   = 1000;
static auto constexpr FILE_OUTPUT = 0;

//============================================================
// 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)
    auto const Kval = std::comp_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)
{
    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,
                   std::int32_t gauss_order = 32)
{
    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);

#pragma omp parallel for
    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 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 const & Evals = diag.eigenvalues;
    auto const & 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++) {
        auto const E_exact = std::numbers::pi * std::numbers::pi * 0.5 * static_cast<double>((k + 1) * (k + 1));
        std::cout << std::format("E[{:>1}]        = {:>10.6f}\nE(exact)[{:>1}] = {:>10.6f}\n\n", k, Evals[k], k, E_exact);
    }

    // ---------------------------
    // 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)
    // ---------------------------

    if constexpr (FILE_OUTPUT == 1) {
        std::ofstream ofs("wavefunctions_compare_50.txt");
        if (!ofs) {
            std::cerr << "Cannot open wavefunctions_compare_50.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 (auto 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 (auto 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に書いてもらいました。以下に示します。

```julia
#!/usr/bin/env julia

using LinearAlgebra
using FastGaussQuadrature
using Elliptic
using Base.Threads  # スレッド並列用
using Printf

# ============================================
# 定数
# ============================================
const N           = 1000   # 基底数
const GAUSS_ORDER = 6 * N # Gauss–Legendreの次数
const HBAR        = 1.0
const MASS        = 1.0
const OUTPUTNUM   = 1000  # 波動関数出力用の分割数 (0～1を何等分するか)

# ============================================
# Gauss–Legendre積分
#  [x_min, x_max] 上で関数func(x)を積分
# ============================================
function numeric_integral(func::Function, x_min::Float64, x_max::Float64;
                          order::Int = GAUSS_ORDER)
    # FastGaussQuadrature.jl の gausslegendre(order) は
    # [-1, +1] 区間の節点 x[] と重み w[] を返す
    x, w = gausslegendre(order)
    # [-1, +1] を [x_min, x_max] に写像
    #   t = ((x_max - x_min)/2)*x + (x_max + x_min)/2
    #   dt/dx = (x_max - x_min)/2
    half_length = 0.5 * (x_max - x_min)
    mid_point   = 0.5 * (x_max + x_min)
    result = 0.0
    for i in eachindex(x)
        t = half_length * x[i] + mid_point
        result += w[i] * func(t)
    end
    return half_length * result
end

# ============================================
# Jacobi elliptic sn, cn, dn をまとめて返す
#
# Elliptic.jl の ellipj(u, m) は
#   sn, cn, dn, φ = ellipj(u, m)
# の形で返してくれる。
# ここではsn, cn, dn のみ使用。
# ============================================
struct JacobiElliptic
    sn::Float64
    cn::Float64
    dn::Float64
end

function jacobi_elliptic(u::Float64, m::Float64)
    sn_, cn_, dn_ = ellipj(u, m)  # Elliptic.jl
    return JacobiElliptic(sn_, cn_, dn_)
end

# ============================================
# \omega_n = 2*(n+1)*K(m)
#   (n は 0-based, n+1 が本来の"モード番号")
# ============================================
function omega_n(n::Int, m::Float64)
    # 完全楕円積分K(m): Ellipticのellipk(m) で計算
    Kval = Elliptic.K(m)
    return 2.0 * (n+1) * Kval
end

# ============================================
# phi_n_raw(x, n, m):
#   sn( ωₙ * x, m )  [正規化定数はまだ掛けない版]
# ============================================
function phi_n_raw(x::Float64, n::Int, m::Float64)
    w = omega_n(n, m)
    je = jacobi_elliptic(w*x, m)
    return je.sn  # sn(...)
end

# ============================================
# 正規化係数 Nₙ = 1 / sqrt( ∫₀¹ sn²(ωₙ x, m) dx )
# ============================================
function normalization_factor(n::Int, m::Float64)
    integrand(x) = phi_n_raw(x, 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(x) = (phi_n_raw(x, i, m) * Norms[i+1]) *
                   (phi_n_raw(x, j, m) * Norms[j+1])
    return numeric_integral(integrand, 0.0, 1.0)
end

# ============================================
# ハミルトニアン行列 Hᵢⱼ = (ℏ² / 2m) ∫₀¹ φᵢ'(x)* φⱼ'(x) dx
#
#   φₙ'(x) = Nₙ * d/dx[ sn(ωₙ x, m) ]
#           = Nₙ * [ ωₙ * cn(ωₙ x, m) * dn(ωₙ x, m) ]
# ============================================
function hamiltonian_Hij(i::Int, j::Int, m::Float64, Norms::Vector{Float64})
    # φₙ'(x) の無名関数を定義
    phi_n_prime(x, nIndex) = begin
        w  = omega_n(nIndex, m)
        je = jacobi_elliptic(w*x, m)
        return Norms[nIndex+1] * (w * je.cn * je.dn)
    end

    integrand(x) = phi_n_prime(x, i) * phi_n_prime(x, j)
    val = numeric_integral(integrand, 0.0, 1.0)
    factor = (HBAR^2) / (2.0 * MASS)
    return factor * val
end

# ============================================
# 一般化固有値問題 Hc = E S c を解いて固有値・固有ベクトルを返す
# ============================================
mutable struct DiagResult
    norms::Vector{Float64}
    eigenvalues::Vector{Float64}
    eigenvectors::Matrix{Float64}
end

function solve_eigen_problem(N::Int, m::Float64)
    # 1) 各基底の正規化定数
    norms = Vector{Float64}(undef, N)
    for n in 0:(N-1)
        norms[n+1] = normalization_factor(n, m)
    end

    # 2) 行列 S, H を構築
    S = Matrix{Float64}(I, N, N)  # とりあえず零行列で初期化
    H = Matrix{Float64}(I, N, N)

    # スレッド並列
    @threads for i in 0:(N-1)
        for j in 0:(N-1)
            S[i+1, j+1] = overlap_Sij(i, j, m, norms)
            H[i+1, j+1] = hamiltonian_Hij(i, j, m, norms)
        end
    end

    # 3) 一般化固有値問題を解く:  H v = λ S v
    # Juliaでは eigvalsや eigen が使える (LinearAlgebra.standard).
    # ただしH, Sは対称(エルミート)とみなせるので `eigen(Symmetric(H), Symmetric(S))`
    # のようにSymmetric()でラップして高速化可能 (ただし厳密にはSに正定値性が要る)
    F = eigen(Symmetric(H), Symmetric(S))

    # 固有値は F.values, 固有ベクトルは F.vectors
    eigenvals = F.values
    eigenvecs = F.vectors

    return DiagResult(norms, eigenvals, eigenvecs)
end

# ============================================
# ヤコビ楕円基底で得た k番目(0-based) の固有状態の波動関数 ψₖ(x)
#   ψₖ(x) = Σₙ cₙₖ * [ Nₙ * sn(ωₙ x, m) ]
# ============================================
function wavefunc_k(x::Float64, k::Int,
                    evecs::Matrix{Float64},
                    norms::Vector{Float64}, m::Float64, N::Int)
    s = 0.0
    for n in 0:(N-1)
        c_nk = evecs[n+1, k+1]
        # φₙ(x) = norms[n+1]*sn(ωₙ*x,m)
        s += c_nk * (phi_n_raw(x, n, m) * norms[n+1])
    end
    return s
end

# ============================================
# 無限井戸の厳密解 (節の数=node) に対応する関数
#   node = 0 -> ℓ=1
#   node = 1 -> ℓ=2
#   ...
#   ψ_exact_node(x) = √2 * sin((node+1)*π*x)
# ============================================
function wavefunc_exact(node::Int, x::Float64)
    l = node + 1
    return sqrt(2.0) * sin(l * π * x)
end

# ============================================
# メイン部分
# ============================================
function main()
    # 1) パラメータ (ここでは m=0.5 とし、基底数=250 は定数として上部で定義)
    m = 0.5
    num_states_to_compare = 5

    println("=== Jacobi-sn basis (m=$(round(m, digits=3)), N=$N) ===")

    # 2) 一般化固有値問題を解く
    diag = solve_eigen_problem(N, m)
    Evals = diag.eigenvalues
    Evecs = diag.eigenvectors

    # 3) 固有値を表示 (0～4まで)
    println("\nEigenvalues (lowest $num_states_to_compare):")
    for k in 0:(num_states_to_compare-1)
        # 無限井戸厳密解: E_exact = (π²/2) * (k+1)²
        E_exact = (π^2 * 0.5) * (k+1)^2
        println("E[$k]        = $(round(Evals[k+1], digits=6))")
        println("E(exact)[$k] = $(round(E_exact, digits=6))")
        println()
    end

    # 4) 固有値をファイルに出力して比較
    # open("energy_comparison.txt", "w") do io
    #     # ヘッダ
    #     println(io, "#N  E_numerical           E_exact")
        
    #     # 0-based で k=0..(N-1) の各固有値をファイルへ出力
    #     for k in 1:N
    #         E_num    = Evals[k]
    #         E_exact  = (π^2 * 0.5) * k^2
    #         # 小数点以下7桁、指数表記で整形 (自由に変更OK)
    #         @printf(io, "%3d %16.7f %16.7f\n", k, E_num, E_exact)
    #     end
    # end

    # 5) 波動関数をファイルに出力して比較 (コメントアウトでOFF)
    # --------------------------------------------------
    # using Printf
    # open("wavefunctions_compare.txt", "w") do io
    #     println(io, "# x   ψ_num[0..4]              ψ_exact[0..4]")
    #     for i in 0:OUTPUTNUM
    #         x = i / OUTPUTNUM
    #         # 数値解
    #         num_vals = [wavefunc_k(x, k, Evecs, diag.norms, m, N) for k in 0:(num_states_to_compare-1)]
    #         # 厳密解
    #         exact_vals = [wavefunc_exact(node, x) for node in 0:(num_states_to_compare-1)]
    #
    #         @printf(io, "%.3f", x)
    #         for val in num_vals
    #             @printf(io, " %.6f", val)
    #         end
    #         for val in exact_vals
    #             @printf(io, " %.6f", val)
    #         end
    #         println(io)
    #     end
    # end
    # println("\nWavefunctions (k=0..4) and exact solutions have been saved to wavefunctions_compare.txt")
    # println("Done.")

    return
end

# スクリプトとして実行された場合のみ main() を呼ぶ
if abspath(PROGRAM_FILE) == @__FILE__
    main()
end

```

# 結果と議論
以下に、固有値（固有エネルギー）の計算結果と厳密値の比較を示します。ここで、粒子の質量$M$は$1$とし、また$\hbar = 1$なる単位系を使いました。ここで、基底状態は波動関数の節の数が$0$に対応しており、第$n$励起状態は節の数が$n$に対応しています。

|       状態       |   N = 5   |  N = 20   |   N = 50  |   厳密    |
|:----------------:|:---------:|:---------:|:---------:|:---------:|
|     基底状態     |  4.935957 |  4.934816 |  4.934803 |  4.934802 |
|   第一励起状態   | 20.010264 | 19.740007 |  19.739282 | 19.739209 |
|   第二励起状態   | 45.023034 | 44.423137 |  44.413387 | 44.413320 |
| 　第三励起状態   | 80.041058 | 78.974932 |  78.957499 | 78.956835 |
|   第四励起状態　 | 125.064217| 123.433592| 123.371079 | 123.370055|

結果を見ると分かるとおり、$N = 5$でも十分な精度を持ち、$N = 20$では、厳密値に非常に近くなっています。また$N = 50$ではほぼ厳密値と一致しています。~~$N = 50$では極めて過小評価するという結果になっています。$N = 50$の結果は変分原理にも反しており、なぜこうなってしまうのか私には分かりませんでした。~~**この理由は、Gauss-Legendre積分の次数が足りなかったからのようです。**

次に、第一励起状態の近似波動関数（$N = 5$、$N = 20$、$N = 50$）と厳密な波動関数とを比較したプロットを以下に示します。

![qwell_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58457/ab433505-b0aa-e8f3-50fa-e0650dee72db.png)

このプロットを見ると、$N = 5$でも十分な精度を持っており、$N = 20$および$N=50$だと、プロット上は完全に厳密解に一致することが分かります。~~$N = 50$の場合は、よく分かりませんが何か特殊なことが起きているようです。~~
次に、基底状態から第999励起状態までの固有値の計算値と厳密値を比較したプロットを以下に示します（$N = 1000$としました）。

![energy_comparison.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/58457/00136f5a-08d3-9bde-6dd4-b7316a9994ab.png)

上記のプロットを見ると、変分法の示すとおり、計算された固有値はわずかに厳密な固有値より高いことが分かります。

# 全体のまとめ
本記事では、無限に深い井戸型ポテンシャル中のシュレーディンガー方程式を、ヤコビの楕円関数で展開し、三角関数で表される厳密解と比較しました。使用するヤコビの楕円関数の数（基底の数）$N$が増えると、精度がその分向上するという、変分法に従った結果となりました。~~使用するヤコビの楕円関数の数（基底の数）$N$が増えても、必ずしも精度が向上しない（それどころか変分原理に反してしまう）という興味深い結果になりました。あとChatGPT o1は本当に賢いです。~~

# 参考サイト
[1] [ChatGPT](https://chatgpt.com/ "https://chatgpt.com/")

# 参考文献
[1] J.M.ティッセン 『計算物理学』シュプリンガー・フェアラーク東京（2003）