レーン=エムデン方程式は以下の形で表される
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\qquad\cdots\,\)(1)
境界条件は、\(\xi\to 0\,\)のとき\(\,\theta=1\,,\,\dfrac{d\theta}{d\xi}=0\,\)である。
［１］中心近傍での扱い
この方程式は\(\xi=0\,\)が確定特異点になっている。この場合、特異点近傍でローラン級数展開を考える。\(\xi=0\,\)の周りの形式解を以下のようにする。
\(\quad\theta(\xi)=1+a_2\xi^2+a_4\xi^4+\cdots\qquad\)(2)
式(2)を式(1)の左辺に代入する。まず、\(\xi\,\)で微分すると
\(\quad\dfrac{d\theta}{d\xi}=2a_2\xi+4a_4\xi^3+\cdots\qquad\)(3)
さらに\(\xi^2\,\)をかけて微分すると
\(\quad\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=6a_2\xi^2+20a_4\xi^4+\cdots\qquad\)(4)
\(\therefore\,\,\)左辺\(\,\,=6a_2+20a_4\xi^2+\cdots\qquad\)(5)
右辺は、まず\(B_2=a_2\xi^2\,,\,B_4=a_4\xi^4\,,\cdots\,\)として、
\(\quad(1+B_2+B_4+\cdots)^N\,\)の展開を考える。
最初の\(\,1\,\)の項は\(\,1^NB_2^0B_4^0\cdots\,\)であるので、
係数は\(\,\dfrac{N!}{N!0!0!\cdots}\,\)となる。
次の項は\(1^{N-1}B_2^1B_4^0\cdots\,\)で、
係数は\(\,\dfrac{N!}{(N-1)!1!0!\cdots}\,\)となり、
次の\(1^{N-1}B_2^0B_4^1B_6^0\cdots\,\)の係数は
\(\,\dfrac{N!}{(N-1)!0!1!0!\cdots}\,\)となる。因って
\(\quad(1+B_2+B_4+\cdots)^N=1+NB_2+NB_4+\cdots\qquad\)(6)
右辺\(\,=-\theta^N=-1-Na_2\xi^2-Na_4\xi^4-\cdots\,\)となるので
\(\,\xi^m\,\)の係数を比較すると
\(\quad 6a_2=-1\,,\qquad 20a_2=-Na_4\qquad\)となり
\(\quad a_2=-\dfrac{1}{6}\,,\qquad a_4=\dfrac{N}{120}\quad\)を得る。
よって中心近傍での解は
\(\quad \theta(\xi)\simeq 1-\dfrac{1}{6}\xi^2+\dfrac{N}{120}\xi^4\quad\cdots\,\)(7)
と表すことが出来る。
［２］数値計算
この方程式は(その2)で求めた解析解以外の解は数値計算に因らざる得ない。
\(\bigstar\,\)常微分方程式の解法
\(\quad\dfrac{dy_i}{dx}=f_i[x,y_i(x)]\quad i=1,2,\dots\quad\)(8)
を初期値問題として数値的に解く。
まず、\(x=x_0\,\)において初期値\(\,y_i(x_0)=y_{i0}\,\)が与えられたとき、\(x=x_0+h\,\)での\(\,y_i\,\)の値\(\,y_i(x_0+h)\,\)をもとめ、次のステップでその\(\,(x,y_i)\,\)を新たな初期値として処理を繰り返す。
\(x_0\,\)と\(\,x_0+h\,\)の間の中間の\(\,x\,\)における導関数、\(f(x,y)\,\)の関数値を何回か計算し、それらの適当な平均を用いることによって、高次のTaylor展開を達成するのが Runge-Kutta 法である。
\(\bigstar\,\)Runge-Kutta 法
\(\quad k_1=hf(x_0,y_0)\qquad\qquad\qquad\cdots\,\)(9)
\(\quad k_2=hf(x_0+\frac{h}{2},y_0+\frac{k_1}{2})\qquad\cdots\,\)(10)
\(\quad k_3=hf(x_0+\frac{h}{2},y_0+\frac{k_2}{2})\qquad\cdots\,\)(11)
\(\quad k_4=hf(x_0+h,y_0+k_3)\qquad\cdots\,\)(12)
\(\quad y(x_0+h)=y(x_0)+\frac{1}{6}[K_1+2(K_2+k_3)+k_4]\,\cdots\,\)(13)
まず点(\(\,x_i,y_i\,\))での微係数から\(\,y\,\)の増分\(\,k_1\,\)を推定する(式9)。次に(\(\,x_i,y_i\,\))と今得た点(\(\,x_i+\Delta x,y_i+k_1\,)\)の中間点での微係数から、これに平行な破線で示す直線によって増分\(\,k_2\,\)の推定を行う(式10)。これより第３の推定点\((\,x_i+\Delta x,y_i+k_3\,)\)を得る。第４の推定点まで求めたところで、これら４つの増分推定を\(\,1:2:2:1\,\)の重みを付けて平均化し、これを最終的な増分推定として\(\,y_{i+1}\,\)を求める(式13)。この式は Taylor展開の４次の精度を持つ。
\(\bigstar\,\) 初期値とプログラミング
解くべき式は(1)式だが、\(\xi=x\,,\,\theta=y\,\)と書き換えて見やすくする。
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-y^N\quad x\!=\!0\quad y\!=\!0\,,\,\dfrac{dy}{dx}=0\,\cdots\,\)(14)
この常微分方程式を連立１階微分方程式に変換する。
\(y1=y\,,\quad y2=x^2\dfrac{dy}{dx}\,\)とおくと
\(\quad \dfrac{dy1}{dx}=\dfrac{y2}{x^2}\qquad\cdots\,\)(15)
\(\quad \dfrac{dy2}{dx}=-x^2y1^N\qquad\cdots\,\)(16)
初期値としては\(\,xi=x=0\,\)でなく、演算の刻み\(\,\xi=x=h\,\)とし、\(y1\,\)は式(7)を使い
\(\quad y1=1-\dfrac{1}{6}h^2+\dfrac{N}{120}h^4\quad\cdots\,\)(17)
\(y2\,\)は、\(\xi^2\dfrac{d\theta}{d\xi}=y2\quad\)としたので
\(\quad y2=h^2\left(-\dfrac{1}{3}h+\dfrac{N}{30}h^3\right)=-\dfrac{1}{3}h^3+\dfrac{N}{120}h^4\quad\cdots\,\)(18)
で設定する。
計算は中心\(\,\,\xi=x=h\,\,\)から始めて、密度係数がゼロ\(\,\,\,\theta=y1=0\,\,\)となる表面までを積分する。また\(\,\,h=1.0\times10^{-6}\,\,\)とした。
［３］数値計算による結果
\(N\,\)を\(0\sim 4.5\,\)まで計算した結果をグラフと表に示す。
この表は、\(N=n\,\)でガス球の表面（圧力ゼロ\(\,\theta=0\,\)）の時の\(\,\xi_n\,\)とそこでの\(\,-1/\xi\,(d\theta/d\xi)\,\)で、中心密度\(\rho_c\,\)を計算するためのパラメータである。
| Nの値 | \(\xi_n\) | \(\left(-\frac{1}{\xi}\frac{d\theta}{d\xi}\right)_{\xi=\xi_n}\) |
| \(0.0\) | \(2.549\) | \(3.291\times10^{-1}\) |
| \(0.5\) | \(2.753\) | \(1.816\times10^{-1}\) |
| \(1.0\) | \(3.142\) | \(1.013\times10^{-1}\) |
| \(1.5\) | \(3.654\) | \(5.564\times10^{-2}\) |
| \(2.0\) | \(4.353\) | \(2.923\times10^{-2}\) |
| \(2.5\) | \(5.355\) | \(1.424\times10^{-2}\) |
| \(3.0\) | \(6.897\) | \(6.152\times10^{-3}\) |
| \(3.5\) | \(9.536\) | \(2.180\times10^{-3}\) |
| \(4.0\) | \(14.972\) | \(5.356\times10^{-4}\) |
| \(4.5\) | \(31.836\) | \(5.386\times10^{-5}\) |
\(\bigstar\,\)結果の見方
\(\xi\,\)は \(\,R=\alpha\xi=\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{1/2}\xi\,\cdots\)(21) である。
また中心から\(\,\xi\,\)より内側に含まれる質量\(\,m(\xi)\,\)は
\(\quad m(\xi)=\displaystyle\int_0^{\alpha\xi}4\pi\rho r^2dr\)
\(\qquad\quad =4\pi\alpha^3\rho_c\displaystyle\int_0^{\xi}\xi^2\theta^Nd\xi\)
\(\qquad\quad =-4\pi\alpha^3\rho_c\displaystyle\int_0^{\xi}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)d\xi\)
\(\qquad\quad =-4\pi\alpha^3\rho_c\xi^2\dfrac{d\theta}{d\xi}\quad\cdots\,\)(22)
で与えられる。したがって星の全質量\(\,M\,\)は
\(\quad M=-4\pi\alpha^3\rho_c\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\)
\(\qquad =-4\pi\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{3/2}\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\,\cdots\,\)(23)
ここで星の平均密度\(\,\bar{\rho}\,\)を考える。
\(\quad\bar{\rho}=\dfrac{M}{4\pi R^3/3}=\dfrac{3M}{4\pi R^3}\)
\(\qquad =-\dfrac{3}{4\pi\alpha^3\xi_n^3}\cdot4\pi\alpha^3\rho_c\left(\xi^2\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\)
\(\qquad =-\dfrac{3}{\xi_n}\rho_c\left(\dfrac{d\theta}{d\xi}\right)_{\xi=\xi_n}\quad\cdots\,\)(24)
\(\bigstar\,\)太陽の場合
これらの関係を太陽の場合で計算してみる。太陽の半径\(\,R=6.960\times10^8\)m 、質量\(\,M=1.9884\times10^{30}\)kg と、恒星の場合に適応できるポリトロープ指数\(\,N=3\,\,(\Gamma=\frac{4}{3}\,)\,\)の計算結果を用いる。
\(\quad\rho_c=\dfrac{M}{4\pi R^3}\left(-\dfrac{1}{\xi}\dfrac{d\theta}{d/xi}\right)_{\xi=\xi_n}^{-1}\quad\cdots\,\)(25)
\(\qquad=\dfrac{1.9884\times10^{30}}{4\pi\cdot (6.960\times10^8)^3\cdot6.152\times10^{-3}}\)
\(\qquad =7.629\times10^4\,\mathrm{kg}\cdot\mathrm{m}^{-3}\)
理科年表によると\(\,\rho_c=1.56\times10^5\,\mathrm{kg}\cdot\mathrm{m}^{-3}\,\)*となっており、計算値は約半分である。まあ、天文のざっくりモデルでの計算としては、まあ合っているといえる。
*[参考]理科年表の太陽の構造の部分
 |
 |
 |
 |
半径\(\,r\,\)、ガス密度は半径の関数とし、\(\,\rho(r)\,\) と表すと、その半径より内側に含まれる質量\(\,M_r\,\)は
星は自分自身の重力による収縮しようとする力とガス圧による対抗する力がバランスして、力学平衡を保っている。