6.2惑星の運動

惑星の運動を考える。導入する値は以下の通り。
質量：M
位置ベクトル：\bf{r}=\bf{r}(t)
万有引力定数：\gamma
太陽の質量：M_s
\bf{r}方向の単位ベクトル：\bf{e}_r
こられを用いると、万有引力の法則より、惑星の運動方程式は
M\frac{d^2\bf{r}}{dt^2}=-\frac{\gamma M_s M}{r^2}\bf{e}_r
となる。k=\gamma M_sと書き、極座標(r,\hspace{3}\theta)を持つ平面を考えると、\bf{r}=r\bf{e}_rであるから、上式は
\frac{d^2\bf{r}}{dt^2}=-\frac{k}{r^2}\bf{e}_r
と書き改められる。
　
ドットを時間微分とすると、
\frac{d\bf{r}}{dt}=\dot r \bf{e}_r + r \dot \theta \bf{e}_{\theta}
\frac{d^2\bf{r}}{dt^2}=(\ddot r- r\dot {\theta^2})\bf{e}_r + (2\dot r \dot \theta + r\ddot \theta)\bf{e}_{\theta}
なので、
(\ddot r- r(\dot \theta)^2)\bf{e}_r + (2\dot r \dot \theta + r\ddot \theta)\bf{e}_{\theta}=-\frac{k}{r^2}\bf{e}_r
が得られる。係数比較により、
\ddot r- r(\dot \theta)^2=-\frac{k}{r^2}…(1)
2\dot r \dot \theta + r\ddot \theta=0…(2)
となる。さて、(2)は
\frac{1}{r}\frac{d}{dt}(r^2 \dot\theta)=0
であるから、積分して
r^2\dot \theta=h
が得られる。hは定数である。\dot\thetaを(1)に代入すると、
\ddot r-\frac{h^2}{r^3}=-\frac{k}{r^2}…(3)
という、非線型常微分方程式が得られる。
p=\frac{1}{r}という変換を用いると、
\begin{align}\frac{dp}{d\theta}&=\frac{dp}{dt}\frac{dt}{d\theta}&=\frac{d}{dt}(\frac{1}{r})\frac{1}{\dot\theta}&=-\frac{\dot r}{r^2\dot\theta}&=-\frac{\dot r}{h}\end{align}
\frac{d^2p}{d\theta ^2}=-\frac{r^2\ddot r}{h^2}
となるので、(3)に戻ると、
\frac{d^2p}{d\theta^2}+p=\frac{k}{h^2}
と、線型第2階微分方程式に帰着できる。
補助解はp=A\cos\theta + B\sin\theta、ひとつの特解はp=\frac{k}{h^2}である。B=0となるようにすると、
r=\frac{1}{p}=\frac{\frac{h^2}{k}}{1+\frac{Ah^2}{k}\cos\theta}
となる。
　
さて、極座標表示される楕円
r=\frac{l}{1+e\cos\theta}
は、離心率e=\frac{Ah^2}{k}、半通径l=\frac{h^2}{k}で形づけられるらしい。そういう運動らしい。
　
たいていの惑星にとって、周回軌道はほぼ円に近似できるらしいので、r=\frac{h^2}{k}の円軌道と仮定する。\dot \theta=\frac{h}{r^2}だったので、
\frac{d\theta}{dt}=\frac{h}{r^2}=\frac{k^{\frac{1}{2}}}{r^{\frac{3}{2}}}
を、全軌道で積分すると、惑星の時間周期Tは
2\pi=(\frac{k^{\frac{1}{2}}}{r^{\frac{3}{2}}})T
T=(\frac{2\pi}{k^{\frac{1}{2}}})r^{\frac{3}{2}}
\log T =K + \frac{3}{2}\log r
と書ける。
本中の値を用いると

earth <- 149.5*10^6 #太陽と地球の距離
distance <- c(0.387, 0.723, 1, 1.524, 5.204, 9.5, 19.2, 30.1, 39.5) * earth #各惑星の太陽間との距離
T.time <- c(0.24, 0.61, 1, 1.91, 11.9, 29.5, 84, 165, 248) #公転周期
plot(distance, T.time, log="xy")
abline(glm(log(T.time, 10) ~ log(distance, 10))$coefficients)
glm(log(T.time) ~ log(distance)) #回帰

Call:  glm(formula = log(T.time) ~ log(distance))

Coefficients:
  (Intercept)  log(distance)  
       -28.24           1.50  

Degrees of Freedom: 8 Total (i.e. Null);  7 Residual
Null Deviance:	    53.54 
Residual Deviance: 0.0003495 	AIC: -59.87

と、傾きが\frac{3}{2}になっている!!