月別アーカイブ: 2019年10月

レーン=エムデン方程式(3)

レーン=エムデン方程式は以下の形で表される

\(\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\,\)である。

[1]中心近傍での扱い
この方程式は\(\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]数値計算
この方程式は(その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)。これより第3の推定点\((\,x_i+\Delta x,y_i+k_3\,)\)を得る。第4の推定点まで求めたところで、これら4つの増分推定を\(\,1:2:2:1\,\)の重みを付けて平均化し、これを最終的な増分推定として\(\,y_{i+1}\,\)を求める(式13)。この式は Taylor展開の4次の精度を持つ。

\(\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)
この常微分方程式を連立1階微分方程式に変換する。

\(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}\,\,\)とした。

[3]数値計算による結果
\(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}\,\)*となっており、計算値は約半分である。まあ、天文のざっくりモデルでの計算としては、まあ合っているといえる。
*[参考]理科年表の太陽の構造の部分







レーン=エムデン方程式(2)


その1で導出したレーン=エムデン方程式は
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\quad\cdots\,\)(1)
で、\(\,\theta\,\)は無次元化した密度、\(\,\xi\,\)は無次元化した半径である。中心での境界条件は、\(\xi\to 0\,\)のとき\(\,\theta=1\,,\,\dfrac{d\theta}{d\xi}=0\,\)である。

[1]\(N=0\,\)の時の解析解
式(1)を一般の\(\,N\,\)について解析解は得られないが、\(N=0\,,\,1\,,\,5\,\)の時は得ることが出来る。
なお、以降は見やすくするために\(\,\xi\,,\,\theta\,\)を\(\,x\,,\,y\,\)と記述する。式(1)は
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-x^N\quad\cdots\,\)(2)
となる。ここに\(\,N=0\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-1\quad\cdots\,\)(3)
両辺に\(\,x^2\,\)をかけると
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-x^2\quad\cdots\,\)(4)

この式の両辺を\(\,x\,\)で積分すると
\(\quad x^2\dfrac{dy}{dx}=-\dfrac{1}{3}x^3+C_1\,\cdots\,\)(5)
となる。\(C_1\,\)は積分定数。両辺を\(\,x^2\,\)で割ると
\(\quad\dfrac{dy}{dx}=-\dfrac{1}{3}x+\dfrac{C_1}{x^2}\quad\cdots\,\)(6)

再度\(\,x\,\)で積分すると
\(\quad y=-\dfrac{1}{6}x^2-\dfrac{C_1}{x}+C_2\quad\cdots\,\)(7)
ここで、\(\,x=0\,\)のとき\(\,y=1\,\)の条件を満たすためには、\(C_1=0\,\)で\(\,C_2=1\,\)とする必要がある。よって、\(N=0\,\)の時は
\(\quad y=1-\dfrac{1}{6}x^2\quad\cdots\,\)(8)
となる。

[2]\(N=1\,\)の時の解析解
式(2)に\(\,N=1\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)+y=0\quad\cdots\,\)(9)
ここで、\(z=xy\,\)とすると
\(\quad\dfrac{dy}{dx}=\dfrac{d}{dx}\left(\dfrac{z}{x}\right)=\dfrac{1}{x}\dfrac{dz}{dx}-\dfrac{z}{x^2}\quad\cdots\,\)(10)

式(10)の両辺に\(\,x^2\,\)をかけると
\(\quad x^2\dfrac{dy}{dx}=x\dfrac{dz}{dx}-z\quad\cdots\,\)(11)
となり、これを\(\,x\,\)で微分すると
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=\dfrac{d}{dx}\left(x\dfrac{dz}{dx}-z\right)=x\dfrac{d^2z}{dx^2}\quad\cdots\,\)(12)
式(12)を式(9)に代入すると
\(\quad\dfrac{1}{x}\dfrac{d^2z}{dx^2}+y=\dfrac{d^2z}{dx^2}+z=0\quad\cdots\,\)(13)
という2階微分方程式になる。

この方程式の解は
\(\quad z=A\cos x+B\sin x\)
で、\(z=xy\,\)と変数を元に戻すと
\(\quad y=\dfrac{A\cos x}{x}+\dfrac{B\sin x}{x}\quad\cdots\,\)(14)

ここで、\(x\to 0\,\)の時、\(y\to 1\,\)なので、\(A=0\,,B=1\,\)となり
\(\quad y=\dfrac{\sin x}{x}\quad\cdots\,\)(15)
となる。

[3]\(N=5\,\)の時の解析解
式(2)に\(\,N=5\,\)を代入すると
\(\quad\dfrac{1}{x^2}\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-y^5\quad\cdots\,\)(16)

ここで、特殊な変数変換を行う。
\(\quad x=e^{-t}\quad\cdots\,\)(17)
\(\quad y=\sqrt{\dfrac{1}{2}e^t}\,Z\quad\cdots\,\)(18)
式(17)より
\(\quad\dfrac{dx}{dt}=-x\quad\cdots\,\)(17)

また、\(Z\,\)の\(\,x\,\)微分は
\(\quad\dfrac{dZ}{dx}=\dfrac{dZ}{dt}\dfrac{dt}{dx}=-\dfrac{1}{x}\dfrac{dZ}{dt}\quad\cdots\,\)(18)

微分演算子には以下の関係がある
\(\quad\dfrac{d}{dx}=-\dfrac{1}{x}\dfrac{d}{dt}\quad\cdots\,\)(19)

まず、式(18)を\(\,x\,\)で微分すると
\(\quad\dfrac{dy}{dx}=-\dfrac{1}{2\sqrt{2}}x^{-\frac{3}{2}}Z+\dfrac{1}{\sqrt{2}}x^{-\frac{1}{2}}\dfrac{dZ}{dx}\)
\(\qquad=-\dfrac{1}{2\sqrt{2}}x^{-\frac{3}{2}}Z-\dfrac{1}{\sqrt{2}}x^{-\frac{3}{2}}\dfrac{dZ}{dt}\quad\cdots\,\)(20)

式(20)の両辺に\(\,x^2\,\)をかけると
\(\quad x^2\dfrac{dy}{dx}=-\dfrac{1}{2\sqrt{2}}x^{\frac{1}{2}}Z-\dfrac{1}{\sqrt{2}}x^{\frac{1}{2}}\dfrac{dZ}{dt}\quad\cdots\,\)(21)

この式(21)を\(\,x\,\)で微分する
\(\quad\dfrac{d}{dx}\left(x^2\dfrac{dy}{dx}\right)=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z-\dfrac{1}{2\sqrt{2}}x^{\frac{1}{2}}\dfrac{dZ}{dt}\)
\(\qquad\qquad-\dfrac{1}{2\sqrt{2}}x^{-\frac{1}{2}}\dfrac{dZ}{dt}-\dfrac{1}{\sqrt{2}}x^{\frac{1}{2}}\dfrac{d}{dx}\left(\dfrac{dZ}{dt}\right)\)
\(\qquad=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z+\dfrac{1}{\sqrt{2}}x^{-\frac{1}{2}}\dfrac{d^2Z}{dt^2}\quad\cdots\,\)(22)

ここで
\(\quad -x^2y^5=-x^2\dfrac{1}{4\sqrt{2}}x^{-\frac{5}{2}}Z^5=-\dfrac{1}{4\sqrt{2}}x^{-\frac{1}{2}}Z^5\quad\cdots\,\)(23)
なので式(22)(23)を(16)に代入すると
\(\quad-\dfrac{1}{4}Z+\dfrac{d^2Z}{dt}=-\dfrac{1}{4}Z^5\quad\)より
\(\quad\dfrac{d^2Z}{dt^2}=\dfrac{1}{4}Z(1-Z^4)\quad\cdots\,\)(24)
という形に書き直すことができる。式(24)の両辺に\(\dfrac{dZ}{dt}\,\)をかけると
\(\quad\left(\dfrac{dZ}{dt}\right)\left(\dfrac{d^2Z}{dt^2}\right)=\dfrac{1}{4}Z\left(\dfrac{dZ}{dt}\right)-\dfrac{1}{4}Z^5\left(\dfrac{dZ}{dt}\right)\,\cdots\,\)(25)

ここで左辺は以下のように変形でき
\(\quad\left(\dfrac{dZ}{dt}\right)\left(\dfrac{d^2Z}{dt^2}\right)=\dfrac{1}{2}\dfrac{d}{dt}\left[\left(\dfrac{dZ}{dt}\right)^2\right]\quad\cdots\,\)(26)

右辺の計算も
\(\quad Z\left(\dfrac{dZ}{dt}\right)=\dfrac{1}{2}\dfrac{d}{dt}\left(Z^2\right)\quad\cdots\,\)(27)
と変形できるので、式(25)は
\(\quad\dfrac{d}{dt}\left[\left(\dfrac{dZ}{dt}\right)^2\right]=\dfrac{1}{4}\dfrac{d}{dt}\left(Z^2\right)-\dfrac{1}{12}\dfrac{d}{dt}\left(Z^6\right)\quad\cdots\,\)(28)

\(\,t\,\)で積分すると
\(\quad\left(\dfrac{dZ}{dt}\right)^2=\dfrac{Z^2}{4}-\dfrac{Z^6}{12}+D\quad\cdots\,\)(29)
となる。\(\,D\,\)は積分定数である。\(x\to 0\,\)のとき、\(Z=0\,dZ/dt=0\,\)なので、\(D=0\,\)である。よって次の式が得られる。

\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=\pm\dfrac{1}{2}dt\quad\cdots\,\)(30)
左辺の正負の符号は境界条件より負を選択する。
\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=-\dfrac{1}{2}dt\quad\cdots\,\)(31)

この方程式を解くために、さらなる変数変換を行う。
\(\quad\dfrac{1}{3}Z^4=\sin^2\phi\quad\cdots\,\)(32)
この式の微分形は以下のようになる。
\(\quad\dfrac{4}{3}Z^3\,dZ=2\sin\phi\,\cos\phi\,d\phi\quad\cdots\,\)(33)

式(32)(33)を(31)の左辺に代入すると
\(\quad\dfrac{dZ}{Z\sqrt{1-\frac{1}{3}Z^4}}=\dfrac{1}{Z\sqrt{1-\sin^2\phi}}\dfrac{3}{4Z^3}2\sin\phi\,\cos\phi\,d\phi\)
\(\qquad=\dfrac{3}{2}\dfrac{1}{Z^4}\dfrac{\sin\phi\,\cos\phi\,d\phi}{\sqrt{1-\sin^2\phi}}=\dfrac{1}{2\sin^2\phi}\dfrac{\sin\phi\,\cos\phi}{\cos\phi}d\phi\)
\(\qquad=\dfrac{1}{2}\dfrac{1}{\sin\phi}\,d\phi=\dfrac{1}{2}\mathrm{cosec}\,\phi\,d\phi\quad\cdots\,\)(34)
より(31)の方程式は以下の形に変形できる。

\(\quad\mathrm{cosec}\,\phi=-dt\quad\cdots\,\)(35)
両辺を積分すると
\(\quad\ln\left|\tan\dfrac{\phi}{2}\right|=-t+C’\quad\cdots\,\)(36)
より
\(\quad\tan\dfrac{\phi}{2}=Ce^{-t}\quad\cdots\,\)(37)

ここで、式(32)から
\(\quad Z^4=3\sin^2\phi=\dfrac{12\tan^2\frac{\phi}{2}}{\left(1+\tan^2\frac{\phi}{2}\right)^2}\quad\cdots\,\)(38)
この式(38)に式(37)を代入すると
\(\quad Z^4=\dfrac{12C^2e^{-2t}}{(1+C^2e^{-2t})^2}\)
\(\quad Z=\pm\left[\dfrac{12C^2e^{-2t}}{(1+C^2e^{-2t})^2}\right]^{\frac{1}{4}}\quad\cdots\,\)(39)

ここで、\(Z\,\)は式(18)より正で、式(17)(18)を代入すると
\(\quad y=\left[\dfrac{3C^2}{(1+C^2e^{-2t})^2}\right]^{\frac{1}{4}}=\left[\dfrac{3C^2}{(1+C^2x^2)^2}\right]^{\frac{1}{4}}\quad\cdots\,\)(40)

\(x=0\,\)で\(\,y=0\,\)という条件なので、
\(\quad 1=(3C^2)^{\frac{1}{4}}\quad\)より\(\quad C^2=\dfrac{1}{3}\quad\)が得られ
\(\quad y=\left[\dfrac{1}{(1+\frac{1}{3}x^2)^2}\right]^{\frac{1}{4}}=\dfrac{1}{\sqrt{1+\frac{1}{3}x^2}}\quad\cdots\,\)(41)
となる。







これまでの経緯

やってみたかった事

・高校では『天文研』に所属して、太陽の黒点の観測や夏の奥多摩での観測合宿などやってそれなりに、天文してました。完成はしませんでしたが、放課後交代で反射鏡を磨いていました。その後、大学と生業は『化学』という天文とは多少距離がある分野で長らく過ごしていました。
・やってみたかったと言うほどではないですが、いつかは...的な思いがあったのは、確かです。

放送大学で

・一応、定年(60歳)後を見据えて50代から、放送大学の専科履修生として勉強していました。定年後に放送大学の「自然と環境コース」に全科履修生として入学しました。(2013年度)
・それまで、専科履修生で取得した単位があったので、3年からの編入としました。昔の大学の単位を認定してくれる制度はあったようですが、面倒だと考えて語学も含めて全ての単位を放送大学で取得することにしました。その夏には卒業研究の申請を提出しています。
・それまで学習センターは試験や面接授業でしか利用していませんでしたが、いろいろなサークルなどがあり、先輩諸兄が指導してくれる環境がありがたかったです。卒論はどうせなら\(\,\mathrm{\TeX}\,\)を使いたいと思っていたところ、サークルにて\(\,\mathrm{\TeX}\,\)の講座を開催してくれたので、渡りに船と勉強させてもらいました。
・自然科学系の卒研に必要な科目と言ったら、まず数学です。一応、工学部だったのですが、昔取った杵柄が全く役に立ちませんでした。その数学を親切に教えてくださる先輩がいて、初歩から(再?)勉強させてもらっています。

天文での卒研

・卒業研究は天文を選択しました。報告書(卒論)はもちろん\(\,\mathrm{\TeX}\,\)で記述しました。放送大学のWAKABAの(資料)の「卒業研究資料_2014年度_06_自然と環境」に掲載されています。テーマはKeplar望遠鏡の観測データを使用した「連星の黒点挙動」です。

\(\,\mathrm{\TeX}\,\)について

・そんな感じで、備考録では\(\,\mathrm{\TeX}\,\)を使える環境をということで、ロリポップのレンタルサーバーの上で WordPress を使っています。
参考までに MathJaxの表示設定はこちらです。
<script type=”text/javascript” async=”” src=”https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_CHTML”></script>
この方式でなくともプラグインして[mathjax](括弧は半角)と文頭で宣言すればよいらしいですが...今のところ、うまく行かないので、文頭にスクリプトを書いています。

ポリトロープについて

・恒星は天文の勉強の基本です。ポリトロープは恒星の入門コースとなっています。これから導かれるレーン=エムデン方程式は、恒星物理学の基本です。約6年前勉強した内容を忘れないようにということで、ここに残します。
・振り子は最近までちゃんとした計算を知りませんでした。そんなんで、勉強した楕円積分とともに、備考録の最初に残しておきます。

レーン=エムデン方程式(1)

[1]Lane Emden 方程式とは

・宇宙物理学や流体力学において、レーン=エムデン方程式(Lane-Emden equation)は、球対称な密度分布を示す力学平衡にある自己重力流体を記述する微分方程式である。名称は宇宙物理学者のジョナサン・ホーマー・レーンとロバート・エムデンに由来する。
・Jonathan Homer Lane(1819-1880,ニューヨーク州のゲネシー生まれ)エール大学でオルムステッド(Denison Olmsted)に師事し、数学と自然科学を学ぶ。1857 年にワシントン市内に特許事務所を開き、特許関係の顧問や特許弁護士などの業務を始めた。1860 年に南北戦争が勃発すると、ペンシルバニア州北部の兄弟の住む村に移り、冷却用ガス膨張装置を用いた極低温装置の開発や天文学にも興味をもち、ひとりで文献に当たりながら太陽内部構造の研究に思いを深めていた。
・太陽の進化についてはジョン・ハーシェルやケルビン卿の収縮説に興味をもっていたが、内部構造を数学的に扱っていない点がレーンには不満であった。そこで彼は太陽をガス球と仮定し、重力平衡と熱力学過程を基本にして内部構造に関する基本方程式を書き始めた。
「太陽の理論的温度」(Lane,J.H.1870,Amer.J.Sci.&Arts,50,57-74) 自己重力のもとで球状に保たれるガス体の全体的平衡状態を取り扱った歴史的論文である。彼は太陽内部を完全気体と仮定したが、これは大胆な発想であった。完全気体を考えると\(\,p=R\rho T\)、静力学的平衡は\(\,dp=-\dfrac{Gm(r)}{r^2}\rho\,dr\,\)、内部は断熱的に対流平衡となるので、内部構造は\(\,1-\left(\dfrac{\rho}{\rho_c}\right)^{\gamma-1}=\displaystyle\int_0^x\dfrac{\mu(x)}{x^2}dx\,\)にて与えられる。

[2]アウグスト・リターのガス球論

・August Ritter(1826-1908)1853年にゲッチンゲン大学卒業。アーヘン工業大学(現RWTH-Aachen)の構造力学と一般力学の担当教授。「大気の高さに関する考察とガス状天体の内部構造について」ガス球の平衡・振動・安定性の問題に取り組み、自己重力にあるポリトロープガス球の平衡状態についての先駆的研究。
・ポリトロープとは圧力と密度に\(\,pV^n=K\,\) または\(\,p=K\rho^{1+1/n}\,\)の関係を満たす過程のこと。
\(\,T/T_c=y\,,\quad r/a=x\,\)とおいて、重力平衡式と組み合わせると以下の基本微分方程式が得られる。
\(\quad\dfrac{d^2y}{dx^2}+\dfrac{2}{x}\dfrac{dy}{dx}+y^n=0\,\quad\cdots\,\)(1)
この式は ポリトロープガス球の内部構造に関する基本方程式として、エムデンやエディントンに大きな影響を与えた。

[3]エムデンのガス球論

・ロベルト・エムデン(Jacob Robert Emden,1862-1940)スイスの天体物理学および気象学の専門家。1907 年にミュンヘン工科大学の教授。「ガス球論」出版。リターの基本方程式もそのままの形で再度導かれており、通常、エムデンの微分方程式と呼ばれている。リターらと異なり、エムデンはガス球半径を有限の場合と無限の場合を併行して考察している。エムデンの基本的考え方はレーンやリターと変わっていないが、理論が整備され、応用例が多いことからガス球理論の集大成と呼ばれ、4 版を重ねて、星の内部構造論の基本的テキストとなって広く知られるようになった。なお、エムデンの甥、マーチン・シュワルツシルドは 1958 に著書「星の構造と進化」を著わしているが、この書で歴史的なマイルストーンの研究者として、順にエムデン、エディントン、チャンドラセカールの 3 人の名前を挙げている。

[4]星の力学的安定

・連続の式:
 半径\(\,r\,\)、ガス密度は半径の関数とし、\(\,\rho(r)\,\) と表すと、その半径より内側に含まれる質量\(\,M_r\,\)は
\(\quad M_r=\displaystyle\int_0^r4\pi r’^2\rho(r)dr’\quad\cdots\,\)(2)
と書くことができる。
この式の両辺を\(\,r\,\)で微分して微分形式で書き表すと
\(\quad\dfrac{dM_r}{dr}=4\pi r^2\rho(r)\quad\cdots\,\)(3)
となる。この式は流体力学の連続の式(continuity equation)の一種である。

・静水圧平衡:
 星は自分自身の重力による収縮しようとする力とガス圧による対抗する力がバランスして、力学平衡を保っている。
\(\,r\,\)を中心からの距離、\(\,P\,\)をそこでの圧力、\(\,G\,\)を重力定数として
\(\quad\dfrac{dP}{dr}=-\dfrac{GM_r}{r^2}\rho\quad\cdots\,\)(4)
と表せる。
連続の式(3)、静水圧平衡の式(4)を解くことで、天体の内部構造を明らかにすることができる。 しかしこの段階では、独立変数が半径\(\,r\,\)で、従属変数は密度\(\,\rho\)、圧力\(\,P\)、質量\(\,M_r\,\)の3つで方程式は2つしかない。ためもうひとつ変数を結びつける式が必要となる。

・ポリトロープ関係
ここで、前出のポリトロープ関係と呼ばれる関係式を導入する。
\(\quad P=K\rho^{1+\frac{1}{N}}\quad\cdots\,\)(5)
\(\,K\,,\,\,N\,\)は定数のパラメータであり、\(\,N\,\)のことをポリトロープ指数と呼ぶ。

[5]レーン=エムデン方程式の導出

天体内部の構造を決めるためには、式(3)、(4)、(5)を解けば良いことになる。まず式(4)から
\(\quad\dfrac{r^2}{\rho}\dfrac{dP}{dr}=-GM_r\quad\cdots\,\)(6)
となり、両辺を\(\,r\,\)で微分すると
\(\quad\dfrac{d}{dr}\left(\dfrac{r^2}{\rho}\dfrac{dP}{dr}\right)=-G\dfrac{dM_r}{dr}\quad\cdots\,\)(7)
となる。これを式(3)に代入すると
\(\quad\dfrac{d}{dr}\left(\dfrac{r^2}{\rho}\dfrac{dP}{dr}\right)=-4\pi Gr^2\rho\quad\cdots\,\)(8)
となる。また、式(5)を\(\,r\,\)で微分すると
\(\quad\dfrac{dP}{dr}=K\left(1+\dfrac{1}{N}\right)\rho^{\frac{1}{N}}\dfrac{d\rho}{dr}\quad\cdots\,\)(9)
が得られる。式(9)を式(8)に代入すると
\(\quad K\left(1+\dfrac{1}{N}\right)\dfrac{d}{dr}\left(r^2\rho^{\frac{1}{N}-1}\dfrac{d\rho}{dr}\right)=-4\pi Gr^2\rho\quad\cdots\,\)(10)

ここで密度\(\,\rho(r)\,\)を星の中心密度\(\,\rho_c\,\)として
\(\quad\rho=\rho_c\theta^N\quad\cdots\,\)(11)
の形で表すこととする。\(\theta\,\)は無次元化した密度
両辺を\(\,r\,\)で微分すると
\(\quad\dfrac{d\rho}{dr}=\dfrac{d}{dr}\left(\rho_c\theta^N\right)=\rho_cN\theta^{N-1}\dfrac{d\theta}{dr}\quad\cdots\,\)(12)
また、式(11)から以下の関係が導かれる
\(\quad\rho^{\frac{1}{N}-1}=\rho_c^{\frac{1}{N}-1}\,\theta^{1-N}\quad\cdots\,\)(13)

式(10)に式(12)(13)を代入すると
\(\quad K\left(1+\dfrac{1}{N}\right)\dfrac{d}{dr}\left(r^2\rho_c^{\frac{1}{N}-1}\theta^{1-N}\rho_cN\theta^{N-1}\dfrac{d\theta}{dr}\right)\!=\!-4\pi Gr^2\rho_c\theta^N\quad\)
\(\quad\dfrac{K(N+1)\rho_c^{\frac{1}{N}-1}}{4\pi G}\dfrac{1}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\cdots\,\)(14)

また式(5)を星の中心にあてはめると
\(\quad P_c=K\rho_c^{1+\frac{1}{N}}\quad\cdots\,\)(15)
この式(15)を式(14)に代入して\(\,K\,\)を消去すると
\(\quad\dfrac{(N+1)P_c}{4\pi G\rho_c^2}\dfrac{1}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\cdots\,\)(16)

左辺の係数部分を以下のように定義すると
\(\quad\alpha\equiv\left[\dfrac{(N+1)K}{4\pi G}\rho_c^{\frac{1}{N}-1}\right]^{1/2}\quad\cdots\,\)(17)
\(\quad\dfrac{\alpha^2}{r^2}\dfrac{d}{dr}\left(r^2\dfrac{d\theta}{dr}\right)=-\theta^N\quad\quad\cdots\,\)(18)

ここで、半径を以下の形に無次元化する。
\(\quad\xi=\dfrac{r}{\alpha}\quad\cdots\,\)(19)
これより無次元化した密度\(\,\theta\,\)の\(\,r\,\)の微分は以下のようになる
\(\quad\dfrac{d\theta}{dr}=\dfrac{d\theta}{d\xi}\dfrac{d\xi}{dr}=\dfrac{1}{\alpha}\dfrac{d\theta}{d\xi}\quad\cdots\,\)(20)

微分演算子は以下の形になる
\(\quad\dfrac{d}{dr}=\dfrac{1}{\alpha}\dfrac{d}{d\xi}\quad\cdots\,\)(21)

式(18)にこれらを用いると
\(\quad\dfrac{1}{\xi^2}\dfrac{d}{d\xi}\left(\xi^2\dfrac{d\theta}{d\xi}\right)=-\theta^N\quad\cdots\,\)(22)
この式は無次元化した密度\(\,\theta\)、無次元化した半径\(\,\xi\,\)を結びつけるもので、
レーン=エムデン方程式と呼ばれている。