0. はじめに
筆者は最近まで古代中国の暦法についての記事を書いていた。その中で驚いたのは，古代の中国人達は現代のレベルで見ても驚くべきほどの精度で天体観測をしており，それに基づく暦法を運用していたことである。
- 儀鳳暦のマジックナンバー1340の謎を解く
- 儀鳳暦のマジックナンバー1340の謎を解く(2)天体物理の基礎用語
- 儀鳳暦のマジックナンバー1340の謎を解く(3)李淳風は調日法を本当に使用したのか？
- 儀鳳暦のマジックナンバー1340の謎を解く(4)共通分母って本当に使い易い？
中国では唐代の末期（822～891）の71年間しかわれなかった宣命暦であるが，日本では貞観四年（862）から貞享元年（1684）まで実に823年間も用いられた。中国では受命改制の思想により王朝交代の度に暦法を変えていたが，日本では政権の実体は武家側（幕府）にあっても朝廷自体は存続しており，朝廷側に暦法の権限が残されていたからだ。もう一つの理由として，日本でもようやく平和な江戸時代になって天文学が成熟し，独自の暦法を制定できるだけの実力がついたということだろう。
欧州ではヨハネス・ケプラー（1571～1630）がケプラーの第三法則を記した「宇宙の調和」を1619年に出版し，アイザック・ニュートン（1642～1727）がいわゆるニュートン力学を用いて天体の運動を説明したプリンキピア（自然哲学の数学的諸原理）を1687年に出版した。一方，日本では麻田剛立（1734～1799）がそれらとは独自にケプラーの第三法則を発見したとされる。ちなみに麻田の弟子の高橋至時（1764～1804）のさらに弟子がかの有名な伊能忠敬（1745～1818）である。
ということで，自分もケプラーの第三法則を導いてみよう。江戸時代の人には負けられない。
ちなみに惑星軌道を円軌道に近似してよいのなら，高校物理でも第三法則を直接導くことができるが，真面目に楕円軌道で計算しようすると周到な準備が必要となる。このため初歩的なレベルからケプラーの法則全てを導いていくが，使える武器は筆者の高校レベルの数学である。
1. ケプラーの法則
念のため，ケプラーの法則を以下に示す。
- 第一法則：惑星の軌道は楕円であること。また太陽は焦点の位置に存在すること
- 第二法則：惑星と太陽とを結ぶ線分が単位時間に描く面積は一定になること
- 第三法則：惑星の公転周期の二乗は軌道長半径の三乗に比例すること
2. 楕円の式を極座標に変換する
まず直交座標における楕円の式を示す。
\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \tag{1}
$a \ge b$ のとき離心率 $\varepsilon = \sqrt{1 - (b/a)^2}$ とおくと楕円の二つ焦点 $F,F'$ の座標は $(\pm\varepsilon\hspace1px a,0)$ となる。このとき楕円の軌道上の点 $P$ から二つの焦点 $F,F'$ までの距離の和は一定となる。
r + r' = 2a \tag{2}
楕円の円周上の点 $P$ の座標を $(\varepsilon\hspace1px a + r\cos \theta, r\sin \theta)$ とおくと余弦定理より
(r')^2 = r^2 + (2\varepsilon a)^2 - 2 r (2\varepsilon a) \cos(\pi - \theta) \tag{3}
となる。式 $(2),(3)$ より $r'$ を消去して $r$ について求めると
r = \frac{b^2}{a(1 + \varepsilon \cos \theta)} \tag{4}
となる。
一般的に離心率（eccentricity）は $e$ という記号で表すことが多いが，本記事ではネイピア数 $e$ も使用しているため，区別のため離心率は $\varepsilon$ で表す。
3. 極座標系における運動方程式
極座標系におけるベクトル ${\bf r}$ を単位ベクトル ${\bf e}^{j\hspace1px\theta}$ を用いて次のように表す。
\def\sp{\hspace{1px}}
{\bf r} = r\,{\bf e}^{j\sp\theta} \tag{5}
時間微分 $\dot{\bf r}$ は次のようになる。※$\dot{\bf e}^{j\hspace1px\theta} = j\hspace2px\dot{\theta}\hspace2px{\bf e}^{j\hspace1px\theta}$
\def\sp{\hspace{1px}}
\dot{\bf r} = \dot{r}\,{\bf e}^{{j\sp\theta}} + j\,\dot{\theta}\,r\,{\bf e}^{j\sp\theta} \tag{6}
二階時間微分 $\ddot{\bf r}$ は次のようになる。
\def\sp{\hspace{1px}}
\begin{aligned}
\ddot{\bf r} &= \ddot{r}\,{\bf e}^{j\sp\theta}
+ j\,\dot{\theta}\,\dot{r}\,{\bf e}^{j\sp\theta}
+ j\,\ddot{\theta}\,r\,{\bf e}^{j\sp\theta}
+ j\,\dot{\theta}\,\dot{r}\,{\bf e}^{j\sp\theta}
+ (j\,\dot{\theta})^2\,r\,{\bf e}^{j\sp\theta} \\
&= (\ddot{r} - \dot{\theta}^2\,r)\,{\bf e}^{j\sp\theta}
+ j\,(\ddot{\theta}\,r + 2\dot{\theta}\,\dot{r})\,{\bf e}^{j\sp\theta} \\
\end{aligned} \tag{7}
ベクトル表記について，高校では一般的に矢印記法が用いられると思うが，本記事では業界標準的な太字記法を採用した。微分記法についても同じくニュートン記法を採用した。本当は黒板太字を採用したかったが，思っていたほど見易いわけではないこと，ドットの位置が離れ過ぎていることから断念した。
| ベクトル記法 | 微分記法 | ||
|---|---|---|---|
| ニュートン | ラグランジュ | ライプニッツ | |
|
矢印 |
$\dot{\vec{e}}$ | $\vec{e}'$ | $\displaystyle \frac{d\hspace1px\vec{e}}{dt}$ |
|
太字 |
$\dot{\bf e}$ | ${\bf e'}$ | $\displaystyle \frac{d\hspace1px\bf e}{dt}$ |
|
黒板太字 |
$\dot{\mathbb{e}}$ | $\mathbb{e'}$ | $\displaystyle \frac{d\hspace1px\mathbb{e}}{dt}$ |
4. 微分方程式を立てる
惑星の質量 $m$，太陽の質量 $M$，重力定数 $G$ とし，太陽の位置を原点におく。惑星に働く力のうち ${\bf e}^{j\hspace{0.1em}\theta}$ 方向に働く力を $f_r$ とし，$j\hspace{0.1em}{\bf e}^{j\hspace{0.1em}\theta}$ 方向に働く力を $f_\theta$ とおくと，${\bf e}^{j\hspace{0.1em}\theta}$ 方向にのみ太陽による引力が働くので
\begin{align}
f_r &= m (\ddot{r} - \dot{\theta}^2\,r) = -\frac{GMm}{r^2} \tag{8} \\
f_\theta &= m(\ddot{\theta}\,r + 2\dot{\theta}\,\dot{r}) = 0 \tag{9}
\end{align}
となる。これより惑星の質量 $m$ とは関係なく
\begin{align}
\ddot{r} - \dot{\theta}^2\,r +\frac{GM}{r^2} &= 0 \tag{10} \\
\ddot{\theta}\,r + 2\dot{\theta}\,\dot{r} &= 0 \tag{11}
\end{align}
という連立微分方程式を解けばよい。
5. ケプラーの第二法則
先にケプラーの第二法則を導く。惑星と太陽とを結ぶ線分が単位時間に描く面積 $\Delta S$ とおくと
\Delta S = \frac{1}{2} \dot{\theta}\,r^2 \tag{12}
となる。面積 $\Delta S$ の時間微分を考えると式 $(11)$ を用いて
\Delta\dot{S} = \frac{1}{2} \ddot{\theta}\,r^2 + \dot{\theta}\,r\,\dot{r}
= \frac{1}{2} r \left( \ddot{\theta}\,r + 2\dot{\theta}\,\dot{r} \right)
= 0
\tag{13}
となる。よって $\Delta S$ は一定（時間不変）となる。以後，面積速度 $\Delta S$ を定数として用いる。
6. ケプラーの第一法則
先に目的を述べると式 $(4)$ を導くことが目標であるが，$q = 1/r$ とおいて
q = \frac{a}{b^2} (1 + \varepsilon\,\cos\theta) \tag{14}
を導くほうが簡単そうである。さらに式 $(12)$ を $\theta$ で二階微分すると
\frac{d^2q}{d\theta^2} = -\frac{a}{b^2}\,\varepsilon\,\cos\theta \tag{15}
となるので式 $(14)$ と式 $(15)$ の和より
\bbox[4pt, border:2px dotted red]{\frac{d^2q}{d\theta^2} + q = \frac{a}{b^2}} \tag{16}
となる。この式 $(16)$ への到達を目標と改めよう。
まず，前項の結果（ケプラーの第二法則）である式 $(12)$ より
\dot{\theta} = \frac{2\Delta S}{r^2} \tag{17}
となる。これを式 $(10)$ に代入すると
\ddot{r} - \frac{4S^2}{r^3} + \frac{GM}{r^2} = 0 \tag{18}
となる。ここで $q = 1/r$ とおくと式 $(17)$ に代入して
\dot{\theta} = 2\Delta S q^2 \tag{19}
となる。また $r$ の時間微分を $q$ の時間微分に置き換え，さらに $dq/d\theta$ に置き換える。式 $(19)$ の関係を用いて $\dot{\theta}$ を消去すると
\dot{r} = \frac{d}{dt}\left(\frac{1}{q}\right) = - \frac{\dot{q}}{q^2} = - \frac{1}{q^2} \frac{dq}{d\theta} \dot{\theta} = -2\Delta S \frac{dq}{d\theta} \tag{20}
となり，さらに $r$ の二階時間微分も同様にして
\ddot{r} = -2\Delta S \frac{d}{dt}\left(\frac{dq}{d\theta}\right) = -2\Delta S \frac{d^2q}{d\theta^2} \dot{\theta} = -4\Delta S^2 q^2\frac{d^2q}{d\theta^2} \tag{21}
となる。これを式 $(18)$ に代入して $r$ を $q$ で置き換えると
-4\Delta S^2 q^2\frac{d^2q}{d\theta^2} - 4\Delta S^2 q^3 + GMq^2 = 0 \tag{22}
\bbox[4pt, border:2px dotted blue]{\frac{d^2q}{d\theta^2} + q = \frac{GM}{4\Delta S^2}} \tag{23}
となり，目標の微分方程式を得られた。
7. ケプラーの第三法則
7.1 円軌道に近似した場合
太陽系の惑星の軌道は表２に示すように極めて円に近い。円軌道に近似してしまうとケプラーの第三法則をとても簡単に導くことができる。
式 $(10)$ より円軌道の場合は $\ddot{r} = 0$ となるから
m\,r\,\dot{\theta}^2 = \frac{GMm}{r^2} \tag{24}
となる。これは遠心力と重力が釣り合っているという意味になる。これより $\dot{\theta}$ について求めると
\dot{\theta} = \sqrt{\frac{GM}{r^3}} \tag{25}
となる。円軌道の場合 $r = a = b$ となり，角速度 $\dot{\theta}$ も一定となるから公転周期 $T$ とおくと
\frac{2\pi}{T} = \dot{\theta} = \sqrt{\frac{GM}{a^3}} \tag{26}
T^2 = \frac{4\pi^2}{GM} a^3 \tag{27}
となり，公転周期の二乗は軌道長半径の三乗に比例することを示した。
| 惑星 | 離心率 $\varepsilon$ | 長短径比 $a/b$ |
|---|---|---|
| 水星 | 0.2056 | 1.021830 |
| 金星 | 0.0068 | 1.000023 |
| 地球 | 0.0167 | 1.000139 |
| 火星 | 0.0934 | 1.004391 |
| 木星 | 0.0485 | 1.001178 |
| 土星 | 0.0555 | 1.001544 |
| 天王星 | 0.0463 | 1.001074 |
| 海王星 | 0.0090 | 1.000041 |
| （冥王星） | 0.2502 | 1.032864 |
7.2 楕円軌道で計算した場合
式 $(16)$ と式 $(23)$ との係数比較より
\frac{a}{b^2} = \frac{GM}{4\Delta S^2} \tag{28}
という関係が得られる。また公転周期 $T$ とおくと楕円の面積 $\pi\hspace1px a\hspace1px b$ であるから面積速度 $\Delta S$ は
\Delta S = \frac{\pi\,a\,b}{T} \tag{29}
となる。これを式 $(28)$ に代入すると
\frac{a}{b^2} = \frac{GM}{4} \left(\frac{T}{\pi\,a\,b}\right)^2 = \frac{GMT^2}{4\pi^2\,a^2\,b^2} \tag{30}
T^2 = \frac{4\pi^2}{GM} a^3 \tag{31}
となり，式 $(27)$ と同じ結果を得られた。
8. まとめ
以上より，高校レベルの数学でも無事ケプラーの法則三つ全て導くことができた。実は十年ほど前にも一度チャレンジしたことがあったが，その当時は実力不足で断念せざるを得なかったのだ。いきなり極座標形式での運動方程式を出されても高校生（レベルの数学能力しか持たない者）には意味不明だよね。以下参考文献を示すが，一番役に立った（お世話になった）のはヨビノリたくみ先生の動画[6]かもしれない。
9. 参考文献
[1] 嘉数次人，江戸時代の天文学【4】江戸幕府の天文学その3，天文教育2007年11月号，日本天文教育普及研究会 
[2] D.L.グッドスティーン，ファインマンの特別講義 惑星運動を語る，岩波書店，2017年
[3] 楕円の極座標表示－オタクof数理の共同ブログ
[4] ケプラーの法則を導出する｜高校生から味わう理論物理入門
[5] 大学物理のフットノート｜力学｜ケプラーの第3法則
[6] 予備校のノリで学ぶ大学の数学・物理【大学物理】力学入門⑤(極座標における運動)【力学】
[7] 山﨑一磨，TeXの入門教材を作りました。ご自由にお使いください。- Qiita