# 例題に見る2体問題について（高校生向け，執筆中）

間違っていたら指摘下さい．よろしくお願いします．by ITMZ．
今日中に書き上げたいんだがどうなるかな．2017/08/08
アクセス権限：Freely(明日になったら制限します．)　Lockedにしました．
単位，次元は今は省略する．文脈でわかるだろう．

あ～これ真面目に書いたらpdf文書で20ページぐらいになるやつや（）
なので今のところ雑に書く．（ごめん）

#### 作業履歴

2017/08/08
22:04休憩なう23:39再開なう2:05疲れたなう2:19今日中はむりポ
2017/08/09
10:49再開15:48休憩なう23:20再開1:58休憩
2017/08/10
14:24再開15:24Gunzi氏が来たので中断．1:44再開2:58休憩
2017/08/11
12:04再開15:19プログラム修正中17:45休憩20:13再開20:49これから講義なので中断
2017/08/13
O氏に進捗を見せたところ何箇所か間違いを指摘されたので修正すること．また理論パートと数値計算パートに分けたほうが良いとのこと．
2017/08/14

#### TODOリスト
- [ ] Todos
- [x] 間違い修正
- [ ] ケプラーの法則
- [ ] 周期を求める
- [ ] ~~RKF45概要~~埋込み型RKの概要
- [ ] タイムステップ幅制御詳細
- [x] ソースコード改良．Sample008~Sample010までは大幅な改良が必要．

## 数値計算例題
### 2体問題

$$
\left\{
\begin{array}{ll}
\dfrac{d^2 x}{dt^2} = - \dfrac{x}{(x^2 + y^2)^\frac{3}{2}} \\
\dfrac{d^2 y}{dt^2} = - \dfrac{y}{(x^2 + y^2)^\frac{3}{2}}
\end{array}
\right.
$$

ただし，

$$
\begin{align}
x(0) &= 1-e,\ \ \dfrac{dx}{dt}(0) = 0,\\
y(0) &= 0, \ \ \ \ \ \ \ \ \ \dfrac{dy}{dt}(0) = \sqrt{\dfrac{1+e}{1-e}}.
\end{align}
$$

例としてパラメータ$e$を$e = 0.9$として数値計算しなさい．（$e$が1に近いほどみょんみょん度が高く，数値的に解きづらいことが知られている．）

## 問題の背景

まず何故この問題が”2体”問題なのか．どこからこの方程式がきたのか？この初期値の意図は？ひとつひとつ紐解いていこう．

### 重心の運動

はじめに万有引力の法則について復習．距離が$r$離れた質点間の力$F$は以下の通りである．ただし，質点の位置ベクトルを$\vec{x}_{1},\vec{x}_{2}$，質量を$m_{1},m_{2}$とする．

$$
F = G \dfrac{m_{1} m_{2}}{r^2}
$$

ここで，$G$は万有引力定数である．この法則により2点のニュートンの運動方程式を立てると，

$$
m_{1} \dfrac{d^2 \vec{x_{1}}}{dt^2} = - \dfrac{G m_{1} m_{2}}{\mid \mid \vec{x}_{1} - \vec{x}_{2}\mid \mid ^2} \dfrac{\vec{x}_{1} - \vec{x}_{2}}{\mid \mid \vec{x}_{1} - \vec{x}_{2}\mid \mid} \\
m_{2} \dfrac{d^2 \vec{x_{2}}}{dt^2} = - \dfrac{G m_{1} m_{2}}{\mid \mid \vec{x}_{2} - \vec{x}_{1}\mid \mid ^2} \dfrac{\vec{x}_{2} - \vec{x}_{1}}{\mid \mid \vec{x}_{2} - \vec{x}_{1}\mid \mid}
$$

となる．（もうすこし丁寧な説明が必要か）

上記2式を足すと，

$$
m_{1} \dfrac{d^2 \vec{x_{1}}}{dt^2} + m_{2} \dfrac{d^2 \vec{x_{2}}}{dt^2} = \vec{0}
$$

少し式変形をすれば，

$$
\dfrac{1}{m_{1} + m_{2}} \left( m_{1} \dfrac{d^2 \vec{x_{1}}}{dt^2} + m_{2} \dfrac{d^2 \vec{x_{2}}}{dt^2} \right) = \vec{0}
$$

すなわち，~~「重心にはたらく力はゼロ」が結論される．~~これより，「重心は等速直線運動する」．

### 相対運動
2式を次のように変形する．
$$
\dfrac{d^2 \vec{x_{1}}}{dt^2} = - \dfrac{G m_{2}}{\mid \mid \vec{x}_{1} - \vec{x}_{2}\mid \mid ^2} \dfrac{\vec{x}_{1} - \vec{x}_{2}}{\mid \mid \vec{x}_{1} - \vec{x}_{2}\mid \mid} \\
\dfrac{d^2 \vec{x_{2}}}{dt^2} = - \dfrac{G m_{1}}{\mid \mid \vec{x}_{2} - \vec{x}_{1}\mid \mid ^2} \dfrac{\vec{x}_{2} - \vec{x}_{1}}{\mid \mid \vec{x}_{2} - \vec{x}_{1}\mid \mid}
$$
次に上の2式の差を考えることにより，相対的（質点1に対する質点2）な運動を考える．
ここで，
$$
\begin{align}
\vec{r} &= \vec{x_{1}} - \vec{x_{2}} \\
r &= \mid \mid \vec{x_{1}} - \vec{x_{2}} \mid \mid
\end{align}
$$
とおくと，
$$
\begin{align}
\dfrac{d^2 \vec{x_{1}}}{dt^2} &= - \dfrac{G m_{2}}{r^{3}} \vec{r} \\
\dfrac{d^2 \vec{x_{2}}}{dt^2} &= \dfrac{G m_{1}}{r^{3}} \vec{r}
\end{align}
$$
差を計算する．
$$
\dfrac{d^2 \vec{x_{1}}}{dt^2} - \dfrac{d^2 \vec{x_{2}}}{dt^2} = - \dfrac{G\left( m_{1} + m_{2} \right)}{r^{3}} \vec{r}
$$
ここで，
$$
\mu = G (m_{1} + m_{2})
$$
とおくと，
$$
\dfrac{d^2}{dt^2} \left( \vec{x_{1}} - \vec{x_{2}} \right) = - \mu \dfrac{\vec{r}}{r^{3}}
$$
$\vec{r} = \vec{x_{1}} - \vec{x_{2}}$なので
$$
\dfrac{d^2 \vec{r}}{dt^2} = - \mu \dfrac{\vec{r}}{r^{3}}
$$
となる．これがどちらかの質点から見た質点の運動方程式である．言ってみれば太陽から見た地球の動きを記述する方程式である．（そして，太陽と地球の重心は等速直線運動をしている．）この方程式の解が楕円軌道になればケプラーの法則の一つは数学という体系の演繹によって立証される．~~2体問題は相対運動を考えることで，1体問題にある意味帰着される．太陽から見れば地球という1体を考えるだけでいい．（本当はちょっと違うがだいたいこんなところだ．詳細は二体問題に関するWikipediaでも良いだろう．ここの説明眠くて曖昧なので明日書き直すこと．）~~
（$\mu$は換算質量でよく使うギリシャ文字らしい．違う記号に書き換えたほうが良いだろうか．）

### タイムスケール変換
時刻$t$を以下のようにスケール変換する．（本当は無次元化の話をしたいがまあ今はいいだろう．）
$$
\tau := \alpha t
$$
方程式
$$
\dfrac{d^2 \vec{r}}{dt^2} = - \mu \dfrac{\vec{r}}{r^{3}}
$$
は，
$$
\dfrac{1}{\alpha^{2}} \dfrac{d^2 \vec{r}}{d\tau^2} = - \mu \dfrac{\vec{r}}{r^{3}}
$$
変形すると，
$$
\dfrac{d^2 \vec{r}}{d\tau^2} = - \alpha^2 \mu \dfrac{\vec{r}}{r^{3}}
$$
ここで，
$$
\alpha = \sqrt{\dfrac{1}{\mu}}
$$
とおけば，
$$
\dfrac{d^2 \vec{r}}{d\tau^2} = - \dfrac{\vec{r}}{r^{3}}
$$
変数をおき直せば
$$
\dfrac{d^2 \vec{r}}{dt^2} = - \dfrac{\vec{r}}{r^{3}}
$$
となる．ここで，
$$
\vec{r} = \left[
\begin{array}{c}
x \\
y
\end{array}
\right]
$$
とおけば，
$$
r = \mid \mid \vec{r} \mid \mid = \sqrt{x^{2} + y^{2} }
$$
まとめると，微分方程式
$$
\dfrac{d^2 \vec{r}}{dt^2} = - \dfrac{\vec{r}}{r^{3}}
$$
は，以下のように
$$
\dfrac{d^2}{dt^2}
\left[
\begin{array}{c}
x \\
y
\end{array}
\right] = - \dfrac{1}{\sqrt{x^2 + y^2}^3}\left[
\begin{array}{c}
x \\
y
\end{array}
\right]
$$
となる．
ベクトルの成分ごとに書き下せば，
$$
\left\{
\begin{array}{ll}
\dfrac{d^2 x}{dt^2} = - \dfrac{x}{(x^2 + y^2)^\frac{3}{2}} \\
\dfrac{d^2 y}{dt^2} = - \dfrac{y}{(x^2 + y^2)^\frac{3}{2}}
\end{array}
\right.
$$
問題の方程式を得る．

#### （みーくんメモ）
$$
\alpha = \sqrt{\dfrac{1}{\mu}} = \sqrt{\dfrac{1}{G\left(m_{1} + m_{2}\right)}}
$$
太陽質量：　1.989 × 10^30 kg
地球の質量：5.972 × 10^24 kg

万有引力定数：6.67408 × 10^-11 m^3 kg^-1 s^-2

((1.98900 × (10^30)) + (5.97200 × (10^24))) * 6.67408 × (10^(-11)) =
1.3274785e+20

sqrt(1 / 1.3274785e+20) =
8.679331e-11

## 軌道の厳密解を求める

次の初期条件の謎について．なぜあの式なのか考えよう．

初等的な解き方，ちょっと高級な解き方．いろいろある．
まずは，$e$が離心率であることを確認したい．前提知識は高校数学＋（行列計算）までの範囲にしたい．
しかし，新たな知識も必要になるかもしれない．その時は説明する．

解くべき方程式（スケール変換する前）は，

$$
\dfrac{d^2}{dt^2}
\left[
\begin{array}{c}
x \\
y
\end{array}
\right] = - \dfrac{\mu}{\left(x^2 + y^2\right)^\frac{3}{2}}\left[
\begin{array}{c}
x \\
y
\end{array}
\right]
$$

であった．これをそのままの直交座標系で解くのは難しい．極座標で計算する．すなわち，

$$
x = r \cos \theta\\
y = r \sin \theta
$$

と変数変換する．（$r$は上の$r$とはもちろん違うもの．）正確に書けば，$r$も$\theta$も$t$の関数なので，

$$
x(t) = r(t) \cos \theta(t)\\
y(t) = r(t) \sin \theta(t)
$$

と書いたほうが正確だが，上のように略記することにする．

問題の方程式の左辺，

$$
\dfrac{d^2}{dt^2}
\left[
\begin{array}{c}
x \\
y
\end{array}
\right] = \dfrac{d^2}{dt^2}
\left[
\begin{array}{c}
r \cos \theta \\
r \sin \theta
\end{array}
\right]
$$

を計算する．まず，$t$の1階微分から計算する．

$$
\begin{align}
\dfrac{d}{dt}
\left[
\begin{array}{c}
r \cos \theta \\
r \sin \theta
\end{array}
\right] &=
\left[
\begin{array}{c}
\dfrac{dr}{dt} \cos \theta - r \sin \theta \dfrac{d \theta}{dt} \\
\dfrac{dr}{dt} \sin \theta + r \cos \theta \dfrac{d \theta}{dt}
\end{array}
\right]\\
& =
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
\dfrac{dr}{dt}\\
r \dfrac{d \theta}{dt}
\end{array}
\right]
\end{align}
$$

さらに$t$で微分する．（2階微分）

$$
\begin{align}
\dfrac{d^2}{dt^2}
\left[
\begin{array}{c}
r \cos \theta \\
r \sin \theta
\end{array}
\right]
& =
\dfrac{d}{dt}
\left(
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
\dfrac{dr}{dt}\\
r \dfrac{d \theta}{dt}
\end{array}
\right]
\right) \\
& =
\dfrac{d}{dt}
\left(
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right] \right)
\left[
\begin{array}{c}
\dfrac{dr}{dt}\\
r \dfrac{d \theta}{dt}
\end{array}
\right]
+
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\dfrac{d}{dt}
\left( \left[
\begin{array}{c}
\dfrac{dr}{dt}\\
r \dfrac{d \theta}{dt}
\end{array}
\right]
\right)
\end{align}
$$

1項目について計算すれば，（ここれへんの計算もっと（初等的に易しく）楽にならないか？）

$$
\begin{align}
\dfrac{d}{dt}
\left(
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right] \right)
\left[
\begin{array}{c}
\dfrac{dr}{dt}\\
r \dfrac{d \theta}{dt}
\end{array}
\right] &=
\left[
\begin{array}{rr}
-\sin \theta \dfrac{d \theta}{dt} & - \cos \theta \dfrac{d \theta}{dt} \\
\cos \theta \dfrac{d \theta}{dt} & - \sin \theta \dfrac{d \theta}{dt}
\end{array}
\right]
\left[
\begin{array}{c}
\dfrac{dr}{dt}\\
r \dfrac{d \theta}{dt}
\end{array}
\right] \\
&=
\left[
\begin{array}{rr}
-\sin \theta & - \cos \theta \\
\cos \theta & - \sin \theta
\end{array}
\right]
\left[
\begin{array}{c}
\dfrac{dr}{dt} \dfrac{d \theta}{dt}\\
r \left( \dfrac{d \theta}{dt} \right)^2
\end{array}
\right] \\
&=
\left[
\begin{array}{rr}
\sin \theta & \cos \theta \\
-\cos \theta & \sin \theta
\end{array}
\right]
\left[
\begin{array}{c}
-\dfrac{dr}{dt} \dfrac{d \theta}{dt}\\
-r \left( \dfrac{d \theta}{dt} \right)^2
\end{array}
\right] \\
&=
\left[
\begin{array}{rr}
\cos \theta & \sin \theta \\
\sin \theta & -\cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
-r \left( \dfrac{d \theta}{dt} \right)^2\\
-\dfrac{dr}{dt} \dfrac{d \theta}{dt}
\end{array}
\right] \\
&=
\left[
\begin{array}{rr}
\cos \theta & -\sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
-r \left( \dfrac{d \theta}{dt} \right)^2\\
\dfrac{dr}{dt} \dfrac{d \theta}{dt}
\end{array}
\right]
\end{align}
$$

2項目に関して計算すると，

$$
\begin{align}
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\dfrac{d}{dt}
\left( \left[
\begin{array}{c}
\dfrac{dr}{dt}\\
r \dfrac{d \theta}{dt}
\end{array}
\right]
\right) = \left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
\dfrac{d^2r}{dt^2}\\
\dfrac{dr}{dt} \dfrac{d \theta}{dt} + r \dfrac{d^2 \theta}{dt^2}
\end{array}
\right]
\end{align}
$$

1項目と2項目を足すと，

$$
\begin{align}
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
\dfrac{d^2r}{dt^2} -r \left( \dfrac{d \theta}{dt} \right)^2\\
2\dfrac{dr}{dt} \dfrac{d \theta}{dt} + r \dfrac{d^2 \theta}{dt^2}
\end{array}
\right]
\end{align}
$$

右辺について変数変換すると，

$$
\begin{align}
-\dfrac{\mu}{\left(x^2 + y^2\right)^\frac{3}{2}}\left[
\begin{array}{c}
x \\
y
\end{array}
\right] &= - \dfrac{\mu}{r^3}
\left[
\begin{array}{c}
r \cos \theta \\
r \sin \theta
\end{array}
\right] \\
&= - \dfrac{\mu}{r^2}
\left[
\begin{array}{c}
\cos \theta \\
\sin \theta
\end{array}
\right]
\end{align}
$$
以上より，
$$
\begin{align}
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
\dfrac{d^2r}{dt^2} -r \left( \dfrac{d \theta}{dt} \right)^2\\
2\dfrac{dr}{dt} \dfrac{d \theta}{dt} + r \dfrac{d^2 \theta}{dt^2}
\end{array}
\right] = - \dfrac{\mu}{r^2}
\left[
\begin{array}{c}
\cos \theta \\
\sin \theta
\end{array}
\right]
\end{align}
$$

両辺を左から

$$
\begin{align}
\left[
\begin{array}{rr}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta
\end{array}
\right]^{-1} =
\left[
\begin{array}{rr}
\cos \theta & \sin \theta \\
-\sin \theta & \cos \theta
\end{array}
\right]
\end{align}
$$

をかけると，

$$
\begin{align}
\left[
\begin{array}{c}
\dfrac{d^2r}{dt^2} -r \left( \dfrac{d \theta}{dt} \right)^2\\
2\dfrac{dr}{dt} \dfrac{d \theta}{dt} + r \dfrac{d^2 \theta}{dt^2}
\end{array}
\right] &= - \dfrac{\mu}{r^2}
\left[
\begin{array}{rr}
\cos \theta & \sin \theta \\
-\sin \theta & \cos \theta
\end{array}
\right]
\left[
\begin{array}{c}
\cos \theta \\
\sin \theta
\end{array}
\right]\\
&= - \dfrac{\mu}{r^2}
\left[
\begin{array}{c}
1 \\
0
\end{array}
\right]
\end{align}
$$

成分ごとに書くと，

$$
\begin{align}
\left\{
\begin{array}{ll}
\dfrac{d^2r}{dt^2} -r \left( \dfrac{d \theta}{dt} \right)^2 = - \dfrac{\mu}{r^2}\\
2\dfrac{dr}{dt} \dfrac{d \theta}{dt} + r \dfrac{d^2 \theta}{dt^2} = 0
\end{array}
\right.
\end{align}
$$

となる．2式目を変形すると，

$$
\begin{align}
2\dfrac{dr}{dt} \dfrac{d \theta}{dt} + r \dfrac{d^2 \theta}{dt^2} &= 0 \\
\dfrac{1}{r} \dfrac{d}{dt} \left( r^{2} \dfrac{d \theta}{dt} \right) &= 0 \\
\dfrac{d}{dt} \left( r^{2} \dfrac{d \theta}{dt} \right) &= 0
\end{align}
$$

よって，$C$を定数として次の式が得られる．（面積速度一定の法則）

$$
r^{2} \dfrac{d \theta}{dt} = C
$$

:::info
**面積速度一定の法則（ケプラーの第2法則）**

時刻$t$までに軌道と原点を繋ぐ直線（動径）が掃く面積を$S(t)$としよう．面積速度一定の法則とは「一定時間$T$の間に動径が掃く面積は一定である」というものである．定式化すると，

$$
\begin{align}
\forall t,\ \forall T\in \mathbb{R},\ \dfrac{S(t+T) - S(t)}{T} = \rm{const.}
\end{align}
$$

となる．

**証明**

時刻$t_{1}$から$\Delta t$秒間の$S(t)$の変化率を考える．

まず，$r_{min}, r_{max}$を次のように定義する．

$$
\begin{align}
r_{min} &:= \min \{ r(t) \mid t_{1} < t < t_{1} + \Delta t \} \\
r_{max} &:= \max \{ r(t) \mid t_{1} < t < t_{1} + \Delta t \}
\end{align}
$$

また，この$\Delta t$秒間で動径のなす角が$\Delta \theta$だけ変化したとすると，S(t)のこの間の変化は，

$$
\begin{align}
\dfrac{1}{2} r_{min}^{2} \Delta \theta \leq S(t_{1} + \Delta t) - S(t_{1}) \leq \dfrac{1}{2} r_{max}^{2} \Delta \theta
\end{align}
$$

両辺を$\Delta t$でわると，

$$
\begin{align}
\dfrac{1}{2} r_{min}^{2} \dfrac{\Delta \theta}{\Delta t} \leq \dfrac{S(t_{1} + \Delta t) - S(t_{1})}{\Delta t} \leq \dfrac{1}{2} r_{max}^{2} \dfrac{\Delta \theta}{\Delta t}
\end{align}
$$

となり，$\Delta t \to 0$を考えると，定義から，$r_{min} \to t_{1},\ t_{max} \to t_{1}$となり，

$$
\begin{align}
\dfrac{1}{2} r_{1}^{2} \dfrac{d \theta}{d t} \leq \dfrac{d S}{d t}(t_{1}) \leq \dfrac{1}{2} r_{1}^{2} \dfrac{d \theta}{dt}
\end{align}
$$

はさみうちの原理より，

$$
\begin{align}
\dfrac{d S}{d t}(t_{1}) = \dfrac{1}{2} r_{1}^{2} \dfrac{d \theta}{dt}
\end{align}
$$

これは任意の時刻で成立するので，

$$
\begin{align}
\dfrac{d S}{d t} = \dfrac{1}{2} r^{2} \dfrac{d \theta}{dt}
\end{align}
$$

（上の$\frac{dS}{dt}$を面積速度という．）ここで，先程の式を思い出すと，

$$
r^{2} \dfrac{d \theta}{dt} = C
$$

であった，ゆえに定数$D$を用いて，

$$
\begin{align}
\dfrac{d S}{d t} = D
\end{align}
$$

が成立する．よって$S(t)$は定数$F$を用いて，

$$
S(t) = Dt + F
$$

と書ける．

よって，任意の時刻$T$秒間の間に掃く面積は，

$$
\begin{align}
\dfrac{S(t+T) - S(t)}{T} = \dfrac{D(t+T) + F - Dt -F}{T} = \dfrac{DT}{T} = D = \rm{const.}
\end{align}
$$

となり，一定である．□
:::

:::success
実は一般に中心力で運動する系であれば，この法則は成り立つ．
:::

先の式を変形して，

$$
\dfrac{d \theta}{dt} = \dfrac{C}{r^2}
$$

これを先の1式目

$$
\dfrac{d^2r}{dt^2} -r \left( \dfrac{d \theta}{dt} \right)^2 = - \dfrac{\mu}{r^2}
$$

に代入して，

$$
\begin{align}
\dfrac{d^2r}{dt^2} -r \left( \dfrac{C}{r^2} \right)^2 = - \dfrac{\mu}{r^2}\\
\dfrac{d^2r}{dt^2} - \dfrac{C^2}{r^3} = - \dfrac{\mu}{r^2}
\end{align}
$$

$r$についての微分方程式を得た．

これからこの微分方程式を解く．いま，

$$
s := \dfrac{1}{r}
$$

とおくと，
$$
\begin{align}
\dfrac{d r}{dt} &= - \dfrac{1}{s^2} \dfrac{d s}{dt} \\
&= - \dfrac{1}{s^2} \dfrac{d s}{d \theta} \dfrac{d \theta}{dt} \\
&= - \dfrac{1}{s^2} s^{2} C \dfrac{d s}{d \theta} \\
&= - C \dfrac{d s}{d \theta}
\end{align}
$$

よって，

$$
\begin{align}
\dfrac{d^2 r}{dt^2} &= - C \dfrac{d^2 s}{d \theta^2} \dfrac{d \theta}{dt} \\
&= - s^{2} C^{2} \dfrac{d^2 s}{d \theta^2}
\end{align}
$$

これより，

$$
\begin{align}
\dfrac{d^2r}{dt^2} - s^{3} C^{2} = - \mu s^{2}
\end{align}
$$

に代入して，

$$
\begin{align}
-s^{2} C^{2} \dfrac{d^2 s}{d \theta^2} - s^{3} C^{2} &= - \mu s^{2} \\
\dfrac{d^2 s}{d \theta^2} &= -s + \dfrac{\mu}{C^{2}} \\
\dfrac{d^2}{d \theta^2}\left( s - \dfrac{\mu}{C^{2}} \right) &= - \left(s - \dfrac{\mu}{C^{2}} \right)
\end{align}
$$

を得る．これを解くと，

$$
\begin{align}
s - \dfrac{\mu}{C^{2}} = A \cos (\theta + \alpha)
\end{align}
$$

$s$を$\dfrac{1}{r}$に戻すと，

$$
\dfrac{1}{r} = \dfrac{\mu}{C^{2}} + A \cos (\theta + \alpha)
$$

よって，

$$
\begin{align}
r &= \dfrac{1}{\dfrac{\mu}{C^{2}} + A \cos (\theta + \alpha)} \\
r &= \dfrac{\dfrac{C^{2}}{\mu}}{1 + \dfrac{A C^{2}}{\mu} \cos (\theta + \alpha)}
\end{align}
$$

まとめると，

$$
\begin{align}
r(\theta) = \dfrac{a}{1 + b \cos \theta}
\end{align}
$$

ただし，（軌道を調べるためには$\alpha = 0$として良い．）
$$
a = \dfrac{C^{2}}{\mu},\ \ \ b = \dfrac{A C^{2}}{\mu}
$$

これが楕円軌道になりえることを確かめる．

$$
\begin{align}
r(1 + b \cos \theta ) &= a \\
r + b r \cos \theta &= a \\
r + b x &= a \\
r &= a - bx \\
r^{2} &= (a - bx)^{2} \\
x^{2} + y^{2} &= (a - bx)^{2} \\
(1 - b^{2}) \left( x^{2} + \dfrac{2ab}{1 - b^{2}}x \right) + y^{2} &= a^{2} \\
(1 - b^{2}) \left( x + \dfrac{ab}{1 - b^{2}} \right)^{2} + y^{2} &= a^{2} + \dfrac{a^{2}b^{2}}{1-b^{2}} \\
\left( x + \dfrac{ab}{1 - b^{2}} \right)^{2} + \dfrac{y^{2}}{1 - b^{2}} &= \dfrac{a^{2}}{1-b^{2}} + \dfrac{a^{2}b^{2}}{(1-b^{2})^2}
\end{align}
$$

よって，楕円軌道になりうる．（上の計算に自信が無い．(´・ω・｀)　後でしっかり確認する．）

楕円軌道になるための条件は，

$$
1-b^2>0.
$$

だろうか？

楕円の標準形は，
$$
\dfrac{x^2}{a^2} + \dfrac{y^2}{b^2} = 1.
$$
極座標では
$$
r = \dfrac{l}{1 + e \cos \theta}
$$
である．ここで$e$は離心率である．

### この問題を初期条件について考える

なぜ問題の$e$が離心率になるのだろう．
考えてみる．

$$
r^{2} \dfrac{d \theta}{dt} = C
$$

（ただし$C$は定数．）

が成立していた．$C$を初期条件から求めよう．

初期条件を振り返ってみよう．

$$
\begin{align}
x(0) &= 1-e,\ \ \dfrac{dx}{dt}(0) = 0,\\
y(0) &= 0, \ \ \ \ \ \ \ \ \ \dfrac{dy}{dt}(0) = \sqrt{\dfrac{1+e}{1-e}}.
\end{align}
$$

$x(0)=1-e,\ \dfrac{d y}{dt}(0)=\sqrt{\dfrac{1+e}{1-e}}$と以上の条件より，

$$
\begin{align}
r^{2} \dfrac{d \theta}{dt} &= (1-e)^{2} \sqrt{\dfrac{1+e}{1-e}} \\
&=\sqrt{(1-e)(1+e)}
\end{align}
$$

なので，

$$
C = \sqrt{(1-e)(1+e)}
$$

となる．

これにより，未定だった定数，

$$
\dfrac{1}{r} = \dfrac{\mu}{C^{2}} + A \cos \theta
$$

$A$が求まる．（問題より$\mu=1$．）$\theta = 0$のとき初期条件を用いて，

$$
\begin{align}
\dfrac{1}{1-e} &= \dfrac{1}{(1-e)(1+e)} + A \\
A &= \dfrac{1}{1-e} - \dfrac{1}{(1-e)(1+e)} \\
A &= \dfrac{e}{(1-e)(1+e)}
\end{align}
$$

よって，以下の求めた軌道の方程式

$$
\begin{align}
r = \dfrac{C^{2}}{1 + A C^{2} \cos \theta}
\end{align}
$$

に対して，

$$
C = \sqrt{(1-e)(1+e)}
$$

$$
A = \dfrac{e}{(1-e)(1+e)}
$$

を代入すると

$$
r = \dfrac{(1-e)(1+e)}{1 + e \cos \theta}
$$

よって，$e$は離心率である．

わかりやすく直交座標で書き直す．

$$
\begin{align}
r &= \dfrac{(1-e)(1+e)}{1 + e \cos \theta}\\
r + e r\cos \theta &= (1-e)(1+e) \\
r + e x &= 1 - e^{2} \\
r^2 &=(1-e^{2}-ex)^2 \\
x^{2} + y^{2} &= (1-e^{2})^{2} - 2e(1-e^{2})x + e^{2}x^{2} \\
(1 - e^{2})x^{2} + 2e(1-e^{2})x + y^{2} &= (1-e^{2})^{2} \\
(1-e^{2})(x^{2} + 2ex) + y^{2} &= (1-e^{2})^2 \\
(1-e^{2})(x + e)^{2} - e^{2}(1-e^{2}) + y^{2} &= (1-e^{2})^{2} \\
(1-e^{2})(x + e)^{2} + y^{2} &= (1-e^{2})^{2} + e^{2}(1-e^{2}) \\
(1-e^{2})(x + e)^{2} + y^{2} &= (1-e^{2})
\end{align}
$$

よって，~~(ここ計算ミスあるな．)~~

$$
\begin{align}
(x + e)^{2} + \dfrac{y^{2}}{1-e^{2}} = 1
\end{align}
$$

となり，

$$
\begin{align}
0\leq e < 1
\end{align}
$$

のとき，楕円である．よって，数値計算例題では楕円軌道をえがく．

## 例題の楕円軌道の周期を求める

面積速度は，

$$
\dfrac{dS}{dt}= \dfrac{1}{2} r^{2} \dfrac{d \theta}{dt}
$$

であり，

$$
\dfrac{1}{2} r^{2} \dfrac{d \theta}{dt} = \dfrac{C}{2} = \dfrac{1}{2} \sqrt{(1-e)(1+e)}
$$

であったので，

$$
\dfrac{dS}{dt} = \dfrac{1}{2} \sqrt{(1-e)(1+e)}
$$

また，楕円

$$
\begin{align}
(x + e)^{2} + \dfrac{y^{2}}{1-e^{2}} = 1
\end{align}
$$

の面積$A$は，

$$
A = \pi \cdot 1 \cdot \sqrt{1-e^2} = \pi \sqrt{(1-e)(1+e)}
$$

よって，周期$T_p$は

$$
T_{p} = \dfrac{A}{\dfrac{dS}{dt}} = \dfrac{\pi \sqrt{(1-e)(1+e)}}{\dfrac{1}{2} \sqrt{(1-e)(1+e)}} = 2 \pi
$$

である．（偶然？定数だった！！！！）

:::success
**ケプラーの第3法則**
惑星の公転周期$T_{p}$の2乗は，軌道の長半径$a$の3乗に比例する．すなわち，

$$
\dfrac{T^{2}}{a^{3}} = \rm{const.}
$$

##### 数値計算例題の場合

自明．

#### 次のような一般の場合

書きかけ．

:::

## 数値計算法
問題再掲．
$$
\left\{
\begin{array}{ll}
\dfrac{d^2 x}{dt^2} = - \dfrac{x}{(x^2 + y^2)^\frac{3}{2}} \\
\dfrac{d^2 y}{dt^2} = - \dfrac{y}{(x^2 + y^2)^\frac{3}{2}}
\end{array}
\right.
$$
ただし，初期条件として，
$$
\begin{align}
x(0) &= 1-e,\ \ \dfrac{dx}{dt}(0) = 0,\\
y(0) &= 0, \ \ \ \ \ \ \ \ \ \dfrac{dy}{dt}(0) = \sqrt{\dfrac{1+e}{1-e}}.
\end{align}
$$
を課す．
これをコンピュータで解く．

解き方としてRK4やRKF45で解くことを目指す．

解き方は今までやった方法と同様．2階の微分方程式を1解の微分方程式に落とす．正確に言うと，2変数2階連立常微分方程式を4変数1階連立常微分方程式に変形する．

$$
\dfrac{dx}{dt} = z,\\
\dfrac{dy}{dt} = w.
$$
とでも置いてやれば良い．そうすることで，
$$
\left\{
\begin{array}{ll}
\dfrac{d x}{dt} = z \\
\dfrac{d z}{dt} = - \dfrac{x}{(x^2 + y^2)^\frac{3}{2}} \\
\dfrac{d y}{dt} = w\\
\dfrac{d w}{dt} = - \dfrac{y}{(x^2 + y^2)^\frac{3}{2}}
\end{array}
\right.
$$
と1階の連立常微分方程式に変形できる．ここで見やすさのために変数名を$x_{1},x_{2},y_{1},y_{2}$を用いて置き換えれば，
$$
\left\{
\begin{array}{ll}
\dfrac{d x_{1}}{dt} = x_{2} \\
\dfrac{d x_{2}}{dt} = - \dfrac{x_{1}}{(x_{1}^2 + y_{1}^2)^\frac{3}{2}} \\
\dfrac{d y_{1}}{dt} = y_{2}\\
\dfrac{d y_{2}}{dt} = - \dfrac{y_{1}}{(x_{1}^2 + y_{1}^2)^\frac{3}{2}}
\end{array}
\right.
$$
を解くことと同等．ここで，
$$
\vec{x} := \left[
\begin{array}{c}
x_1 \\
x_2 \\
y_1 \\
y_2
\end{array}
\right],
$$
として，このベクトル$\vec{x}$からベクトルを返す関数として次のように$\vec{f}(\vec{x})$を定義して，
$$
\vec{f}(\vec{x}) := \left[
\begin{array}{c}
x_2 \\
-\dfrac{x_{1}}{(x_{1}^2 + y_{1}^2)^\frac{3}{2}} \\
y_2 \\
-\dfrac{x_{1}}{(x_{1}^2 + y_{1}^2)^\frac{3}{2}}
\end{array},
\right]
$$
とすれば，方程式は，
$$
\dfrac{d \vec{x}}{dt} = \vec{f}(\vec{x})
$$
と簡潔に書くことができる．$\vec{f}$のC++風擬似コードとしては以下の通り，

//R^4からR^4への関数f．
vector4dim f(vector4dim u){
    Float x       = u(0);
    Float x_dot   = u(1);
    Float y       = u(2);
    Float y_dot   = u(3);

    Float R = x * x + y * y;

    return vector4dim{
        x_dot,
        - x / mp::sqrt(R * R * R),
        y_dot,
        - y / mp::sqrt(R * R * R)
    };
}

ここまで書ければ，あとは例えばオイラー法であれば，以前教えたように（時間があれば詳細を書くこと）離散化すれば
$$
\dfrac{\vec{x}_{n+1} - \vec{x}_{n}}{\Delta t} = \vec{f}(\vec{x}_{n})
$$
のようにあまりベクトルに慣れていなくても形式的にこうしてもらって構わないし，これで正しいメソッドである．見慣れた形に変形すると，
$$
\vec{x}_{n+1} = \vec{x}_{n} + \Delta t \vec{f}(\vec{x}_{n})
$$
これに初期条件を加えて計算すれば良い．

古典的ルンゲ・クッタ法についても同様である．

（RKF45等については上の章を書いたら．）
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
みょ～ん
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
a
　／￣＼
○ ／￣￣￣＼ヘ
　 /・　 ・　 ＼>
／￣￣￣＼　　 Ｖ|
｜ ――― ｜　　｜|
＼＿＿＿／　　 ∧|
　　＼　　　　／〉
　　　￣￣￣￣￣
　　ヽ(´･ω･)ﾉ
　　　 |　 /
　　　 UU
（わからない）かつ（質問しないorググらない）なら、アンコウ投げんぞ!