二元連立一階微分方程式と二階微分方程式のコラボ

Question

前提・実現したいこと

二元連立一階微分方程式を関係式として持つ変数を含む、二階微分方程式を四元連立一階微分方程式として解きたい。

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

二元連立一階微分方程式は解かれているが、二階微分方程式だけ解けていないようです。エラーは出ていませんが二階微分方程式の方の変数（コード中ではx1）と時間tのグラフも一直線で、x1の値も全て同じ値になっていて明らかにおかしいのです。

何が起こっているのか、どうすればよいのか、どなたかアドバイスお願い致します。

エラーメッセージ

該当のソースコード

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



t_min = 0.0
t_max = 10**4


N = 10**4


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]

    
    
if (__name__ == '__main__'):
    
    t = np.linspace(t_min, t_max, N)  # 時間刻み
    
    n = 200/360       # 平均運動 n≡2π/T ～200, 100, 50, etc...[°/日]
    a = 4*10**6        # 衛星の軌道長半径 [km]
    e = 0.1
    k = n*a*pow(1-e*e, -1/2)
    
    u0 = [a, 0.0, 0.0, 1.0]     # r,f,γ,dγ/dtの初期値
    
    u = odeint(func, u0, t, args=(k, n, a, e))
    

試したこと

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

Python3をAnacondaで使っています。

Answer 1

t_max = 10**2
N = 10**6

刻み幅が大きいことが原因

Comments

[HaruKasu]

GradMaisonTokyo　様
またも助けて頂きありがとうございます！見たいものによってどれほどのスケールにするかが大事だということは漠然とわかってはいるのですが、結局具体的な値はどうするかとなると経験不足もあるのか見当がつきません。今後もいろいろと試しながら慣れていこうと思います。

[WathMorks]

計算時間幅 t_maxを出来る限り小さくとることがポイントです。

[HaruKasu]

GradMaisonTokyo　様
そうなんですか…！！このコードは惑星・衛星の相互作用による衛星の周回運動への影響を見るというの課題だったので時間もつい長くとらないといけないと思ってしまったのですが、そういうわけでもなかったのかも知れません（ケースバイケースかも知れませんが）。勉強になります。