オイラー法を用いた天体の運動シミュレーション

# はじめに
太陽の周りをまわる地球のような1点の天体の動きは、近似的に中心力を用いた運動方程式である微分方程式で近似的に叙述することができる。この場合、力学的エネルギーと角運動量は保存される。今回は、オイラー法を用いた数値解析でこのことを示したい。具体的には、定点からの中心力を与えることで、以下のような2次元空間における天体の動きをシミュレーションする。

![x_y.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/eeb82d3a-5f5f-4977-e84a-618637b77d60.png)

ただし、今回の課題は、以下のリンクにある『数値計算入門』の課題5.3を参考にした。

https://opac.ll.chiba-u.jp/da/curator/900121236/compute.pdf

# 問題設定
定点$0$から距離$r$だけ離れた物体に働く中心力を$F(r)=-r^\alpha$とすると、$xy$平面における微分方程式は以下のように表される。ただし質量$m=1$とする。

```math
\frac{d^2 x}{d t^2}= -r^\alpha \frac{x}{r},\frac{d^2 y}{d t^2}= -r^\alpha \frac{y}{r},r=\sqrt{x^2+y^2}
```

この場合において、角運動量は、速度と変位を用いて以下の式で与えられる。

```math
L=x v_y-y v_x
```

一方で、ポテンシャルエネルギーは以下のように与えられる。

```math
V(r)=-\int_{r=\infty}^{r=r}F(r)dr=-\int_{r=\infty}^{r=r}-r^\alpha dr=

\begin{equation}
\left\{ \,
    \begin{aligned}
    & \frac{r^\alpha}{\alpha+1} (\alpha \ne -1) \\
    & log(r) (\alpha=-1) \\
    \end{aligned}
\right.
\end{equation}

```

一方で、微小区間の物体のポテンシャルエネルギーを分割して考えると以下のようになる。

```math
\Delta V =  -\textbf{F} \cdot \Delta \textbf{r}=-\textbf{a} \cdot \Delta \textbf{r}=-\textbf{a} \cdot \textbf{v} \Delta t=-(a_x v_x+a_y v_y)\Delta t
```


これをオイラー法を用いることで、ポテンシャルエネルギーの推定値を求めることができる。（理論値との比較ができる）
# オイラー法
ある物体の位置は、速度を用いて近似的に以下のように表すことができる。
```math
x(t+\Delta t)=x(t)+v(t)\Delta t
```

一方で、速度$v(t)$は加速度$a(t)$を用いて近似的に、以下のように表すことができる。

```math
v(t+\Delta t)=v(t)+a(t)\Delta t
```

このことを用いて、$x(t)$を求めていく方法をオイラー法をいう。

# プログラム
上記のアルゴリズムや問題設定を用いてプログラムを書くと以下のようになる。（実行の際は、適宜コメントを外したり付けたりして欲しい。）ただし、$\alpha =-1.8$とした。

```python Euler_Planet.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#微小時間
h=1.0*10**-5
# 中心力の定数
alpha=-1.8

# 初期値
## 時間
t=0
## 位置
x=2
y=0
## 速度
v_x=0
v_y=0.5

# 計算結果を入れるためのリスト
t_ary=[]
a_x_ary=[]
a_y_ary=[]
v_x_ary=[]
v_y_ary=[]
x_ary=[]
y_ary=[]
V_ary=[]
E_ary=[]
E_0_ary=[]
L_ary=[]
V_0_ary=[]
V_1_ary=[]
V_x=0
V_y=0
V_0=0
while t<=50:
  r=math.sqrt(x**2+y**2)
  a_x=-r**alpha*x/r
  a_y=-r**alpha*y/r
  v_x=v_x+a_x*h
  v_y=v_y+a_y*h
  x=x+v_x*h
  y=y+v_y*h
  t=t+h
  # オイラー法で求めたポテンシャルエネルギー
  #V_0=-math.sqrt(V_x**2+V_y**2)
  #V_0=V_0-np.array([a_x,a_y])@(np.array([v_x,v_y])*h)
  V_0=V_0-(a_x*v_x+a_y*v_y)*h
  # 理論値のポテンシャルエネルギー
  if alpha != -1:
    V_1=r**(alpha+1)/(alpha+1)
  else:
    V_1=np.log(r)
  E=0.5*(v_x**2+v_y**2)+V_1
  E_0=0.5*(v_x**2+v_y**2)+V_0
  L=x*v_y-y*v_x
  t_ary.append(t)
  a_x_ary.append(a_x)
  a_y_ary.append(a_y)
  v_x_ary.append(v_x)
  v_y_ary.append(v_y)
  x_ary.append(x)
  y_ary.append(y)
  V_0_ary.append(V_0)
  V_1_ary.append(V_1)
  E_ary.append(E)
  E_0_ary.append(E_0)
  L_ary.append(L)

## 時間と変位の関係
# plt.plot(t_ary,x_ary,label="x")
# plt.plot(t_ary,y_ary,label="y")
# plt.legend()
# plt.xlabel("t")
# plt.ylabel("x,y")
# plt.savefig("x_t_y_t.png")
# plt.show()

## 位置関係
# plt.plot(x_ary,y_ary)
# plt.xlabel("x")
# plt.ylabel("y")
# plt.savefig("x_y.png")
# plt.show()

# 力学的エネルギーの推移
# plt.plot(t_ary,E_ary)
# plt.xlabel("t")
# plt.ylabel("E")
# plt.savefig("エネルギー2.png")
# plt.show()

## 角運動量の推移
# plt.plot(t_ary,L_ary)
# plt.xlabel("t")
# plt.ylabel("L")
# plt.savefig("角運動量.png")
# plt.show()

## ポテンシャルエネルギーの理論値とオイラー法で求めた値の比較
# plt.plot(V_1_ary,V_0_ary)
# plt.xlabel("V_1:理論値")
# plt.ylabel("V_0:オイラー法")
# plt.savefig("V_0_V_1.png")
# plt.show()

# # 力学的エネルギーの推移
plt.plot(t_ary,E_0_ary)
plt.xlabel("t")
plt.ylabel("E_0")
plt.savefig("エネルギー_推定値.png")
plt.show()
```

これを実行する（コメントを外したりは適宜行う）と以下のようなグラフらが出力される。

# 結果

## 位置と時間

![x_t_y_t.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/859d74c8-d3d0-4441-13a6-eda5454e0ad9.png)

$x,y$座標ともに複雑な振動をしているということがわかる。


## 2次元座標

![x_y.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/3f4092d1-3447-9d8b-4f97-d8792c6c46e8.png)

かなり複雑な天体の軌道を示している。

## 角速度

![角運動量.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/5e355736-a3ca-33d3-f33b-c7039c4e48f8.png)

ほとんど$0$となり保存しているようだが、オイラー法を用いているため、わずかながら誤差が生じている。


## 力学的エネルギー
### 理論値

![エネルギー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/6cd983b4-6800-219b-7359-42b0ddd73e8c.png)


### 推定値
オイラー法を用いて求めた位置エネルギーの推定値を用いたものを以下に示す。

![エネルギー_推定値.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/51cf97e8-b4c4-9fb4-62dc-9ee9cf24915e.png)

このように、系全体のエネルギーの総量が徐々に変化しているように推定されてしまっている。これは、オイラー法による誤差によるもの推測される。


## 位置エネルギー
推定値と理論値の比較を以下に示す。
![V_0_V_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/da654e99-527a-9107-ca6a-558b4aa95f73.png)


このように、$y=x$に傾きは漸近していることから、オイラー法は数値解析手法として有効であると考えられる。$y=x$からずれているものは、積分定数の違いによるものが考えられる。

## 力学的エネルギーの推定値と実際の値
最後に、力学的エネルギーの推定値と実際の値について比較する。

![エネルギー_推定値E.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/d3770a58-a87f-3e22-a863-cd985e58cfba.png)

このように、振動的な挙動を示すことがわかった。また、推定値も0.125からほぼ一定の値を示していることから、これがオイラー法による精度の限界ではないかと考えられる。

# まとめ
今回は、オイラー法を用いて天体の運動のシミュレーションを行った。オイラー法では、ある程度の誤差はでてしまうが、今回のような1天体の軌道ならば比較的誤差は少なくなると考えられる。