前提・実現したいこと
衛星が惑星の周りを軌道運動している時、潮汐力の影響による衛星のひょうどう運動を、連立微分方程式を解くことで数値計算したいです。
発生している問題・エラーメッセージ
(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で使っています。
回答1件
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。
2020/01/21 02:07