４元連立微分方程式で衛星の軌道計算

Question

前提・実現したいこと

衛星が惑星の周りを軌道運動している時、潮汐力の影響による衛星のひょうどう運動を、連立微分方程式を解くことで数値計算したいです。

発生している問題・エラーメッセージ

(t,r)グラフを表示してrがtに対して周期的であれば良いのですが、直線的なグラフになってしまい、修正の仕方がわかりません。エラーが出て困っているのではなく、微分方程式が正しく解けていないことはグラフから明らかだが修正箇所及び方法がわからない、という質問です。こんな時間ですが、宜しければどなたかお願い致します…！

エラーメッセージ

該当のソースコード

import math
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


day = 3000

t_min = 0.0
t_max = day*24*60*60


N = 10**5

n = 2*np.pi
a = 1.0

e = 0.1

k = n*a*pow(1-e*e, -1/2)


def func(u, t, k, n, a, e):
    
    r = u[0]
    f = u[1]
    x1 = u[2]             # γ
    x2 = u[3]             # dγ/dt
    
    drdt = k*e*np.sin(f)
    dfdt = k*(1+e*np.cos(f))/r
    dx1dt = x2
    dx2dt = -0.02*n*n*pow(a/r, 3)*np.sin(2*x1-2*f)
    
    
    return ([drdt, dfdt, dx1dt, dx2dt])

    


t = np.linspace(t_min, t_max, N)  # 時間刻み


u0 = [1.0, 10.0, 10.0, 10.0]     # r,f,γ,dγ/dtの初期値

u_list = odeint(func, u0, t, args=(k, n, a, e))


plt.plot(t, u_list[:,0])

試したこと

そもそも方程式を書き間違っていないかなどは確認済みです。

補足情報（FW/ツールのバージョンなど）

Python3をAnacondaで使っています。

Answer 1

day = 0.1
N = 10**4
に変更する。

オリジナルは時間区間の設定が広すぎます。

Comments

[HaruKasu]

GradMisonTokyo 様
ありがとうございます！グラフがうまく出ない理由をいろいろ検討したつもりでしたが、そういったことが原因のひとつになりうるということを知ることができてとても感謝しています！今後も大学の研究でプログラミングが必須の日々が続くので、上手くいかない時の検討材料のひとつに組み込みます。