Juliaで特殊関数を使う
量子力学にはやたらと特殊関数が登場する. 調和振動子や水素原子などの初歩的な系でさえ可視化に苦労したので, Juliaでの特殊関数の利用方法をまとめておくことにした. コメントに要望などを書き込んで頂ければ追記する.
パッケージ
Juliaの特殊関数のパッケージはSpecialFunctions.jlである. ガンマ関数, ベッセル関数, 誤差関数, 楕円積分をはじめとし, 多くの特殊関数はこちらのパッケージで利用可能である. しかし, 量子力学に登場するエルミート多項式, ラゲールの陪多項式, ルジャンドルの陪多項式などはサポートされていないので, SpecialPolynomials.jlを利用する. これらの多項式はPolynomials.jlと連携することで, 微分演算や積分演算が可能であり, 文字列として出力することも可能であるから, さらにLatexify.jlと連携することでLaTeX表示も可能である. 今回はLaTeX表示に少し拘り, LaTeXStringsも利用した. 他にもPrintf.jl, Plots.jlを使用した.
# using Pkg
# Pkg.add("Plots")
# Pkg.add("Printf")
# Pkg.add("Latexify")
# Pkg.add("Polynomials")
# Pkg.add("SpecialPolynomials")
# Pkg.add("SpecialFunctions")
using Plots
using Printf
using Latexify
using LaTeXStrings
using Polynomials
using SpecialPolynomials
using SpecialFunctions
LaTeX表示
Latexify.jlに文字列として数式を渡すことで, LaTeX出力してくれる.
latexify("ax^2+bx+c") |> display
latexify("f(x) = 1/2*x") |> display
latexify("f_1(x) = 1/2*x") |> display
残念なことに
L"f_1(x) = 1/2*x" |> display
この問題はLaTeXStrings.jlとLatexify.jlをうまく組み合わせると解決できる. env=:raw
で1/2*x
を\frac{1}{2} x
に変換し, latexstring()
に渡すことで結合できる. ついでにcdot=false
で掛け算の
latexify("1/2*x", env=:raw) |> display
latexstring("f_1(x) = ", latexify("1/2*x", env=:raw)) |> display
latexstring("f_1(x) = ", latexify("1/2*x", env=:raw, cdot=false)) |> display
\frac{1}{2} \cdot x
エルミート多項式
エルミート多項式
と表される. SpecialPolynomials.jlでサポートされており, basis(Hermite, n)(x)
のように書けばよい.
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for n in 0:11
p = basis(Hermite, n)
q = convert(Polynomial, p)
latexstring("H_{$n}(x) = ", latexify("$q", env=:raw, cdot=false)) |> display
end
例えば
# using SpecialPolynomials
p = basis(Hermite, 2)
p(5)
98.0
応用例
ラゲール多項式
注意せよ. 文献によって
に対応する式であり, basis(Laguerre{0}, n)(x)
のように書けばよい. Wikipedia(英語版)と見比べて頂きたい.
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for n in 0:7
p = basis(Laguerre{0}, n)
q = convert(Polynomial, p)
latexstring("L_{$n}(x) = ", latexify("$q", env=:raw, cdot=false)) |> display
end
を満たす解
を利用したいのであればfactorial(n)*basis(Laguerre{0}, n)(x)
とすればよい. Wikipedia(日本語版)の
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for n in 0:7
p = factorial(n)*basis(Laguerre{0}, n)
q = convert(Polynomial, p)
latexstring("\\tilde{L}_{$n}(x) = ", latexify("$q", env=:raw, cdot=false)) |> display
end
ラゲールの陪多項式
ラゲールの陪多項式
が知られている. SpecialPolynomials.jlでは一般化ラゲール多項式basis(Laguerre{k}, n-k)(x)
である.
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for n in 0:4
for k in 0:n
p = basis(Laguerre{k}, n-k)
q = convert(Polynomial, p)
latexstring("L_$n^$k(x) = ", latexify("$q", env=:raw, cdot=false)) |> display
end
end
Polynomials.jlのderivative(q,k)
を用いてラゲール多項式
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for n in 0:4
for k in 0:n
p = basis(Laguerre{0}, n)
q = convert(Polynomial, p)
r = derivative(q,k)
latexstring("L_$n^$k(x) = ", latexify("$r", env=:raw, cdot=false)) |> display
end
end
応用例
一般化ラゲール多項式
一般化ラゲール多項式
Generalized Laguerre polynomialsの和訳. 少なくともGoogle翻訳ではラゲールの陪多項式と訳されるが, 日本で使われているラゲールの陪多項式とは定義が違うようなので, 注意されたい. 調査したところ, かなり紛らわしいので, 一般化ラゲール多項式
前提
まず, 一般化ラゲール多項式
多項式 | 満たすべき方程式 | 出典 |
---|---|---|
ラゲール多項式 |
||
ラゲールの陪多項式 |
日本語のWikipedia | |
Generalized Laguerre polynomials 一般化ラゲール多項式 |
英語のWikipedia |
共通点
ラゲールの陪多項式
相違点
一般化ラゲール多項式
これは, そもそも満たすべき方程式が違うためである. 逆に言うと, ラゲールの陪多項式
対応関係
見ての通り, ラゲールの陪多項式
このことから, 水素原子の動径関数は多くが
この違いを抜きにしても規格化定数が間違っていたりするので, 数値積分などで簡単にでも確認した方がよい.
使用法
SpecialPolynomials.jlでサポートされており, basis(Laguerre{α}, n)(x)
である. もちろん
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for n in 0:4
for k in 0:n
p = basis(Laguerre{k}, n)
q = convert(Polynomial, p)
latexstring("L_$n^{($k)}(x) = ", latexify("$q", env=:raw, cdot=false)) |> display
end
end
応用例
ルジャンドル多項式
ルジャンドル多項式
で与えられる. basis(Legendre, n)(x)
のように書けばよい.
こちらと見比べて頂きたい.
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for n in 0:11
p = basis(Legendre, n)
q = convert(Polynomial, p)
latexstring("P_{$n}(x) = ", latexify("$q", env=:raw, cdot=false)) |> display
end
ルジャンドルの陪多項式
ルジャンドルの陪多項式
なお, Wikipediaの英語版だと符号がつけ加えられているが, 結局は線形微分方程式の解なので定数倍の違いは著者の好みに委ねられてしまう. 定義が違えば球面調和関数の規格化定数も変わるので注意されたい.
ここでは英語版の定義を採用する. こちらと比較して頂きたい.
# using SpecialPolynomials
# using Polynomials
# using Latexify
# using LaTeXStrings
for k in 0:4
for m in 0:k
p = basis(Legendre, k)
q = convert(Polynomial, p)
r = derivative(q, abs(m))
if m % 2 == 0
s = (-1)^m*(Polynomial([1,-2]))^(Int(m/2))*r
latexstring("P_{$k}^{$m}(x) = ", latexify("$s", env=:raw, cdot=false)) |> display
else
s = (-1)^m*(Polynomial([1,-2]))^(Int((m-1)/2))*r
latexstring("P_{$k}^{$m}(x) = \\left(", latexify("$s", env=:raw, cdot=false), "\\right)\\sqrt{1-x^2}") |> display
end
end
end
関数として利用する場合は次のようにする.
# using SpecialPolynomials
# using Polynomials
function P(k,m,x)
if m<0
return (-1)^m * factorial(l-m) / factorial(l+m) * P(k,m,x)
end
p = basis(Legendre, k)
q = convert(Polynomial, p)
r = derivative(q, abs(m))
return (-1)^m * (1-x^2)^(m/2) * r(x)
end
# Plots
for m in 0:2
plot(title=latexstring("m=$m"), xlabel=L"x", ylabel=latexstring("P_k^m(x)"), xlims=(-1,1), ylims=(-[1.2,3,15][m+1],[1.2,3,23][m+1]), legend=:bottomright)
for k in 0:[4,5,6][m+1]
plot!(x->P(k,m,x), label=latexstring("k=$k"), lw=2)
end
plot!() |> display
end
応用例
ガンマ関数
ガンマ関数factorial(x)
で利用できる. ガンマ関数はSpecialFunctions.jlでサポートされており, gamma(x)
である.
# using SpecialFunctions
println("n\tn!\tΓ(n)")
for i in 1:9
x = factorial(i-1)
y = gamma(i)
println("$i\t$x\t$y")
end
n (n-1)! Γ(n)
1 1 1.0
2 1 1.0
3 2 2.0
4 6 6.0
5 24 24.0
6 120 120.0
7 720 720.0
8 5040 5040.0
9 40320 40320.0
非整数についても
# using SpecialFunctions
gamma(1/2)^2
3.1415926535897936
描写すると次のようなグラフになる. こちらと見比べて頂きたい.
plot(zlims=(0,6), xlabel=L"y", ylabel=L"x", zlabel=L"|\Gamma(x+iy)|", camera=(45,30))
wireframe!(-5:0.2:5, -5:0.05:5, (y,x)->abs(gamma(x+im*y)))
Γ(z) = try; gamma(z); catch; NaN; end
plot(xlims=(-5,5), ylims=(-10,10), xlabel=L"x", ylabel=L"\Gamma(x)")
plot!(-5:0.01:5, x->Γ(x), label="", lw=2)
誤差関数
誤差関数はHartree-Fock法の核-電子間引力の積分などに登場する. SpecialFunctions.jlでサポートされており, 誤差関数はerf(x)
, 相補誤算関数はerfc(x)
である. たまにお世話になるので一応, 簡単な数値積分を用いて検証しておく.
# 数値積分のルーチン
function integral(func, x_min, x_max, dx)
sum = 0.0
for x in x_min:dx:x_max
sum += func(x)*dx
end
return sum
end
# using SpecialFunctions
# using Printf
@printf("%-9s\t%-9s\t%-9s\t%-9s\n", " x", " erf(x)", " 1-erfc(x)", " integral")
for x in 0:0.5:3
@printf("%9f\t%9f\t%9f\t%9f\n", x, erf(x), 1-erfc(x), 2/sqrt(pi)*integral(t->exp(-t^2), 0, x, 0.0000001))
end
x erf(x) 1-erfc(x) integral
0.000000 0.000000 0.000000 0.000000
0.500000 0.520500 0.520500 0.520500
1.000000 0.842701 0.842701 0.842701
1.500000 0.966105 0.966105 0.966105
2.000000 0.995322 0.995322 0.995322
2.500000 0.999593 0.999593 0.999593
3.000000 0.999978 0.999978 0.999978
plot(-3:0.1:3, erf, xlabel="x", ylabel="erf(x)", label="", lw=2) |> display
plot(-3:0.1:3, erfc, xlabel="x", ylabel="erfc(x)", label="", lw=2) |> display
動作環境
versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, haswell)
参考文献
- SpecialFunctions.jl
- SpecialPolynomials.jl
- Polynomials.jl, 微分演算, 積分演算
- Latexify.jl
- LaTeXStrings
- ガンマ関数
- ラゲールの陪多項式
- 水素原子のSchrödinger方程式の固有関数の動径部分
- Generalized Laguerre polynomials
- ルジャンドル多項式
- ルジャンドルの陪多項式
- 球面調和関数
- Plots.jl入門
この記事の元になったJupyter Notebookのデータは下記のリンクにある.
Discussion