JuliaLang Advent Calendar 2022 17日目です🎄
Juliaの数式処理パッケージ Symbolics.jl を使って, 量子力学の教科書によく登場するRodriguesの公式を展開していきます. 下表の直交多項式はRodriguesの公式で定義されることが多いものの, 微分が入っていて扱いづらいため, 実装には漸化式や閉形式が用いられます. Rodriguesの公式は文献によって異なることがあり, どの定義に基づく実装なのか確かめるのはなかなか大変です. そこで, Symbolics.jlに頑張ってもらいましょう.
名前 |
Rodriguesの公式 |
エルミート多項式(Hermite) |
Hn(x)=(−1)nex2dxndne−x2 |
ルジャンドル多項式(Legendre) |
Pn(x)=2nn!1dxndn(x2−1)n |
ルジャンドル陪多項式(associated Legendre) |
Pkm(x)=(1−x2)m/2dxmdmPk(x) |
ラゲール多項式(Laguerre) |
Ln(x)=n!exdxndn(e−xxn) |
ラゲール陪多項式(associated Laguerre) |
Lnk(x)=dxkdkLn(x) |
一般化ラゲール多項式(generalized Laguerre) |
Ln(α)(x)=n!x−αexdxndn(xn+αe−x) |
※ 古い量子力学の教科書では, Laguerre多項式の定義にn!1の因子を含まないLn(x)=exdxndn(e−xxn)が採用されることが多いように思われ, 当初はこの定義を採用するつもりでした. しかし, DLMFや, 一般化ラゲール多項式の定義との整合性を取るためにn!1の因子を加え, 上表の定義を採用することとしました.
パッケージ・環境
Symbolics.jlに加えて, 表示を整えるためにLatexify.jlとLaTeXStrings.jlを用います.
Symbolics.jlの使い方
公式ドキュメントによれば, Differential()
で微分演算子を定義して, expand_derivatives()
で微分を展開できます. 例えばexp(−x2)の微分は下記のように計算できます.
dxde−x2
−2xe−x2
2階微分は次のように計算できます.
dx2d2e−x2
−2e−x2+4x2e−x2
上記のように計算した導関数にexp(x2)を掛けて整理すればHermite多項式が計算できます. simplify()
のデフォルトだとうまく簡略化できないので, expand=true
のオプションを追加しましょう. すると, 下記のようにn=2のHermite多項式が求められます.
ex2dx2d2e−x2
(−2e−x2+4x2e−x2)ex2
−2+4x2
D = Differential(x)^0
でエラーが出ますので, n=0つまり0階微分のときは恒等写像x->x
になるようにするとn=0とn>0で共通のコードが使えるためすっきります. 三項演算子を使うとさらにシンプルに書けます.
ex2e−x2
ex2e−x2
なお, Int64だと階乗はfactorial(20)
までしか使えないので, factorial(21)
以降を使う予定の場合はfactorial(big(21))
のようにBigIntを使う必要があります. 追加のパッケージ等は不要です.
エルミート多項式
Hermite多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します. 閉形式については, いつも通りに関数を定義し, そこに@variables x
で定義したSymbolics.jlの変数を渡せばSymbolics.jlのオブジェクトとして返ってきます. 展開した結果はDLMF 18.5.18等と見比べてください.
Hn(x)=(−1)nex2dxndne−x2=n!m=0∑⌊n/2⌋m!(n−2m)!(−1)m(2x)n−2m
H0(x)=ex2e−x2=1=1
H1(x)=−ex2dxde−x2=2x=2x
H2(x)=ex2dx2d2e−x2=−2+4x2=−2+4x2
H3(x)=−ex2dx3d3e−x2=8x3−12x=8x3−12x
H4(x)=ex2dx4d4e−x2=12−48x2+16x4=12−48x2+16x4
H5(x)=−ex2dx5d5e−x2=32x5+120x−160x3=32x5+120x−160x3
H6(x)=ex2dx6d6e−x2=−120+720x2+64x6−480x4=−120+720x2+64x6−480x4
ちゃんと一致していますね! 多項式の順番が少し変ですが, 値は合っているので, 閉形式がRodriguesの公式による定義と一致することの確認としては問題ありません.
ルジャンドル多項式
Legendre多項式のRodriguesの公式を展開し, 閉形式と比較します. 閉形式については, Legendre陪多項式のm=0の場合を用います. 展開した結果はDLMF 18.5.16等と見比べてください.
Pn(x)=2nn!1dxndn(x2−1)n=Pn0(x)
P0(x)=1=1=1
P1(x)=21dxd(−1+x2)=x=x
P2(x)=81dx2d2(−1+x2)2=2−1+23x2=2−1+23x2
P3(x)=481dx3d3(−1+x2)3=−23x+25x3=−23x+25x3
P4(x)=3841dx4d4(−1+x2)4=83−415x2+835x4=83−415x2+835x4
P5(x)=38401dx5d5(−1+x2)5=815x+863x5−435x3=815x+863x5−435x3
P6(x)=460801dx6d6(−1+x2)6=16−5+16105x2+16231x6−16315x4=16−5+16105x2+16231x6−16315x4
ルジャンドル陪多項式
Legendre陪多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します. 展開した結果はMcQuarrie&Simon(1997)等と見比べてください.
Pnm(x)=(1−x2)m/2dxmdmPn(x)=2n1(1−x2)m/2j=0∑⌊2n−m⌋(−1)jj!(n−j)!(n−2j−m)!(2n−2j)!x(n−2j−m)
P00(x)=1=1=1
P10(x)=21dxd(−1+x2)=x=x
P11(x)=(1−x2)21dxd21dxd(−1+x2)=(1−x2)21=(1−x2)21
P20(x)=81dx2d2(−1+x2)2=2−1+23x2=2−1+23x2
P21(x)=(1−x2)21dxd81dx2d2(−1+x2)2=3(1−x2)21x=3(1−x2)21x
P22(x)=(1−x2)dx2d281dx2d2(−1+x2)2=3−3x2=3−3x2
P30(x)=481dx3d3(−1+x2)3=−23x+25x3=−23x+25x3
P31(x)=(1−x2)21dxd481dx3d3(−1+x2)3=−23(1−x2)21+215(1−x2)21x2=−23(1−x2)21+215(1−x2)21x2
P32(x)=(1−x2)dx2d2481dx3d3(−1+x2)3=15x−15x3=15x−15x3
P33(x)=(1−x2)23dx3d3481dx3d3(−1+x2)3=15(1−x2)23=15(1−x2)23
P40(x)=3841dx4d4(−1+x2)4=83−415x2+835x4=83−415x2+835x4
P41(x)=(1−x2)21dxd3841dx4d4(−1+x2)4=−215(1−x2)21x+235(1−x2)21x3=−215(1−x2)21x+235(1−x2)21x3
P42(x)=(1−x2)dx2d23841dx4d4(−1+x2)4=2−15+60x2−2105x4=2−15+60x2−2105x4
P43(x)=(1−x2)23dx3d33841dx4d4(−1+x2)4=105(1−x2)23x=105(1−x2)23x
P44(x)=(1−x2)2dx4d43841dx4d4(−1+x2)4=105(1−x2)2=105(1−x2)2
ラゲール多項式
Laguerre多項式のRodriguesの公式を展開し, 閉形式と比較します. 閉形式についてはLaguerre陪多項式のk=0および一般化Laguerre多項式のα=0の場合を用います. 展開した結果はDLMF等と見比べてください.
Ln(x)=n!exdxndn(e−xxn)=Ln0(x)=Ln(0)(x)
L0(x)=exe−x=1=1=1
L1(x)=exdxdxe−x=1−x=1−x=1−x
L2(x)=21exdx2d2x2e−x=1+21x2−2x=1+21x2−2x=1+21x2−2x
L3(x)=61exdx3d3x3e−x=1+23x2−61x3−3x=1+23x2−61x3−3x=1+23x2−61x3−3x
L4(x)=241exdx4d4x4e−x=1+3x2+241x4−32x3−4x=1+3x2+241x4−32x3−4x=1+3x2+241x4−32x3−4x
L5(x)=1201exdx5d5x5e−x=1+245x4−1201x5−35x3−5x+5x2=1+245x4−1201x5−35x3−5x+5x2=1+245x4−1201x5−35x3−5x+5x2
L6(x)=7201exdx6d6x6e−x=1+7201x6+85x4−201x5−310x3−6x+215x2=1+7201x6+85x4−201x5−310x3−6x+215x2=1+7201x6+85x4−201x5−310x3−6x+215x2
ラゲール陪多項式
Laguerre陪多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します. 閉形式については, いつも通りに関数を定義し, そこに@variables x
で定義したSymbolics.jlの変数を渡せばSymbolics.jlのオブジェクトとして返ってきます. 展開した結果はWikipedia等と見比べてください.
Lnk(x)=dxkdkLn(x)=m=0∑n−k(−1)m+km!(m+k)!(n−m−k)!n!xm=(−1)kLn−k(k)(x)
L00(x)=exe−x=1=1=1
L10(x)=exdxdxe−x=1−x=1−x=1−x
L11(x)=dxdexdxdxe−x=−1=−1=−1
L20(x)=21exdx2d2x2e−x=1+21x2−2x=1+21x2−2x=1+21x2−2x
L21(x)=dxd21exdx2d2x2e−x=−2+x=−2+x=−2+x
L22(x)=dx2d221exdx2d2x2e−x=1=1=1
L30(x)=61exdx3d3x3e−x=1+23x2−61x3−3x=1+23x2−61x3−3x=1+23x2−61x3−3x
L31(x)=dxd61exdx3d3x3e−x=−3−21x2+3x=−3−21x2+3x=−3−21x2+3x
L32(x)=dx2d261exdx3d3x3e−x=3−x=3−x=3−x
L33(x)=dx3d361exdx3d3x3e−x=−1=−1=−1
L40(x)=241exdx4d4x4e−x=1+3x2+241x4−32x3−4x=1+3x2+241x4−32x3−4x=1+3x2+241x4−32x3−4x
L41(x)=dxd241exdx4d4x4e−x=−4−2x2+61x3+6x=−4−2x2+61x3+6x=−4−2x2+61x3+6x
L42(x)=dx2d2241exdx4d4x4e−x=6+21x2−4x=6+21x2−4x=6+21x2−4x
L43(x)=dx3d3241exdx4d4x4e−x=−4+x=−4+x=−4+x
L44(x)=dx4d4241exdx4d4x4e−x=1=1=1
一般化ラゲール多項式
一般化Laguerre多項式のRodriguesの公式(上段)を展開し, 閉形式(下段)と比較します.
Ln(α)(x)=n!x−αexdxndn(xn+αe−x)=k=0∑n(−1)k(n+αn−k)k!xk=k=0∑n(−1)kΓ(α+k+1)Γ(n−k+1)Γ(α+n+1)k!xk
二項係数(ab)はbinomial(a,b)
で計算できます. ただし, αが整数の場合以外はbinomial(a,b)
は計算できないので, (ab)=b!(a−b)!a!=Γ(b+1)Γ(a−b+1)Γ(a+1) を利用し, 階乗表示を経由してガンマ関数から計算することにします.(そもそも, αが整数なら陪関数との対応関係を使えば十分です. Morse振動子の波動関数ではαは非整数なので, 非整数に対応しておく必要があります.) なお, Laguerre多項式の定義によっては, 一般化Laguerre多項式はLaguerre多項式の一般化にならないことに注意してください. α=0の場合はLn(0)(x)=n!exdxndn(xne−x)となり, Laguerre多項式の定義をLn(x)=exdxndn(xne−x)とした場合とは一致しません. このノートではLaguerre多項式をLn(x)=n!exdxndn(xne−x)と定義したので, 一般化になっています.
L0(0)(x)=exe−x=1=1=1
L1(0)(x)=exdxdxe−x=1−x=1−x=1−x
L1(1)(x)=xexdxdx2e−x=2−x=2−x=2−x
L2(0)(x)=21exdx2d2x2e−x=1+21x2−2x=1+21x2−2x=1+21x2−2x
L2(1)(x)=x21exdx2d2x3e−x=3+21x2−3x=3+21x2−3x=3+21x2−3x
L2(2)(x)=x221exdx2d2x4e−x=6+21x2−4x=6+21x2−4x=6+21x2−4x
L3(0)(x)=61exdx3d3x3e−x=1+23x2−61x3−3x=1+23x2−61x3−3x=1+23x2−61x3−3x
L3(1)(x)=x61exdx3d3x4e−x=4+2x2−61x3−6x=4+2x2−61x3−6x=4+2x2−61x3−6x
L3(2)(x)=x261exdx3d3x5e−x=10+25x2−61x3−10x=10+25x2−61x3−10x=10+25x2−61x3−10x
L3(3)(x)=x361exdx3d3x6e−x=20+3x2−61x3−15x=20+3x2−61x3−15x=20+3x2−61x3−15x
L4(0)(x)=241exdx4d4x4e−x=1+3x2+241x4−32x3−4x=1+3x2+241x4−32x3−4x=1+3x2+241x4−32x3−4x
L4(1)(x)=x241exdx4d4x5e−x=5+5x2+241x4−65x3−10x=5+5x2+241x4−65x3−10x=5+5x2+241x4−65x3−10x
L4(2)(x)=x2241exdx4d4x6e−x=15+215x2+241x4−x3−20x=15+215x2+241x4−x3−20x=15+215x2+241x4−x3−20x
L4(3)(x)=x3241exdx4d4x7e−x=35+221x2+241x4−67x3−35x=35+221x2+241x4−67x3−35x=35+221x2+241x4−67x3−35x
L4(4)(x)=x4241exdx4d4x8e−x=70+14x2+241x4−34x3−56x=70+14x2+241x4−34x3−56x=70+14x2+241x4−34x3−56x
まとめ
おかげで各種の閉形式についてRodriguesの公式との整合性が簡単に確認できるようになりました. Symbolics.jlはまだ発展途中のパッケージですが, かなり活躍してくれましたね. 自分でルールを追加して機能を拡張できるようなので, これから比較にならないほど発展すると思われます. 現状のシステムに強いて注文を付けるとすれば, 多項式の各項を次数の順番でソートできると良いですね. Issueに上がっているので, 対応されるものと思われます.
参考文献
https://symbolics.juliasymbolics.org/dev/manual/derivatives/
文献 |
|
|
|
|
|
Ln(α)(x) |
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