JuliaLang Advent Calendar 2022 17日目です🎄
Juliaの数式処理パッケージ Symbolics.jl を使って, 量子力学の教科書によく登場するRodriguesの公式を展開していきます. 下表の直交多項式はRodriguesの公式で定義されることが多いものの, 微分が入っていて扱いづらいため, 実装には漸化式や閉形式が用いられます. Rodriguesの公式は文献によって異なることがあり, どの定義に基づく実装なのか確かめるのはなかなか大変です. そこで, Symbolics.jlに頑張ってもらいましょう.
| 名前 |
Rodriguesの公式 |
| エルミート多項式（Hermite） |
H_{n}(x) = (-1)^n \mathrm{e}^{x^2} \frac{\mathrm{d}^n}{\mathrm{d}x^n} \mathrm{e}^{-x^2} |
| ルジャンドル多項式（Legendre） |
P_n(x) = \frac{1}{2^n n!} \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( x^2-1 \right)^n |
| ルジャンドル陪多項式（associated Legendre） |
P_k^m(x) = \left( 1-x^2 \right)^{m/2} \frac{\mathrm{d}^m}{\mathrm{d}x^m} P_k(x) |
| ラゲール多項式（Laguerre） |
L_n(x) = \frac{\mathrm{e}^x}{n!} \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right) |
| ラゲール陪多項式（associated Laguerre） |
L_n^{k}(x) = \frac{\mathrm{d}^k}{\mathrm{d}x^k} L_n(x) |
| 一般化ラゲール多項式（generalized Laguerre） |
L_n^{(\alpha)}(x)= \frac{x^{-\alpha}e^x}{n!} \frac{d^n}{dx^n}\left(x^{n+\alpha}e^{-x}\right) |
※ 古い量子力学の教科書では, Laguerre多項式の定義に\frac{1}{n!}の因子を含まないL_n(x) = \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)が採用されることが多いように思われ, 当初はこの定義を採用するつもりでした. しかし, DLMFや, 一般化ラゲール多項式の定義との整合性を取るために\frac{1}{n!}の因子を加え, 上表の定義を採用することとしました.
パッケージ・環境
Symbolics.jlに加えて, 表示を整えるためにLatexify.jlとLaTeXStrings.jlを用います.
入力
using Symbolics
using Latexify
using LaTeXStrings
using SpecialFunctions
versioninfo()
出力
Julia Version 1.8.3
Commit 0434deb161 (2022-11-14 20:14 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1185G7 @ 3.00GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, tigerlake)
Threads: 1 on 8 virtual cores
Symbolics.jlの使い方
公式ドキュメントによれば, Differential()で微分演算子を定義して, expand_derivatives()で微分を展開できます. 例えば\exp(-x^2)の微分は下記のように計算できます.
@variables x
D = Differential(x)
a = exp(-x^2)
b = D(a)
c = expand_derivatives(b)
display(a)
display(b)
display(c)
2階微分は次のように計算できます.
@variables x
D = Differential(x)^2
a = exp(-x^2)
b = D(a)
c = expand_derivatives(b)
display(a)
display(b)
display(c)
上記のように計算した導関数に\exp(x^2)を掛けて整理すればHermite多項式が計算できます. simplify()のデフォルトだとうまく簡略化できないので, expand=trueのオプションを追加しましょう. すると, 下記のようにn=2のHermite多項式が求められます.
@variables x
n = 2
D = Differential(x)^n
a = (-1)^n * exp(x^2)
b = exp(-x^2)
c = a * D(b)
d = expand_derivatives(c)
e = simplify(d, expand=true)
display(a)
display(b)
display(c)
display(d)
display(e)
D = Differential(x)^0でエラーが出ますので, n=0つまり0階微分のときは恒等写像x->xになるようにするとn=0とn>0で共通のコードが使えるためすっきります. 三項演算子を使うとさらにシンプルに書けます.
@variables x
n = 0
D = n==0 ? x->x : Differential(x)^n
a = (-1)^n * exp(x^2)
b = exp(-x^2)
c = a * D(b)
d = expand_derivatives(c)
e = simplify(d, expand=true)
display(a)
display(b)
display(c)
display(d)
display(e)
なお, Int64だと階乗はfactorial(20)までしか使えないので, factorial(21)以降を使う予定の場合はfactorial(big(21))のようにBigIntを使う必要があります. 追加のパッケージ等は不要です.
julia> factorial(20)
2432902008176640000
julia> factorial(21)
OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead
Stacktrace:...
エルミート多項式
Hermite多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します. 閉形式については, いつも通りに関数を定義し, そこに@variables xで定義したSymbolics.jlの変数を渡せばSymbolics.jlのオブジェクトとして返ってきます. 展開した結果はDLMF 18.5.18等と見比べてください.
\begin{aligned}
H_{n}(x)
&= (-1)^n \mathrm{e}^{x^2} \frac{\mathrm{d}^n}{\mathrm{d}x^n} \mathrm{e}^{-x^2} \\
&= n! \sum_{m=0}^{\lfloor n/2 \rfloor} \frac{(-1)^m}{m! (n-2m)!}(2 x)^{n-2m}
\end{aligned}
for n in 0:6
@variables x
D = n==0 ? x->x : Differential(x)^n
a = (-1)^n * exp(x^2)
b = exp(-x^2)
c = a * D(b)
d = expand_derivatives(c)
e = simplify(d, expand=true)
H(x; n=0) = factorial(n) * sum(m -> (-1)^m // (factorial(m) * factorial(n-2*m)) * (2*x)^(n-2*m), 0:Int(floor(n/2)))
display(latexstring(
"""
\\begin{aligned}
H_{$n}(x)
= $(latexify(c, env=:raw))
&= $(latexify(e, env=:raw)) \\\\
&= $(latexify(H(x, n=n), env=:raw))
\\end{aligned}
"""))
end
ちゃんと一致していますね！ 多項式の順番が少し変ですが, 値は合っているので, 閉形式がRodriguesの公式による定義と一致することの確認としては問題ありません.
ルジャンドル多項式
Legendre多項式のRodriguesの公式を展開し, 閉形式と比較します. 閉形式については, Legendre陪多項式のm=0の場合を用います. 展開した結果はDLMF 18.5.16等と見比べてください.
\begin{aligned}
P_n(x)
&= \frac{1}{2^n n!} \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left(x^2-1 \right)^n \\
&= P_n^0(x)
\end{aligned}
for n in 0:6
@variables x
D = n==0 ? x->x : Differential(x)^n
a = 1 // (2^n * factorial(n))
b = (x^2 - 1)^n
c = a * D(b)
d = expand_derivatives(c)
e = simplify(d, expand=true)
P(x; n=0, m=0) = (1//2)^n * (1-x^2)^(m//2) * sum(j -> (-1)^j * factorial(2*n-2*j) // (factorial(j) * factorial(n-j) * factorial(n-2*j-m)) * x^(n-2*j-m), 0:Int(floor((n-m)/2)))
display(latexstring(
"""
\\begin{aligned}
P_{$(n)}(x)
= $(latexify(c, env=:raw))
&= $(latexify(e, env=:raw)) \\\\
&= $(latexify(P(x, n=n), env=:raw))
\\end{aligned}
"""))
end
ルジャンドル陪多項式
Legendre陪多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します. 展開した結果はMcQuarrie&Simon(1997)等と見比べてください.
\begin{aligned}
P_n^m(x)
&= \left( 1-x^2 \right)^{m/2} \frac{\mathrm{d}^m}{\mathrm{d}x^m} P_n(x) \\
&= \frac{1}{2^n} (1-x^2)^{m/2} \sum_{j=0}^{\left\lfloor\frac{n-m}{2}\right\rfloor} (-1)^j \frac{(2n-2j)!}{j! (n-j)! (n-2j-m)!} x^{(n-2j-m)}
\end{aligned}
for n in 0:4
for m in 0:n
@variables x
Dn = n==0 ? x->x : Differential(x)^n
Dm = m==0 ? x->x : Differential(x)^m
a = 1 // (2^n * factorial(n))
b = (x^2 - 1)^n
c = (1 - x^2)^(m//2) * Dm(a * Dn(b))
d = expand_derivatives(c)
e = simplify(d, expand=true)
P(x; n=0, m=0) = (1//2)^n * (1-x^2)^(m//2) * sum(j -> (-1)^j * factorial(2*n-2*j) // (factorial(j) * factorial(n-j) * factorial(n-2*j-m)) * x^(n-2*j-m), 0:Int(floor((n-m)/2)))
display(latexstring(
"""
\\begin{aligned}
P_{$n}^{$m}(x)
= $(latexify(c, env=:raw))
&= $(latexify(e, env=:raw)) \\\\
&= $(latexify(simplify(P(x, n=n, m=m), expand=true), env=:raw))
\\end{aligned}
"""))
end
end
ラゲール多項式
Laguerre多項式のRodriguesの公式を展開し, 閉形式と比較します. 閉形式についてはLaguerre陪多項式のk=0および一般化Laguerre多項式の\alpha=0の場合を用います. 展開した結果はDLMF等と見比べてください.
\begin{aligned}
L_n(x)
&= \frac{\mathrm{e}^x}{n!} \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right) \\
&= L_n^{0}(x) \\
&= L_n^{(0)}(x)
\end{aligned}
for n in 0:6
@variables x
D = n==0 ? x->x : Differential(x)^n
a = exp(x) / factorial(n)
b = exp(-x) * x^n
c = a * D(b)
d = expand_derivatives(c)
e = simplify(d, expand=true)
AL(x; n=0, k=0) = sum(m -> (-1)^(m+k) * factorial(n) // (factorial(m) * factorial(m+k) * factorial(n-m-k)) * x^m, 0:n-k)
GL(x; n=0, α=0) = sum(k -> (-1)^(k) * binomial(n+α,n-k) * x^k // factorial(k), 0:n)
display(latexstring(
"""
\\begin{aligned}
L_{$n}(x)
= $(latexify(c, env=:raw))
&= $(latexify(e, env=:raw)) \\\\
&= $(latexify(AL(x, n=n, k=0), env=:raw)) \\\\
&= $(latexify(GL(x, n=n, α=0), env=:raw))
\\end{aligned}
"""))
end
ラゲール陪多項式
Laguerre陪多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します. 閉形式については, いつも通りに関数を定義し, そこに@variables xで定義したSymbolics.jlの変数を渡せばSymbolics.jlのオブジェクトとして返ってきます. 展開した結果はWikipedia等と見比べてください.
\begin{aligned}
L_n^{k}(x)
&= \frac{\mathrm{d}^k}{\mathrm{d}x^k} L_n(x) \\
&= \sum_{m=0}^{n-k} (-1)^{m+k} \frac{n!}{m!(m+k)!(n-m-k)!} x^m \\
&= (-1)^k L_{n-k}^{(k)}(x)
\end{aligned}
for n in 0:4
for k in 0:n
@variables x
Dn = n==0 ? x->x : Differential(x)^n
Dk = k==0 ? x->x : Differential(x)^k
a = exp(x) / factorial(n)
b = exp(-x) * x^n
c = Dk(a * Dn(b))
d = expand_derivatives(c)
e = simplify(d, expand=true)
AL(x; n=0, k=0) = sum(m -> (-1)^(m+k) * factorial(n) // (factorial(m) * factorial(m+k) * factorial(n-m-k)) * x^m, 0:n-k)
GL(x; n=0, α=0) = sum(k -> (-1)^(k) * binomial(n+α,n-k) * x^k // factorial(k), 0:n)
display(latexstring(
"""
\\begin{aligned}
L_{$n}^{$k}(x)
= $(latexify(c, env=:raw))
&= $(latexify(e, env=:raw)) \\\\
&= $(latexify(AL(x, n=n, k=k), env=:raw)) \\\\
&= $(latexify((-1)^k * GL(x, n=n-k, α=k), env=:raw))
\\end{aligned}
"""))
end
end
一般化ラゲール多項式
一般化Laguerre多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します.
\begin{aligned}
L_n^{(\alpha)}(x)
&= \frac{x^{-\alpha}e^x}{n!} \frac{d^n}{dx^n}\left(x^{n+\alpha}e^{-x}\right) \\
&= \sum_{k=0}^n(-1)^k \left(\begin{array}{l} n+\alpha \\ n-k \end{array}\right) \frac{x^k}{k !} \\
&= \sum_{k=0}^n(-1)^k \frac{\Gamma(\alpha+n+1)}{\Gamma(\alpha+k+1)\Gamma(n-k+1)} \frac{x^k}{k !}
\end{aligned}
二項係数\left(\begin{array}{l} a \\ b \end{array}\right)はbinomial(a,b)で計算できます. ただし, \alphaが整数の場合以外はbinomial(a,b)は計算できないので, \left(\begin{array}{l} a \\ b \end{array}\right) = \frac{a !}{b !(a-b) !} = \frac{\Gamma(a+1)}{\Gamma(b+1)\Gamma(a-b+1)} を利用し, 階乗表示を経由してガンマ関数から計算することにします.（そもそも, \alphaが整数なら陪関数との対応関係を使えば十分です. Morse振動子の波動関数では\alphaは非整数なので, 非整数に対応しておく必要があります.） なお, Laguerre多項式の定義によっては, 一般化Laguerre多項式はLaguerre多項式の一般化にならないことに注意してください. \alpha=0の場合はL_n^{(0)}(x) = \frac{e^x}{n!} \frac{d^n}{dx^n}\left(x^{n}e^{-x}\right)となり, Laguerre多項式の定義をL_n(x) = e^x \frac{d^n}{dx^n}\left(x^{n}e^{-x}\right)とした場合とは一致しません. このノートではLaguerre多項式をL_n(x) = \frac{e^x}{n!} \frac{d^n}{dx^n}\left(x^{n}e^{-x}\right)と定義したので, 一般化になっています.
for n in 0:4
for α in 0:n
@variables x
D = n==0 ? x->x : Differential(x)^n
a = exp(x) * x^(-α) / factorial(n)
b = exp(-x) * x^(n+α)
c = a * D(b)
d = expand_derivatives(c)
e = simplify(d, expand=true)
GL1(x; n=0, α=0) = sum(k -> (-1)^(k) * binomial(n+α,n-k) * x^k // factorial(k), 0:n)
GL2(x; n=0, α=0) = sum(k -> (-1)^(k) * (Int(gamma(α+n+1)) // Int((gamma(α+1+k)*gamma(n-k+1)))) * x^k // factorial(k), 0:n)
display(latexstring(
"""
\\begin{aligned}
L_{$n}^{($α)}(x)
= $(latexify(c, env=:raw))
&= $(latexify(e, env=:raw)) \\\\
&= $(latexify(GL1(x, n=n, α=α), env=:raw)) \\\\
&= $(latexify(GL2(x, n=n, α=α), env=:raw))
\\end{aligned}
"""))
end
end
まとめ
おかげで各種の閉形式についてRodriguesの公式との整合性が簡単に確認できるようになりました. Symbolics.jlはまだ発展途中のパッケージですが, かなり活躍してくれましたね. 自分でルールを追加して機能を拡張できるようなので, これから比較にならないほど発展すると思われます. 現状のシステムに強いて注文を付けるとすれば, 多項式の各項を次数の順番でソートできると良いですね. Issueに上がっているので, 対応されるものと思われます.
参考文献
https://symbolics.juliasymbolics.org/dev/manual/derivatives/
| 文献 |
|
|
|
|
|
|
| cpprefjp |
hermite |
legendre |
assoc_legendre |
laguerre |
assoc_laguerre |
|
| The Digital Library of Mathematical Functions (DLMF) |
18.3 Table1, 18.5 Table1, 18.5.13, 18.5.18
|
18.3 Table1, 18.5 Table1, 18.5.16
|
|
18.3 Table1, 18.5 Table1, 18.5.17
|
|
18.3 Table1, 18.5 Table1, 18.5.12
|
| L. D. Landau, E. M. Lifshitz, Quantum Mechanics (Pergamon Press, 1965) |
p.595 (a.4), (a.6) |
p.598 (c.1) |
p.598 (c.4) |
p.603 (d.13) |
p.603 (d.13) |
|
| L. I. Schiff, Quantum Mechanics (McGraw-Hill Book Company, 1968) |
p.71 (13.12) |
|
p.79 (14.12) |
|
p.93 (16.19) |
|
| A. Messiah, Quanfum Mechanics (Dover Publications, 1999) |
p.491 (B.59) |
p.493 (B.72),p.494 Table
|
p.493 (B.72) |
p.483 (B.12) |
p.483 (B.12) |
|
| W. Greiner, Quantum Mechanics: An Introduction Third Edition (Springer, 1994) |
p.152 (7.22) |
p.83 (4) |
p.83 (5) |
p.149 (21) |
|
|
| D. J. Griffiths, Introduction to Quantum Mechanics (Prentice Hall, 1995) |
p.41 Table 2.1, p.43 (2.70)
|
p.126 (4.28), p.96 Table3.1
|
p.126 (4.27) |
p.139 (4.88), p.140 Table4.4
|
p.139 (4.87), p.140 Table4.5
|
|
| D. A. McQuarrie, J. D. Simon, Physical Chemistry: A Molecular Approach (University Science Books, 1997) |
p.170 Table 5.2 |
p.195 Table6.1 |
p.196 (6.26), p.196 Table6.2
|
|
p.207 Table6.4 |
|
| P. W. Atkins, J. De Paula, Atkins' Physical Chemistry, 8th edition (W. H. Freeman, 2008) |
p.293 Table 9.1 |
|
|
|
|
|
| J. J. Sakurai, J. Napolitano, Modern Quantum Mechanics Third Edition (Cambridge University Press, 2021) |
p.524 (B.29) |
|
|
|
p.245 Problem 3.30.b |
|
| J. P. Dahl, M. Springborg, J. Chem. Phys. 88, 4535 (1988) |
|
|
|
|
|
(62), (63) |
| W. K. Shao, Y. He, J. Pan, J. Nonlinear Sci. Appl., 9, 5, 3388 (2016) |
|
|
|
|
|
(1.6) |
https://gist.github.com/ohno/287799f1f591fb3e6541c7ecee69cd03
Discussion