万有引力の法則と惑星の公転について

生きているうちに一度は納得しておきたい事柄はいろいろありますが、万有引力の法則はその一つでしょう。人類は力学法則を発見することでミサイルを飛ばしたりミサイルを迎撃したりすることができるようになったわけですが、その法則をちゃんと説明できるかと言われると心もとないところがあります。というわけで今回は惑星の公転についていろいろやっていきたいと思います。

ニュートンの時代には観測事実からケプラーの法則と呼ばれる次の3つが知られていました。

第一法則（楕円軌道の法則）: 惑星は太陽を一つの焦点とする楕円軌道上を動く
第二法則（面積速度一定の法則）: 惑星と太陽を結ぶ線が一定時間に描く面積は一定である
第三法則（調和の法則）: 惑星の公転周期の2乗は軌道長半径の3乗に比例する

ニュートンのプリンキピアという本では運動の3法則から出発してケプラーの法則を含む様々な現象を説明しています。力学の教科書に出てくるような微分方程式を使ったものとは異なり図形的な証明によって構成されています。

そこでこの記事では万有引力が距離の2乗に反比例することと、万有引力の逆2乗則を仮定したときに惑星が楕円軌道を描くことについて、プリンキピア流の図形的な説明と教科書的な微分方程式による説明を行って、死ぬ時の心残りを一つ減らしておきたいと思います。

運動の3法則

ニュートンによって体系化された古典的な力学は次の3つの法則をもとに組み立てられています。

第一法則（慣性の法則）: 物体は力が加えられない限り、静止または等速直線運動を続ける
第二法則（運動方程式）: 物体に力が加わると、運動量が力の方向に、力の大きさに比例して変化する( $\frac{dp}{dt}=F$ )
第三法則（作用・反作用の法則）: 2つの物体の間に相互に働く力は、大きさが等しく逆向きになる

この3法則から出発して古典力学の様々な法則が導かれていきます。

逆2乗則の証明（プリンキピア版）

まずはじめにプリンキピアでの幾何学的な証明についてまとめておきたいと思います。説明は

プリンキピアを読む―ニュートンはいかにして「万有引力」を証明したのか? (ブルーバックス)

作者: 和田純夫
出版社/メーカー: 講談社
発売日: 2009/05/21
メディア: 新書
クリック: 10回
この商品を含むブログ (6件) を見る

を参考にしながら自分が納得できるように今回のテーマに関わる部分を抜き出しました。詳細は本を参照してください

まずは証明に必要な補題をいくつか証明します。

面積速度一定の法則

まずは「公転する物体がある1点からの力に引かれて運動している場合、一定時間に物体と力の中心を結んだ線分が描く面積は一定である」ことを示します。この時点では引力が距離とどのような関係にあるかは何も仮定していません。引力が連続的に働いていると考えると難しいので、瞬間的な力が一定の時間間隔で断続的に働くものとし時間間隔がゼロになる極限を考えることにします。

力の中心をSとします。A点にあった物体はt秒後にB点に到達し、B点でS方向の引力を受けて方向転換し、またt秒後にはC点に到達したとします。引力を受けなかったらt秒後に到達していたはずの点をC'とします。引力の影響をBVとすれば、CはBC'とBVを平行四辺形で合成した点になります。
f:id:kivantium:20160207210205p:plain:w600
ここで、 $\triangle{SAB}$ と $\triangle{SBC'}$ はAB=BC'より底辺の長さが等しく高さが共通なので面積が等しいです。
また、 $\triangle{SBC'}$ と $\triangle{SBC}$ は底辺SBが共通で、SB // CC'より高さが等しいので面積が等しくなります。

従って、 $\triangle{SAB}$ と $\triangle{SBC}$ の面積が等しいので、tを0に近づけた極限を考えれば一定時間に物体と力の中心を結んだ線分が描く面積は一定であることが分かります。

等加速度運動における時間と移動距離の関係

加速度 $a(=F/m)$ で移動している物体が $t$ 秒間に移動する距離は $\frac{1}{2}at^2$ となります。これは積分やグラフを書くなどすれば分かります。

接線からのずれの長さ

曲線ABを考え、BからAにおける法線に下ろした垂線の足をC、Bから曲線のAにおける接線に下ろした垂線の足をDとする。このときDをAに近づけた極限で線分DBの長さは線分ABの長さの2乗に比例する。
f:id:kivantium:20160207214108p:plain:w600
Bを通りABに垂直な直線とACの交点をGとする。このとき $\triangle{ABD}\sim \triangle{GAB}$ なので $AB:BD=GA:AB$ 。よって $BD=\frac{AB^2}{GA}$ 。
素直な曲線であればBをAに近づけていくとGはある一点に近づいていくのでGAは固定されBDの長さはABの2乗に比例することが分かる。

ここで怪しいのは「Gはある一点に近づいていく」というところですが、曲線が円である場合には角ABGが垂直であることからGは常に円の中心になります。円でない場合にもよほど曲率が小さくない限りA付近でこの曲線に最も近い円の中心に近づいていくようです。

楕円の性質

S, Hを楕円の焦点とし、楕円上の点Pで接線を引く。また接線上のP以外の一点をR、接線を挟んでHの反対側にある点をH'とする。
f:id:kivantium:20160207224344p:plain:w600
このとき∠RPS = ∠HPZであることを示す。

Rは楕円の外にあることから
$SP + PH < SR + RH$
$PH = PH', \ \ RH = RH'$ より、
$SP + PH' < SR + RH'$
SからH'を結ぶ最短ルートはPを経由するものであることが分かる。
従って、S, P, H'は一直線上に並ぶ。
よって $\angle RPS = \angle H'PZ$
$\triangle{H'PZ}$ と $\triangle{HPZ}$ は合同なので $\angle H'PZ = \angle HPZ$
以上より $\angle RPS = \angle HPZ$

逆2乗則の証明

ここまでで準備が整ったので、「惑星の軌道が楕円であり向心力の中心が楕円の焦点である場合、向心力は距離の2乗に反比例する」ことの証明を行います。

物体が楕円軌道上のPからQに動いたとし、向心力の中心となる楕円の焦点をSとします。PRをPでの楕円の接線、DKは楕円の中心を通りPRに平行な弦、FはPからDKに下ろした垂線の足、TはQからPSに下ろした垂線の足、RはQRがPSと平行になるように取り、QXはRPと平行な直線とPSの交点、QVはQXの延長とPCの交点とします。
f:id:kivantium:20160211233841p:plain:w800

(1) PE=ACの証明

Hを楕円のもうひとつの焦点とすると、SC = HC。
点IをDK // IHとなるようにPS上に取ると、三角形の相似よりSE:EI = SC:CHとなるのでSE = EI。
$\begin{align} 2PE &= 2(PI + IE) \\ &= PI + (PI + IE + ES) \ \ \ (\because IE = ES) \\ &= PI + PS \\ &= PH + PS \\ \end{align}$

これは先述した楕円の性質より、∠RPI = ∠ZPH。
PR // HIより∠RPI = ∠PIH、∠ZPH = ∠PHI。
よって∠PIH = ∠PHIなので△PHIが二等辺三角形であることによります。

$\begin{align} PH + PS &= AS + AH \ \ \ (\because \textrm{楕円上の点は焦点からの距離の和が等しい}) \\ &= 2AC \ \ \ (\because \textrm{楕円は左右対称}) \\ \end{align}$
以上よりPE = ACが示されました。

(2) QR:PV = AC:PCの証明

XV // ECより△PEC∽△PXV。
よってPE : PC = PX : PV

□PXQRは平行四辺形なのでPX=QR。
よってQR : PV = PE : PC。

PE = ACを使って
PE : PC = AC : PC

以上よりQR : PV = AC : PC。

(3) $GV\cdot PV : QV^2 = PC^2 : CD ^2$ の証明

軌道が円の場合を考えると、CD = CPより、 $\frac{CP^2}{CD^2} = 1$
円の接線は、接点と中心を結ぶ直線と垂直なので
PR ⊥ PC。
QV // PRよりQV ⊥ PC。
よって△PQV ∽ △QGV
$\frac{PV}{QV} = \frac{QV}{VG}$ より、 $\frac{PV \cdot VG}{QV^2}=1$

楕円は円をCA方向に引き伸ばしたものであり、PV//VG//CP、QV//CDが保たれることを考えると分子分母は楕円でも同じ割合で変化するので、楕円の場合もこの式が成立します。
よって $GV\cdot PV : QV^2 = PC^2 : CD ^2$ 。

(4) (QがPに一致する極限で)QV:QT = CD:CBの証明

QがPに近づくとき、「接線からのずれの長さ」よりVXはQR（QPの2乗に比例して変化する）に比例して変化し、QXはQPの1乗に比例して変化するため、VXの変化はQXの変化に比べて急速です。
そのためQがPに近づく極限で、QV : QT = QX : QT。

QX // EFより∠QXT = ∠PEF。
∠PFE = ∠QTX = 90°
よって△QTX∽△PFEなのでQX : QT = PE : PF。
PE = ACよりPE : PF = AC : PF。

A, Bで楕円に接する長方形とP, D, C, Kで楕円に接する長方形は、円に接する同じ面積の正方形を同じ比率で拡大したものなので楕円でも面積が等しいままです。
よってCB x AC = CD x PF。
よってAC : PF = CD : CB。
以上より、QV:QT = CD:CB。

(5)まとめ

(2), (3), (4)の結果を左辺同士・右辺同士でそれぞれ掛けると、（(4)のみ2乗を掛ける）
$QR\cdot (GV\cdot PV) \cdot (QV)^2 : PV \cdot (QV)^2 \cdot (QT)^2 \\ = AC \cdot (PC)^2 \cdot (CD)^2 : PC\cdot (CD)^2\cdot (CB)^2$
が成立します。

整理して、
$QR\cdot GV : QT^2 = AC \cdot PC : CB^2$
が成立し、GV = 2PCを使うと
$\frac{QT^2}{QR} = \frac{2CB^2}{AC}$
となります。

ここで注目するべきことは、右辺がPにもQにも依存しない値になっていることです。
「等加速度運動における時間と移動距離の関係」より向心力を受けて落下した距離を表すQRは、向心力 x 時間^2に比例します。
よって、向心力はQR/(時間^2)に比例します。
「面積速度一定の法則」よりSP x QTが時間に比例するので、結局、
向心力は $\frac{QR}{(SP\cdot QT)^2}$ に比例することになります。
$\frac{QR}{QT^2}$ が定数であることを考えると、結局向心力はSPの2乗に反比例することが分かります。

SPというのは向心力の中心から惑星までの距離のことなので、これが示したかったことになります。

逆2乗則の証明（微分を使う場合）

ニュートンは微積分学の創始者として知られていますが、プリンキピアの中ではこれまで見てきた通り微分を明示的には使っていません。一方現在大学で教わる一般的な古典力学では微分をバンバン使っていろんなことを証明しています。比較として微分を使った証明も書いておきます。

相対運動の運動方程式

先ほどの証明では引力の中心が固定されていることを前提に話が進められていましたが、実際には引力の中心である太陽も動いています。この場合にも同じような運動方程式が成立することを確認しておきます。

質量が $m_1, m_2$ の物体の位置ベクトルを $\mathbf{r_1}, \mathbf{r_2}$ とし、それぞれが $f_{21}, f_{12}$ の力で引かれているとすると運動方程式は
f:id:kivantium:20160211214807p:plain:w400
$m_1 \ddot{\mathbf{r_1}} = f_{21} \\ m_2 \ddot{\mathbf{r_2}} = f_{12}　= -f_{21}$
となります。

これを整理すると
$\ddot{\mathbf{r_1}} - \ddot{\mathbf{r_2}} = (\frac{1}{m_1}+\frac{1}{m_2})f_21$
より、 $\mu = \frac{m_1m_2}{m_1+m_2}$ とおいて
$\mu \ddot{\mathbf{r_1}} - \ddot{\mathbf{r_2}} = f_{21}$
という運動方程式が成立します。

これは、片方を固定した状況を考えても、質量を少し変えれば相対運動についても同様の運動方程式を考えてよいことを示しています。
プリンキピアでも同様の説明があるようですがここでは省略します。

極座標の運動方程式

xy座標で回転運動を考えるのは難しいので極座標に変換して考えます。
証明は二次元極座標における運動方程式とその導出 | 高校数学の美しい物語を見ていただくとして、極座標での運動方程式は
$m(\ddot{r}-r\dot{\theta}^2)=F_r \\ m\dfrac{1}{r}\dfrac{d}{dt}(r^2\dot{\theta})=F_{\theta}$
のようになります。 $F_r$ は中心に向かう力、 $F_{\theta}$ はそれに垂直な力を表しています。

逆2乗則の証明

http://www.th.phys.titech.ac.jp/~muto/lectures/Gmech08/chap07.pdfを参考にしました。

中心に向けた力しか働かないという仮定より、 $\frac{d}{dt}(r^2\dot{\theta})=0$ が成り立つため面積速度は一定です。
定数 $h$ を用いて $r^2\dot{\theta}=h$ と置けます。
こうすると $F_r = m(\ddot{r}-\frac{h^2}{r^3})$ が成立します。

楕円軌道を描くという仮定より、 $r = \frac{l}{1+\varepsilon \cos \theta}$ と書ける（参考：二次曲線（楕円，放物線，双曲線）の極座標表示 | 高校数学の美しい物語）ので、 $1+\varepsilon \cos \theta = \frac{l}{r}$ 。
時間微分して、 $-\frac{l}{r^2}\dot{r}=-\varepsilon\dot{\theta}\sin \theta$ より、 $\dot{r} = \frac{h}{l}\varepsilon\sin\theta$ 。
これをさらに時間微分すると、
$\begin{align} \ddot{r} &= \frac{h}{l}\varepsilon\dot{\theta}\cos\theta \\ &=\frac{h}{l}\varepsilon\frac{h}{r^2}\cos\theta \\ &= \frac{h^2}{lr^2}(\frac{l}{r}-1) \\ &= \frac{h^2}{r^3}-\frac{h^2}{lr^2} \end{align}$
を得ます。これを運動方程式に代入して
$\begin{align} F_r &= m\left((\frac{h}{r^3}-\frac{h^2}{lr^2})-\frac{h^2}{r^3}\right) \\ &= -\frac{mh^2}{lr^2} \end{align}$
を得ます。これは引力が $r^2$ に反比例することを表しているのでこれで証明は完了です。

逆2乗則を仮定した時に惑星の軌道が円錐曲線になることの証明

プリンキピア版

概略だけ書くと、

軌道上の位置や速度などの初期状態を決めるとその条件を満たす円錐曲線が一つ書ける。（書く方法が具体的に紹介されているがここでは省略）
同じ条件を満たす軌道は一つしかないので惑星の軌道は円錐曲線になる

というような感じです。先ほどの逆2乗則の証明と比べるとあっさりしていました。

微分を使う場合

地球の公転軌道が楕円であることの導出 | 高校数学の美しい物語の説明が必要十分だったのでリンクを紹介しておしまいにします。

ルンゲクッタ法を使ったシミュレーション

ここまで全部数学っぽい感じで惑星の公転と万有引力の法則が得られる経緯を説明してきましたが、このブログの説明文にある「プログラムを使っていろいろやります」を達成するためにプログラムで惑星の軌道をシミュレーションしてみようと思います。

万有引力の法則が与えられると、惑星の軌道は微分方程式を解くことで計算できます。微分方程式を数値的に解く上で強力な方法の一つがルンゲクッタ法（ルンゲ＝クッタ法 - Wikipedia）です。
菊池誠さんによるとシンプレクティック・インテグレーターの方がいいらしいのですが、まずはオーソドックスなルンゲクッタを使ってみることにしました。

ルンゲクッタ法は $\frac{dy}{dt} = f(t, y), \ \ y(t_0)=y_0$ という微分方程式に対して、小さい数 $h$ を用いて
$k_1 = f(t_n, y_n) \\ k_2 = f(t_n+\frac{h}{2}, y_n+\frac{h}{2}k_1) \\ k_3 = f(t_n+\frac{h}{2}, y_n+\frac{h}{2}k_2) \\ k_4 = f(t_n+h, y_n+hk_3)$
を計算し、 $y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$ を順番に求めることで解を求めます。
惑星の運動の数値計算 (ルンゲ・クッタ法)を参考に
${\displaystyle \frac{d}{dt} \left( \begin{array}{c} x \\ y \\ v_x \\ v_y \end{array} \right) = \left( \begin{array}{c} v_x \\ v_y \\ x(x^2+y^2)^\frac{3}{2} \\ yx(x^2+y^2)^\frac{3}{2} \end{array} \right) }$
というベクトルに対する微分方程式に上のルンゲクッタ法を適用することでシミュレーションを行ってみます。

ソースコードは以下の通りです。

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;
	
// 物体の位置・速度を表すクラス
class Pos {
public:
    double x, y, vx, vy;
    Pos(void){};
    Pos(double i, double j, double k, double l) :x(i), y(j), vx(k), vy(l){};

    // Pos同士を+演算子で足せるようにする
    Pos operator +(const Pos &rvalue) {
        Pos tmp;
        tmp.x = x + rvalue.x;
        tmp.y = y + rvalue.y;
        tmp.vx = vx + rvalue.vx;
        tmp.vy = vy + rvalue.vy;
        return tmp;
    }
    // doubleの値を*演算子で掛けられるようにする
    Pos operator *(const double &rvalue) {
        Pos tmp;
        tmp.x = x * rvalue;
        tmp.y = y * rvalue;
        tmp.vx = vx * rvalue;
        tmp.vy = vy * rvalue;
        return tmp;
    }

    // ルンゲクッタのf()を計算する
    Pos f(void) {
        Pos tmp;
        double r = pow(x*x+y*y, -1.5);
        tmp.x = vx;
        tmp.y = vy;
        tmp.vx = -x*r;
        tmp.vy = -y*r;
        return tmp;
    }

};

int main() {
    // 初期値 x=1, y=0, vx=0, vy=1.2
    Pos w(1.0, 0.0, 0.0, 1.2);
    Pos p, k1, k2, k3, k4;

    // 小さめの刻み
    double h = 1e-2;

    // 10万回繰り返す
    for(int i=0; i<100000; ++i) {
        printf("%lf, %lf\n", w.x, w.y);
        // ルンゲクッタ
        p = w;
        k1 = p.f();
        p = w + k1 * (h*0.5);
        k2 = p.f();
        p = w + k2 * (h*0.5);
        k3 = p.f();
        p = w + k3 * h;
        k4 = p.f();
        
        w = w + (k1+k2*2.0+k3*2.0+k4)*(h/6.0);
    }
}

演算子オーバーロードの書き方がよく分かっていないので有識者からの指摘をお待ちしています。

./a.out > cor.txt

のように出力を保存して次のPythonコードでアニメーションを作りました。

matplotlibによるサイクロイドのアニメーションにて - matsulibの日記とDrawing and Animating Shapes with Matplotlib — Nick Charltonを参考にしました。

#! /usr/bin/env python
# -*- coding: utf-8 -*-

import csv
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation

# ファイル読み込み
f = open("cor.txt", "r")
reader = csv.reader(f)

cor_x = []
cor_y = []
for row in reader:
    cor_x.append(float(row[0]))
    cor_y.append(float(row[1]))
f.close()

# プロットするものの準備
fig = plt.figure()
ax = plt.axes(xlim=(-3, 3), ylim=(-3, 3))
earth = plt.Circle((1, 0), 0.1, fc='b')
sun = plt.Circle((0, 0), 0.2, fc='r')
time_text = ax.text(0.05, 0.9, 't=0', transform=ax.transAxes)
plt.gca().set_aspect('equal', adjustable='box')

# 初期化
def init():
    earth.center = (cor_x[0], cor_y[0])
    ax.add_patch(earth)
    ax.add_patch(sun)
    return earth,

# アニメーション描画
def animate(i):
    # 50回に1回だけ表示
    earth.center = (cor_x[i*50], cor_y[i*50])
    time_text.set_text("t={}".format(i*50))
    return earth,

# 軌道の表示
plt.plot(cor_x, cor_y, 'k')
# アニメーション
anim = animation.FuncAnimation(fig, animate, 
                               init_func=init, 
                               frames=200, 
                               interval=50,
                               blit=False,
                               repeat=False)
# アニメーションを保存する場合は次のようにする
#anim.save("earth.mp4")
plt.show()