Blog

ベッセルの微分方程式

December 19, 2019      

今回は，円筒座標系において熱伝導方程式を解くための前準備として，1を参考にベッセルの微分方程式についてまとめておく．

確定特異点

d2udz2+P(z)dudz+Q(z)u=0\begin{equation} % \label{eq:1983Terasawa_754} \frac{d^2 u}{dz^2} + P(z) \frac{du}{dz} + Q(z) u = 0 \end{equation}

(1)において，P(z)P(z)またはQ(z)Q(z)がz=az=aにおいて極(pole)をもち，かつその点において(z−a)P(z)(z-a)P(z)および(z−a)2Q(z)(z-a)^2 Q(z)は正則であるとする．このとき極z=az=aを(1)の確定特異点(Regular Singular Point)と呼ぶ．（複素平面上ある変域内のいずれの点においてもf′(z)f'(z)が存在する場合には、関数f(z)f(z)はその変域内において正則であるといい、このような関数を正則関数という。また、ある変域内で正則な関数は、その辺域内において何回微分しても正則である。これをGoursatの定理という）

p(z−a)=(z−a)P(z)p(z-a)=(z-a)P(z)およびq(z−a)=(z−a)2Q(z)q(z-a)=(z-a)^2 Q(z)とおくと，(1)は次のように書き換えられる．

d2udz2+p(z−a)z−adudz+q(z−a)(z−a)2u=0\begin{equation} \frac{d^2 u}{dz^2} + \frac{p(z-a)}{z-a} \frac{du}{dz} + \frac{q(z-a)}{(z-a)^2} u = 0 \end{equation}

また，p(z−a)p(z-a)およびq(z−a)q(z-a)はz=az=aの周りで正則なので，以下のようにTaylor展開できる．

p(z−a)=p0+(z−a)p1+(z−a)2p2+...q(z−a)=q0+(z−a)q1+(z−a)2q2+...\begin{align*} p(z-a) &= p_0 + (z-a) p_1 + (z-a)^2 p_2 + ... \\ q(z-a) &= q_0 + (z-a) q_1 + (z-a)^2 q_2 + ... \end{align*}

(1)の解として以下の無限級数を仮定し、λ\lambdaおよびAnA_nを定める．

u=(z−a)λ{A0+∑n=1∞An(z−a)n}\begin{equation} u=(z-a)^{\lambda}\left\{ A_0 + \sum_{n=1}^{\infty} A_n (z-a)^n \right\} \end{equation}

ただし，uuの一階微分，二階微分はそれぞれ以下のように表される．

dudz=λ(z−a)λ−1{A0+∑n=1∞An(z−a)n}+(z−a)λ{∑n=1∞nAn(z−a)n−1}=(z−a)λ−1{λA0+∑n=1∞(λ+n)An(z−a)n}d2udz2=(λ−1)(z−a)λ−2{λA0+∑n=1∞(λ+n)An(z−a)n}+(z−a)λ−1{∑n=1∞n(λ+n)An(z−a)n−1}=(z−a)λ−2{(λ−1)λA0+∑n=1∞(λ+n−1)(λ+n)An(z−a)n}\begin{align*} \frac{du}{dz} &= \lambda (z-a)^{\lambda-1} \left\{ A_0 + \sum_{n=1}^{\infty} A_n (z-a)^n \right\} + (z-a)^{\lambda} \left\{ \sum_{n=1}^{\infty} n A_n (z-a)^{n-1} \right\} \\ &= (z - a)^{\lambda-1} \left\{ \lambda A_0 + \sum_{n=1}^{\infty} (\lambda + n) A_n (z-a)^n \right\} \\ \frac{d^2 u}{dz^2} &= (\lambda - 1) (z-a)^{\lambda-2} \left\{ \lambda A_0 + \sum_{n=1}^{\infty} (\lambda + n) A_n (z-a)^n \right\} + (z - a)^{\lambda-1} \left\{ \sum_{n=1}^{\infty} n (\lambda + n) A_n (z-a)^{n-1} \right\} \\ &= (z-a)^{\lambda-2} \left\{ (\lambda - 1) \lambda A_0 + \sum_{n=1}^{\infty} (\lambda + n - 1)(\lambda + n) A_n (z-a)^n \right\} \end{align*}
(z−a)λ−2{(λ−1)λA0+∑n=1∞(λ+n−1)(λ+n)An(z−a)n}    +(z−a)λ−2p(z−a){λA0+∑n=1∞(λ+n)An(z−a)n} +(z−a)λ−2q(z−a){A0+∑n=1∞An(z−a)n}=0\begin{align*} &(z-a)^{\lambda-2} \left\{ (\lambda - 1) \lambda A_0 + \sum_{n=1}^{\infty} (\lambda + n - 1)(\lambda + n) A_n (z-a)^n \right\} \\ &~~~~+ (z - a)^{\lambda-2} p(z-a) \left\{ \lambda A_0 + \sum_{n=1}^{\infty} (\lambda + n) A_n (z-a)^n \right\}  + (z-a)^{\lambda-2} q(z-a) \left\{ A_0 + \sum_{n=1}^{\infty} A_n (z-a)^n \right\} = 0 \end{align*}

(z−a)(z-a)について低次の項から，整理すると以下のように表される．

A0{(λ−1)λ+p0λ+q0}=0A0(p1λ+q1)+A1{λ(λ+1)+p0(λ+1)+q0}=0A0(p2λ+q2)+A1{p1(λ+1)+q1}+A2{(λ+1)(λ+2)+p0(λ+2)+q0}=0⋮\begin{gather*} A_0 \left\{ (\lambda - 1)\lambda + p_0 \lambda + q_0 \right\} = 0 \\ A_0 ( p_1 \lambda + q_1 ) + A_1 \left\{ \lambda(\lambda+1) + p_0 (\lambda+1) + q_0 \right\} = 0 \\ A_0 ( p_2 \lambda + q_2 ) + A_1 \left\{ p_1 (\lambda+1) + q_1 \right\} + A_2 \left\{ (\lambda+1)(\lambda+2) + p_0 (\lambda+2) + q_0 \right\}= 0 \\ \vdots \end{gather*}

これは一般に次のように表すことができる．

A0ϕ(λ)=0A1ϕ(λ+1)+A0ϕ1(λ)=0A2ϕ(λ+2)+A1ϕ1(λ+1)+A0ϕ2(λ)=0⋮Anϕ(λ+n)+⋯+A0ϕn(λ)⋮\begin{gather*} A_0 \phi(\lambda) = 0 \\ A_1 \phi(\lambda+1) + A_0 \phi_1 (\lambda) = 0 \\ A_2 \phi(\lambda+2) + A_1 \phi_1 (\lambda+1) + A_0 \phi_2 (\lambda) = 0 \\ \vdots \\ A_n \phi(\lambda+n) + \cdots + A_0 \phi_n (\lambda) \\ \vdots \end{gather*}
where{ϕ(λ)=λ2+(p0−1)λ+q0ϕn(λ)=λpn+qn(n=1,2,3,⋯ )\begin{gather} % \label{eq:1983Terasawa_757} \mathrm{where} \hspace{10pt} \left\{ \begin{array} {l} \phi(\lambda) = \lambda^2 + (p_0 - 1)\lambda + q_0 \\ \phi_n(\lambda) = \lambda p_n + q_n \hspace{10pt} (n=1,2,3,\cdots) \end{array} \right. \end{gather}

(z−a)(z-a)についてゼロ次の項から，確定特異点z=az=aにおける指数を定める式が得られる．

ϕ(λ)≡λ2+(p0−1)λ+q0=0\begin{equation} % \label{eq:1983Terasawa_758} \phi(\lambda) \equiv \lambda^2 + (p_0 - 1)\lambda + q_0 = 0 \end{equation}

ここでの目的は，(1)の特解を得ることなのでA0=1A_0=1とおいてよい．ϕ(λ+n)\phi(\lambda + n)がゼロでない限りそれぞれのλ\lambdaに対して，A1,A2,...A_1, A_2, ...が一義的に定まり，以下の特解が得られる．一般解はこれらの線形結合として表される．

u1=(z−a)λ1{1+∑n=1∞An1(z−a)n}u2=(z−a)λ2{1+∑n=1∞An2(z−a)n}\begin{align*} u_1 &= (z-a)^{\lambda_1} \left\{ 1+ \sum_{n=1}^{\infty} A_n^{1} (z-a)^n \right\} \\ u_2 &= (z-a)^{\lambda_2} \left\{ 1+ \sum_{n=1}^{\infty} A_n^{2} (z-a)^n \right\} \end{align*}

ベッセルの微分方程式

d2ydz2+1zdydz+(1−n2z2)y=0\begin{equation} % \label{eq:1983Terasawa_771} \frac{d^2 y}{dz^2} + \frac{1}{z} \frac{dy}{dz} + \left( 1 - \frac{n^2}{z^2} \right) y = 0 \end{equation}

ベッセルの微分方程式は(6)のように表され，x=0x=0はその確定特異点である．p(x)=1p(x)=1およびq(x)=x2−n2q(x)=x^2-n^2とおけて，(5)よりλ\lambdaが以下の式より得られる．

ϕ(λ)≡λ2−n2=0\begin{equation} \phi(\lambda) \equiv \lambda^2 - n^2 = 0 \end{equation}

λ1=n\lambda_1=nのとき，以下の形の特解が存在する．

y=xn(1+A1x+⋯+Amxm+⋯ )\begin{equation} y = x^n (1 + A_1 x + \cdots + A_m x^m + \cdots) \end{equation}

(4)よりϕ1=0,ϕ2=1,ϕ3=0,ϕ4=0,⋯\phi_1=0, \phi_2=1, \phi_3=0, \phi_4=0, \cdotsなので，係数A0,A1,⋯A_0, A_1, \cdotsは以下のように求められる．

Am=−Am−2ϕ(λ+m)=−1ϕ(λ+m)−Ai−4ϕ(λ+m−2)=−1ϕ(λ+m)−1ϕ(λ+m−2)⋯−A0ϕ(λ+2)\begin{equation*} A_m = - \frac{A_{m-2}}{\phi(\lambda + m)} = \frac{-1}{\phi(\lambda + m)} \frac{-A_{i-4}}{\phi(\lambda + m - 2) } = \frac{-1}{\phi(\lambda + m)} \frac{-1}{\phi(\lambda + m - 2) } \cdots \frac{-A_0}{\phi(\lambda + 2) } \end{equation*}
A2m=(−1)mm!22m(n+1)(n+2)⋯(n+m)\begin{equation} A_{2m} = \frac{(-1)^m}{m! 2^{2m} (n+1)(n+2) \cdots (n+m)} \end{equation}

このように得られた解に定数1/2nΓ(n+1)1/2^n \Gamma(n+1)をかけたものを，Bessel Function of the first kindと呼ぶ．λ2=−n\lambda_2=-nの解は，nnの符号を変えれば得られて，それぞれ以下のように表される．

Jn(x)=∑m=0∞(−1)mm!Γ(n+m+1)(x2)n+2mJ−n(x)=∑m=0∞(−1)mm!Γ(−n+m+1)(x2)−n+2m\begin{align} J_n(x) &= \sum_{m=0}^{\infty} \frac{(-1)^m}{m! \Gamma (n+m+1)} \left( \frac{x}{2} \right)^{n+2m} \\ J_{-n}(x) &= \sum_{m=0}^{\infty} \frac{(-1)^m}{m! \Gamma (-n+m+1)} \left( \frac{x}{2} \right)^{-n+2m} \end{align}

nnが整数でない場合には，JnJ_nとJ−nJ_{-n}は別の関数となるので，(6)の一般解は次のように表される．

y=C1Jn(x)+C2J−n(x)\begin{equation} y = C_1 J_n(x) + C_2 J_{-n}(x) \end{equation}

変形ベッセル関数

(6)の変数としてizizをとると，以下のように変形できる．

d2v(iz)d(iz)2+1izdv(iz)d(iz)+(1−n2(iz)2)v(iz)=0−d2v(iz)dz2−1zdv(iz)dz+(1+n2z2)v(iz)=0\begin{align*} \frac{d^2 v(iz)}{d(iz)^2} + \frac{1}{iz} \frac{d v(iz)}{d(iz)} + \left( 1 - \frac{n^2}{(iz)^2} \right) v(iz) = 0 \\ - \frac{d^2 v(iz)}{dz^2} - \frac{1}{z} \frac{d v(iz)}{dz} + \left( 1 + \frac{n^2}{z^2} \right) v(iz) = 0 \end{align*}

ただし，iiを含む微分は次のように表すことができる．

d(iz)dzdd(iz)=ddz,   dd(iz)=−iddz,   d2d(iz)2=−d2dz2\begin{gather*} \frac{d(iz)}{dz} \frac{d}{d(iz)} = \frac{d}{dz}, ~~~ \frac{d}{d(iz)} = - i \frac{d}{dz},~~~ \frac{d^2}{d(iz)^2} = - \frac{d^2}{dz^2} \end{gather*}

u(z)=v(iz)u(z)=v(iz)を新しくとると，Modified Bessel’s Differential Equationが得られる．

d2udz2+1zdudz−(1+n2z2)u=0\begin{equation} % \label{eq:1983Terasawa_1163} \frac{d^2 u}{dz^2} + \frac{1}{z} \frac{d u}{dz} - \left( 1 + \frac{n^2}{z^2} \right) u = 0 \end{equation}

Bessel Functionにizizを代入したものが特解となるが，一般的にはe−12νπiJν(iz)e^{-\frac{1}{2}\nu \pi i} J_{\nu} (iz)を使用する．

Iν=∑m=0∞1m!Γ(ν+m+1)(z2)ν+2m\begin{equation} % \label{eq:1983Terasawa_1164} I_{\nu} = \sum_{m=0}^{\infty} \frac{1}{m! \Gamma (\nu+m+1)} \left( \frac{z}{2} \right)^{\nu+2m} \end{equation}

また，e12νπiJν(iz)e^{\frac{1}{2}\nu \pi i} J_{\nu} (iz)も(14)の解となる．

I−ν=∑m=0∞1m!Γ(−ν+m+1)(z2)−ν+2m\begin{equation} I_{-\nu} = \sum_{m=0}^{\infty} \frac{1}{m! \Gamma (-\nu+m+1)} \left( \frac{z}{2} \right)^{-\nu+2m} \end{equation}

応用上は，e12νπiJν(iz)e^{\frac{1}{2}\nu \pi i} J_{\nu} (iz)を(13)の特解としては採用せず，以下の関数を用いる．

Kν(z)=π2I−ν(z)−Iν(z)sin⁡νπ\begin{equation} % \label{eq:1983Terasawa_1165} K_{\nu}(z) = \frac{\pi}{2} \frac{I_{-\nu}(z) - I_{\nu}(z)}{\sin \nu \pi} \end{equation}

プログラム上で第1種変形ベッセル関数(14)，第2種変形ベッセル関数(16)を使用する場合，例えば以下のような関数を用いることが出来る．次数ν\nuは整数でなくてもよいが，実数でなくてはならない．変数zzは複素数をとる．

modified Bessel function of the 1st kind,  MATLAB: besseli(ν,z),  Python: scipy.special.iv(ν,z)modified Bessel function of the 2nd kind,  MATLAB: besselk(ν,z),  Python: scipy.special.kv(ν,z)\begin{align*} \mathrm{modified~Bessel~function~of~the~1st~kind,} ~~ &\mathrm{MATLAB:}\ besseli(\nu,z), ~~\mathrm{Python:}\ scipy.special.iv(\nu, z) \\ \mathrm{modified~Bessel~function~of~the~2nd~kind,} ~~ &\mathrm{MATLAB:}\ besselk(\nu,z), ~~\mathrm{Python:}\ scipy.special.kv(\nu, z) \end{align*}
  1. 寺沢寛一，“自然科学者のための数学概論 増訂版” ，岩波書店，1983↩