💡
FiniteDifferenceMatrices.jlの紹介

2024/12/05に公開
Julia
Julia Advent Calendar 2024 の5日目です🎄

 はじめにFiniteDifferenceMatrices.jl というパッケージを開発・公開しました🎉 このパッケージの宣伝も兼ねて, 調和振動子における変分法, 固有値問題, 虚時間発展の計算例について解説します.

 FiniteDifferenceMatrices.jl以下では, FiniteDifferenceMatrices.jlの定式化と使い方を解説していきます. まずインストールしましょう. 下記を実行してください.
インストール時, 初回のみ
import Pkg; Pkg.add("FiniteDifferenceMatrices.jl")
使用時は下記を実行します.
使用時, 毎回
using FiniteDifferenceMatrices
2階微分の2次精度の中心差分近似は

\frac{\mathrm{d}^{2}f}{\mathrm{d} x^{2}}(x) = \frac{f(x-\Delta x) - 2f(x) + f(x+\Delta x)}{\Delta x^{2}} + O(\Delta x^{2}),
です. これはTaylor展開を使ってゴリゴリ計算していけば求められます. 具体的には, 中心の異なる2つのTaylor級数

\begin{aligned}
f(x+\Delta x)
&= f(x)
+ \frac{\mathrm{d} f}{\mathrm{d} x}(x) \Delta x
+ \frac{1}{2!} \frac{\mathrm{d}^{2} f(x)}{\mathrm{d} x^{2}} \Delta x^{2}
+ \frac{1}{3!} \frac{\mathrm{d}^{3} f(x)}{\mathrm{d} x^{3}} \Delta x^{3}
+ O\left(\Delta x^{4}\right),
\\
f(x-\Delta x)
&= f(x)
- \frac{\mathrm{d} f}{\mathrm{d} x}(x) \Delta x
+ \frac{1}{2!} \frac{\mathrm{d}^{2} f(x)}{\mathrm{d} x^{2}} \Delta x^{2}
- \frac{1}{3!} \frac{\mathrm{d}^{3} f(x)}{\mathrm{d} x^{3}} \Delta x^{3}
+ O\left(\Delta x^{4}\right),
\end{aligned}
の和を取って整理すると,

\begin{aligned}
  f(x-\Delta x)
+ f(x+\Delta x)
&=
2f(x)
+ \frac{\mathrm{d}^{2} f}{\mathrm{d} x^{2}}(x) \Delta x^{2}
+ O\left(\Delta x^{4}\right)
\\
\frac{\mathrm{d}^{2} f}{\mathrm{d} x^{2}}(x) \Delta x^{2}
&=
  f(x-\Delta x)
- 2f(x)
+ f(x+\Delta x)
- O\left(\Delta x^{4}\right)
\\
\frac{\mathrm{d}^{2} f}{\mathrm{d} x^{2}}(x)
&=
\frac{f(x-\Delta x) - 2f(x) + f(x+\Delta x)}{\Delta x^{2}}
- \frac{O\left(\Delta x^{4}\right)}{\Delta x^{2}}
\\
\frac{\mathrm{d}^{2} f}{\mathrm{d} x^{2}}(x)
&=
\frac{f(x-\Delta x) - 2f(x) + f(x+\Delta x)}{\Delta x^{2}}
+ O\left(\Delta x^{2}\right).
\end{aligned}
が得られます. n階微分のm次精度の差分近似は

\frac{\mathrm{d}^n f}{\mathrm{d}x^n}(x) = \frac{1}{\Delta x^n} \sum_{i} c_i f(x+i\Delta x) + O(\Delta x^m),
と書けます. 係数c_iは先ほどの例と同じように求められますが, nやmを変えるたびに計算し直すのはかなり大変です. 大きいn,mについては文献値もありません. そこで, Discourseの投稿を元に, この係数を出力する関数 fdcoefficient を実装しました. 2階微分の4次精度（n=2, m=2）のの中心差分 :c の例を示します.
n=2,m=2,中心差分の例
julia> fdcoefficient(n=2, m=2, d=:c)
Dict{Int64, Rational{Int64}} with 3 entries:
  0  => -2
  -1 => 1
  1  => 1
続いて, 2階微分の4次精度の中心差分の例は次のようになります.
n=2,m=4,中心差分の例
julia> fdcoefficient(n=2, m=4, d=:c)
Dict{Int64, Rational{Int64}} with 5 entries:
  0  => -5//2
  -1 => 4//3
  2  => -1//12
  -2 => -1//12
  1  => 4//3
一応, Wikipediaと一致する結果が得られました. 実装が本当に合っているか不安だったので, 論文やWikipediaに載っている表を愚直に手作業でハードコードし, 出力が一致することを確認するテストを実装してあります.
https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/API/#FiniteDifferenceMatrices.fdcoefficient-Tuple{}
次に, 差分近似を行列で表してきましょう.　上記の近似を行列で表すと次のように書けます.

\left(\begin{array}{c}
  \psi''(x_1) \\
  \psi''(x_2) \\
  \psi''(x_3) \\
  \vdots \\
\end{array}\right)
\simeq
\frac{1}{\Delta x^2}
\left(\begin{array}{ccccccc}
  -2 &  1 &  0 & \ldots \\
   1 & -2 &  1 & \ldots \\
   0 &  1 & -2 & \ldots \\
  \vdots & \vdots & \vdots & \ddots & \vdots \\
\end{array}\right)
\left(\begin{array}{c}
  \psi(x_1) \\
  \psi(x_2) \\
  \psi(x_3) \\
  \vdots \\
\end{array}\right)
行列の記法に慣れない人のために抜き出すと, \psi''(x_2)は次のようになっています.

\psi''(x_2) = \frac{\psi(x_1) - 2f(x_2) + \psi(x_3)}{\Delta x^{2}}
このとき両端では \psi(x_0)=0, \psi(x_{n_\mathrm{max}+1})=0 としていることに相当します. この行列部分を生成する関数 fdmatrix を実装しました. n階, m次精度, 中心差分(:c), \Delta x=1 (h=1//1) のときの差分行列は次のように計算できます.
julia> fdmatrix(5, n=2, m=2, d=:c, h=1//1)
5×5 SparseArrays.SparseMatrixCSC{Rational{Int64}, Int64} with 13 stored entries:
 -2   1   ⋅   ⋅   ⋅
  1  -2   1   ⋅   ⋅
  ⋅   1  -2   1   ⋅
  ⋅   ⋅   1  -2   1
  ⋅   ⋅   ⋅   1  -2
https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/API/#FiniteDifferenceMatrices.fdmatrix-Tuple{Int64}
これと同様に, ハミルトニアンを疎行列として近似できます. 2次精度の中心差分近似における, ハミルトニアン \hat{H} = - \frac{\hbar^2}{2m} \frac{\mathrm{d}^2}{\mathrm{d}x^2} + V(x) , 波動関数 \psi(x), エネルギー E はそれぞれ下記のように書けます.

\pmb{H}
= - \frac{\hbar^2}{2m} \cdot \frac{1}{\Delta x^2}
  \left(\begin{array}{ccccccc}
    -2 & 1 & 0 & \ldots \\
    1 & -2 & 1 & \ldots \\
    0 & 1 & -2 & \ldots \\
    \vdots & \vdots & \vdots & \ddots \\
  \end{array}\right)
  +
  \left(\begin{array}{ccccccc}
    V(x_1) & 0 & 0 & \ldots \\
    0 & V(x_2) & 0 & \ldots \\
    0 & 0 & V(x_3) & \ldots \\
    \vdots & \vdots & \vdots & \ddots \\
  \end{array}\right),

\vec{\psi}
=
\left(\begin{array}{c}
  \psi(x_1) \\
  \psi(x_2) \\
  \psi(x_3) \\
  \vdots \\
\end{array}\right),

E = \frac{\langle\psi|\hat{H}|\psi\rangle}{\langle\psi|\psi\rangle}
\simeq \frac{\vec{\psi}^\dagger \pmb{H} \vec{\psi}}{\vec{\psi}^\dagger \vec{\psi}}.
実際に調和振動子のエネルギーが正しく計算できるか試していきましょう. 調和振動子の無次元化されたハミルトニアン, 基底状態の波動関数, 基底状態のエネルギーはそれぞれ下記のように書けます.

\hat{H} = - \frac{1}{2} \frac{\mathrm{d}^2}{\mathrm{d}x^2} + \frac{1}{2}x^2,

\psi_0(x) = \pi^{-1/4} \exp \left( -\frac{1}{2}x^2 \right),

E_0 = \frac{\langle \psi_0 | \hat{H} | \psi_0 \rangle}{\langle \psi_0 | \psi_0 \rangle} = \frac{1}{2}.
コードはほとんど数式そのままです.
入力
using FiniteDifferenceMatrices
using SparseArrays

Δx = 0.1
X = -10:Δx:10
H = -1/2 * fdmatrix(length(X), n=2, m=2, d=:c, h=Δx) + spdiagm([1/2*x^2 for x in X])
ψ = [π^(-1//4)*exp(-1/2*x^2) for x in X]
E = (ψ' * H * ψ) / (ψ' * ψ)
出力
0.49968776025398826
実際に厳密解0.5に近い値が計算できました.
https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/examples/#Discrete-Approximation-of-Hamiltonian

 変分法エネルギーが計算できるようになりましたので, 変分法も簡単に実装できます. パラメータ \alpha を含んだ試行波動関数 \psi_\mathrm{trial} を次のように定義します.

\psi_\mathrm{trial}(x) = \left( \frac{2\alpha}{\pi} \right)^{\frac{1}{4}} \exp \left( - \alpha x^2 \right).
パラメータ \alpha を変えてこの試行波動関数から計算できるエネルギー E を最小化します. 前述の解析解から, \alpha=\frac{1}{2}でエネルギーが最小になるはずです.
入力
using FiniteDifferenceMatrices
using SparseArrays
using Printf

Δx = 0.1
X = -10:Δx:10
H = -1/2 * fdmatrix(length(X), n=2, m=8, d=:c, h=Δx) + spdiagm([1/2*x^2 for x in X])

println("  α              E")
println("------------------")
for α in 0.1:0.1:1.0
  ψ = [((2*α)/π)^(1//4)*exp(-α*x^2) for x in X]
  E = (ψ' * H * ψ) / (ψ' * ψ)
  @printf("%.1f    %.9f\n", α, E)
end
出力
  α              E
------------------
0.1    1.299999995
0.2    0.725000000
0.3    0.566666667
0.4    0.512500000
0.5    0.500000000
0.6    0.508333333
0.7    0.528571428
0.8    0.556250000
0.9    0.588888888
1.0    0.624999999
実際に α = 0.5 で最小値を取っていることが確かめられました. CairoMakie.jlで波動関数もプロットしてみましょう.
入力
using CairoMakie
CairoMakie.activate!(type = "svg")

f = Figure(size=(420,300), fontsize=11.5)
ax = Axis(f[1,1], xlabel=L"x", ylabel=L"\psi(x)", limits=(-5,5,0.0,0.8), titlesize=16.5, ylabelsize=16.5, xlabelsize=16.5)
for α in 0.3:0.1:0.6
  ψ = [((2*α)/π)^(1//4)*exp(-α*x^2) for x in X]
  lines!(ax, X, ψ, label="α = $α")
end
lines!(ax, X, x -> π^(-1//4)*exp(-1/2*x^2), label="exact", color=:black, linestyle=:dash, linewidth=2)
axislegend(ax, position=:rt, rowgap=0, padding=(0,0,0,0), framevisible=false)
f
試行波動関数が α = 0.5 で厳密解に一致することも確かめられました.
https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/examples/#Variational-Method

 固有値問題上記で計算したハミルトニアン行列 \pmb{H} の固有値問題を解くことで, エネルギーの近似値と, 波動関数の近似値の配列が得られます. LinearAlgebra.jl または ArnoldiMethod.jl のいずれかを用いて, Hを渡すだけです. ベンチマークには Antique.jl を用います.
入力
using FiniteDifferenceMatrices
using SparseArrays
using LinearAlgebra
# using ArnoldiMethod
using Printf
import Antique
HO = Antique.HarmonicOscillator(k=1.0, m=1.0, ℏ=1.0)

Δx = 0.1
X = -5:Δx:5
I = keys(X)
T = -HO.ℏ^2/2/HO.m * fdmatrix(length(X), n=2, m=8, d=:c, h=Δx)
V = Diagonal([Antique.V(HO,x) for x in X])
H = T + V

# decomp, history = partialschur(H, nev=10, tol=1e-9, which=SR())
# E = real.(decomp.eigenvalues)
E, ψ = eigen(Matrix{Float64}(H))

println(" n         numerical        analytical")
println("--    --------------    --------------")
for n in 0:9
  @printf("%2d    %.12f    %.12f\n", n, E[n+1], Antique.E(HO, n=n))
end
出力
 n         numerical        analytical
--    --------------    --------------
 0    0.499999999985    0.500000000000
 1    1.500000001062    1.500000000000
 2    2.500000034551    2.500000000000
 3    3.500000555298    3.500000000000
 4    4.500006084209    4.500000000000
 5    5.500049973732    5.500000000000
 6    6.500321040771    6.500000000000
 7    7.501649561884    7.500000000000
 8    8.506863078481    8.500000000000
 9    9.523313467073    9.500000000000
このように固有値としてエネルギーの近似値が得られました. 基底状態についてはかなり良い近似になっています. 固有ベクトルψ[:,1]が基底状態（n=0）の波動関数に対応しており, ψ[:,n+1]が第n励起状態の波動関数に対応しています. CairoMakie.jlでプロットしてみましょう.
入力
using CairoMakie
CairoMakie.activate!(type = "svg")

f = Figure(size=(840,600), fontsize=11.5)
for i in 1:4
    n = i - 1
    ax = Axis(
      f[Int(ceil(i/2)),rem(i-1,2)+1],
      xlabel = L"x",
      ylabel = L"|\psi(x)|^2",
      titlesize = 16.5,
      ylabelsize = 16.5,
      xlabelsize = 16.5,
    )
  lines!(ax, -5..5, x -> Antique.ψ(HO,x,n=n)^2, label="analytical", color=:black, linestyle=:solid, linewidth=1)
  scatter!(ax, X, (ψ[:,i] / sqrt(Δx)) .^ 2, label="numerical", markersize=6)
  axislegend(ax, "n = $n", position=:rt, rowgap=0, framevisible=false, titlesize=15.0)
end
f
save("eigen.svg", f) # hide
nothing # hide
いずれも精度よく決定できていますが, 高い準位では振動の節が増えていくため, 格子間隔も小さくしないと記述が難しいであろうことが想像できます. 実際, エネルギーの誤差はn=9でn=0よりもかなり大きくなっています.
https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/examples/#Eigenvalue-Problem

 虚時間発展虚時間発展が基底状態への射影に相当することは偏微分方程式の変数分離による解法を知っていれば容易に理解できます. 興味のある方は Kosztin et al. (1996) 等を参照してください. ここでは詳細に踏み込まず, 虚時間のSchrödinger方程式

-\frac{\partial \varPsi(x,\tau)}{\partial \tau}
= \hat{H} \varPsi(x,\tau),
の形式解

\varPsi(x,\tau)
= \exp \left[- \hat{H} \tau \right] \varPsi(x,0),
を直接的に計算する方法を紹介します. Juliaには行列指数関数が最初から組み込まれていますので行列表示したハミルトニアンHを初期値のベクトルΨ₀にかけてΨ = exp(-Matrix{Float64}(H)*τ) * Ψ₀を計算すれば虚時間τに対する近似解Ψが計算できます. 行列指数関数の部分を疎行列として効率的に計算したければFastExpm.jlを使用してください.
入力
using FiniteDifferenceMatrices
using SparseArrays
using CairoMakie
CairoMakie.activate!(type = "svg")

# domain X, Hamiltonian H, initial state ψ₀
Δx = 0.1
X = -10:Δx:10
H = -1/2 * fdmatrix(length(X), n=2, m=4, d=:c, h=Δx) + spdiagm([1/2*x^2 for x in X])
Ψ₀ = [-2.5 ≤ x ≤ 2.5 ? 0.2 : 0.0 for x in X]

# plot
f = Figure(size=(420,300), fontsize=11.5)
ax = Axis(f[1,1], xlabel="x", ylabel="Ψ(x,τ)", limits=(-5,5,0.0,0.8), titlesize=16.5, ylabelsize=16.5, xlabelsize=16.5)
for τ in [0.0, 0.1, 1.0, 10.0]
  Ψ = exp(-Matrix{Float64}(H)*τ) * Ψ₀
  lines!(ax, X, Vector{Float64}(real(Ψ / sqrt(Δx * Ψ' * Ψ))), label="τ = $τ")
end
lines!(ax, X, x -> π^(-1//4)*exp(-1/2*x^2), label="g.s.", color=:black, linestyle=:dash, linewidth=2)
axislegend(ax, position=:rt, rowgap=0, padding=(0,0,0,0), framevisible=false)
f
初期値は水色の台形ですが, 右辺がうまく表示されずに切れてしまっています. 計算上は問題ありませんので気にしないでください. このように, 虚時間発展を計算すると波動関数の励起状態の成分を少なくして基底状態の成分だけを取り出す（基底状態に射影する）ことができます.
https://ohno.github.io/FiniteDifferenceMatrices.jl/stable/examples/#Imaginary-Time-Evolution

 まとめFiniteDifferenceMatrices.jl というパッケージを紹介しました. 差分近似の係数については先行研究と一致することがテストされています. 変分法, 固有値問題, 虚時間発展などについて教育的な計算例を提供しました.

 宣伝下記の要綱で Julia in Physics 2024 を開催します！ 奮ってご参加ください🎉
【日時】2024年12月14日(土) 13:00 - 19:00

【方式】対面 / Zoom ハイブリット

【会場】東京大学 本郷地区キャンパス

【最寄】根津駅（千代田線）

【詳細】https://ohno.github.io/julia_in_physics_2024/
https://ohno.github.io/julia_in_physics_2024/
はじめに

FiniteDifferenceMatrices.jl

変分法

固有値問題

虚時間発展

まとめ

宣伝

Discussion
