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

Last updated at 2025-01-20Posted at 2025-01-13

はじめに

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

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

問題設定

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

\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}

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

L=x v_y-y v_x

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

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}

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

\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

これをオイラー法を用いることで、ポテンシャルエネルギーの推定値を求めることができる。（理論値との比較ができる）

オイラー法

ある物体の位置は、速度を用いて近似的に以下のように表すことができる。

x(t+\Delta t)=x(t)+v(t)\Delta t

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

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()