フロベニウス法でベッセル
第１種 整数m次ベッセル関数 Jm(x)
１。フロベニウス法とは
このワークシートはMath by Codeの一部です。
前回まで、べき級数を利用した微分方程式の解法を見てきた。
今回は、べき級数の拡張版であるフロベニウス法でベッセルの微分方程式を解こう。
＜フロベニウス法はべき級数法の拡張版＞
p(x),q(x)がｘ＝０で解析的なとき、
微分方程式B=の係数は解析的ではない。
だから、ｘ＝０は確定特異点という。
解の１つはy=Σamxm+r=xr(a0+a1x+a2x2+....)=a0xr+a1xr+1+a2xr+2+.....（a0≠0)ｒは複素数。
これとは1次独立な解をもつ。
係数が解析的でなくても、解くことができる、べき級数法の拡張版といえる。
まず、級数の形にしてみよう。
y'=ra0xr-1+(r+1)a1xr+(r+2)a2xr+1+....
y"=r(r-1)a0xr-2+(r+1)ra1xr-1+(r+2)(r+1)a2xr+....
p(x)=p0+p1x+p2x2+.....
q(x)=q0+q1x+q2x2+.....
Bの左辺にx2をかけて展開しよう。
x2(r(r-1)a0xr-2+(r+1)ra1xr-1+(r+2)(r+1)a2xr+....)
+x(p0+p1x+p2x2+....)(ra0xr-1+(r+1)a1xr+(r+2)a2xr+1+....)
+(q0+q1x+q2x2+.....)(a0xr+a1xr+1+a2xr+2+.....)=0
xrの係数の和は、r(r-1)a0＋p0ra0＋q0a0
=[r(r-1)＋p0r＋q0]a0=0 a0≠0だから、r(r-1)＋p0r＋q0=0
rについて整理する。
r2+(p0-1)r+q0=0
このｒについての方程式を決定方程式という。
ｒの解が重複解か？重複解でないときは２解の差が整数かどうかで対処の仕方が変わる。
２．ガンマ関数
ベッセルの微分方程式は熱、光、惑星の運動などの宇宙工学など幅広く使われる大切なものだといわれる。その解でガンマ関数という表現が登場するので、そこから調べてみよう。
＜ガンマ関数の定義と性質とは＞
ガンマ関数は階乗関数の一般化だとよく言われる。どういうことだろうか。
ガンマ関数は正の実数xについて、Γ(x)=∫tx-1e-t dt[0,∞]と定義される。
・x=1のときΓ(1)=∫t1-1e-t dt=∫e-t dt=[-e-t][0,∞]=0-(-1)=1
・部分積分∫Fg=FG-∫fGで、F=tx,g=e-t とするとき、gを積分するとG=-e-t、Fを微分するf=xtx-1
x=x+1のときΓ(x+1)=∫txe-t dt=[-txe-t ][0,∞]-∫-xtx-1e-tdt = x∫tx-1e-tdt=xΓ(x)
xが正の実数のとき、Γ(x+1)=xΓ(x)言い換えると、Γ(x)=(x-1)Γ(x-1)
・とくに、xが整数のときも成り立つから、Γ(2)=1Γ(1), Γ(3)=2Γ(2),....。
つまり、nが正の整数のとき、Γ(n)=(n-1)!（階乗関数）
・x=1/2のときt=y2とおくと、dt/dy = 2yからdt=2y dyで、 y=t1/2
Γ(1/2)=∫t1/2-1e-t dt[0,∞]=∫t-1/2e-t dt[0,∞]=∫y-1e-y^2 (2y dy)[0,∞] =2∫e-y^2dy[0,∞]＝2Iとおく。
I＝∫e-y^2dy[0,∞]は、J=∫e-y^2dy[0,a]=∫e-x^2dx[0,a]とおくと、
J2=(∫e-x^2dx[0,a])(∫e-y^2dy[0,a])=∫∫D(a)e-(x^2+y^2)dxdy、 D(a)={(x,y)|x,yが0以上a以下の正方形}とする。
原点中心の第一象限の半径ｐの領域をE(p)とすると、E(a)≦D(a)≦E(√2 a)。
正方形が大小の四分円の弧にはさまれる。それぞれの領域での2重積分 e-(x^2+y^2)もこの順になる。
ここで、極座標変換x=rcosθ,y=rsinθをして、2重積分の積分変数を置き換えよう。
ヤコビアン|{xr xθ}{yr yθ}|=xr * yθ -xθ* yr=cosθ r cosθ - (-rsinθ sinθ) =r(cos2θ+sin2θ)=rとなるね。
だから、関数f(x,y)=e-(x^2+y^2)のx,yに極座標を代入して、dxdyを r dθdrにおきかえると
∫∫E(a)e-(x^2+y^2)dxdy=∫dθ[π/2, 0] ∫e-r^2 rdr[a,0]=π/4(1-e-a^2 )はa→∞でπ/4に収束する。
同様に∫∫E(√2a)e-(x^2+y^2)dxdyもπ/4に収束する。
だから、挟み撃ちで、J2もπ/4に収束するから、Iは√(π/4)=1/2√πとできるね。
結局、Γ(1/2)=2・1/2√π=√π。
（例）
Γ(2.5)=1.5Γ(1.5)=1.5・0.5Γ(0.5)=3/4√π。
x=[1,2,3,4,5]のとき、Γ(x)=[1,1,2,3,4]
x=[0.5 , 1.5, 2.5, 3.5, 4.5, 5.5 ] =√π [1, 1/2, 3/4,15/8, 105/16, 1155/32 ]
ガンマ関数 y=Γ(x)
３．ベッセルの微分方程式
＜第１種関数＞
mを非負の定数とするとき、
p(x)=x, q(x)=x2-m2が解析的で、ベッセルの微分方程式ｘ2y"+xy'+(x2-m2)y=0を解きたい。
微分方程式B=の係数は解析的ではない。
解の１つをy=xrΣa_0xn=xr(a0+a1x+a2x2+....)=a0xr+a1xr+1+a2xr+2+.....（a0≠0) ｒは複素数としよう。
y'=ra0xr-1+(r+1)a1xr+(r+2)a2xr+1+....
y"=r(r-1)a0xr-2+(r+1)ra1xr-1+(r+2)(r+1)a2xr+....
p(x)=x＝p0+p1x+p2x2+.....だから、p1=1で、他は０。
q(x)=x2-m2=q0+q1x+q2x2+.....だから、q0=-m2, q2=1で、他は０。
x2(r(r-1)a0xr-2+(r+1)ra1xr-1+(r+2)(r+1)a2xr+....)
+x(p0+p1x+p2x2+....)(ra0xr-1+(r+1)a1xr+(r+2)a2xr+1+....)
+(q0+q1x+q2x2+.....)(a0xr+a1xr+1+a2xr+2+.....)=0
この係数を具体化すると、
x2(r(r-1)a0xr-2+(r+1)ra1xr-1+(r+2)(r+1)a2xr+....)
+x(ra0xr-1+(r+1)a1xr+(r+2)a2xr+1+....)
+(x2-m2)(a0xr+a1xr+1+a2xr+2+.....)=0
xrの係数の和は、r(r-1)a0 + ra0-m2a0＝a0(r2-m2)=0から、r=±mの2解をもつ。
帰結C1: xr+1の係数の和は、(r+1)r a1＋(r+1)a1-m2a1＝[(r+1)2 -m2]a1=0
帰結C2:xr+2の係数の和は、(r+2)(r+1)a2＋(r+2)a2-m2a2+a0=0から、[(r+2)2 -m2]a2+a0=0
帰結Ck:xr+kの係数の和は、(r+k)(r+k-1)ak＋(r+k)ak-m2ak+ak-2=0から、[(r +k+m)(r+k-m)]ak+ak-2=0
・r=m（非整数）に対しては、
帰結C１から、mが非負ならa1=0。帰結C2から連鎖的に、奇数番目の係数a_odd=0
a0=Aとおくと、帰結Ckから、ak=-ak-2/[(r +k+m)(r+k-m)]=-ak-2/[(k+2m)k]
ここで、k=2sとおくと、a2s=-a2s-2/[(2s+2m)2s]=-a2s-2/[4(m+s)s]
Aに対して,sに1,2,3,4,5,....pを入れてこの漸化式を連続適用すると、
基本解y=xrΣa_nxnの奇数番目の係数は０で、
偶数番目の係数はa_even=a2p=(-1)pA/[22p p!(m+1)(m+2).....(m+p)]
A=1/[2^mΓ(m+1)]とおくと、a_even=a2p=(-1)pA/[22p+m p!Γ(m+p+1)]（mは非負の実数）
解y=y1=xmΣa_nxn = [p for 0 to ∞] を
m次の第１種ベッセル関数Jm(x)という。
mが整数でないなら、r=-mとおいて同様にして得られる基本解y2=J(-m)(x)はy1と一次独立になる。
y1とy2の１次結合が一般解となる。
・r=m=n(整数）に対して
A=1/[2^n n!]とおくと、a_even=(-1)p/[22p+n p!(n+p)!] (ｎは整数）
r=nとしたとき、xnにp=0から∞の級数をかけて、
解y=xrΣa_nxn = [p for 0 to ∞ ] をｎ次の第１種ベッセル関数Jn(x)という。
rが整数だとJ-n(x)は解Jn(x)に従属します。そのために、ノイマン関数Nn(x)=[Jn(x)cosπn-J-n(x)]/sinπnを
もう一つの基本解として、その1次結合を一般解とします。
このノイマン関数のことを第２種ベッセル関数ということもありますね。
（例）
J0(x)=
（例）
J1(x)=
（例）
r=m=1/2を入れて、
J 1/2 (x) = [p for 0 to ∞ ]
分母のp!は１からｐまでのｐ個の積で、
分母のΓ(p+1+1/2)はΓ(1/2)に1/2,3/2,....(2p+1)のｐ+1個の積をかけたもの。
(x/2)^(2p+1/2) はxの2p+1個の積/（2のp個の積と２のp+1個の積）×√2/xと等しい。
２のp個の積をp!のp個の要素に１対１につけると、2から2pまでのp個の偶数の積
２のp+1個の積を1/2,3/2,....(2p+1)のｐ+1個の要素に１対１につけると、
1から2p+1までのp+1個の奇数の積。
これらの積は、(-1)^p/[ Γ(1/2)×偶数の積×奇数の積] ｘ^(2p+1)*×√2/x
pがΣの変数だから、それ以外を定数扱いにしてΣの前にだそう。
J 1/2 (x) = [p for 0 to ∞ ]
Σの中はsinxのマクローリン展開と同じになったね。
J 1/2 (x) =√(２/πｘ) sin x
（例）
同様にして、...............................
J -1/2 (x) =√(２/πｘ) cos x
第１種1/2次、-1/2次のベッセル関数 J1/2(x),J-1/2(x)
４．実装
質問：第１種ベッセル関数をコードで実装するにはどうしたらよいですか。
やはり、科学計算のライブラリの多いpythonがおすすめです。
ベッセル関数は特殊関数という種類なのでscipy.specialを使います。
scipy.specialのjv(n,x)関数が、第n次の第１種ベッセル関数だね。
figで絵を枠を用意し、fig.add_subplot(111)が枠にプロットしょう。
ｘの値を[-10,10]区間を70等分くらいにすると線がなめらかになるでしょう。
ｙの値がjv(n,x)です。この(x,y)を使い、ax.plot(x, jv(v,x), label=...)でグラフをかきます。
#[IN]Python==============================
# 第1種ベッセル関数
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jv
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(111)
x = np.linspace(-10, 10, 70)
for n in range(5):
ax.plot(x, jv(n, x), label=f"J{n}(x)")
ax.legend()
#[OUT]=================================
