f = ma = GMm/R^2 ってのから。加速度はx,yの二階微分。4次R-K法で数値計算する。概略はこんな感じ
r = sqrt(x * x Y y * y) 一般化を考えている。n = -2, a < 0とすると重力の場合
x'' = a * x * r ^ (n - 1)
y'' = a * y * r ^ (n - 1)
x ' = vx ; vx' = a * x * (x * x + y * y) ^ (n / 2 - 1 / 2)
y ' = yx ; yx' = a * y * (x * x + y * y) ^ (n / 2 - 1 / 2)
何かの本にあった。これをC++で地道に計算する。次のようなコード(一部)
Kx1 = h * func_fx(vx);
Lx1 = h * func_gx(x, y, n, a);
Kx2 = h * func_fx(vx + Lx1 / 2.0);
Lx2 = h * func_gx(x + Kx1 / 2.0, y, n, a);
Kx3 = h * func_fx(vx + Lx2 / 2.0);
Lx3 = h * func_gx(x + Kx2 / 2.0, y, n, a);
Kx4 = h * func_fx(vx + Lx3);
Lx4 = h * func_gx(x + Kx3, y, n, a);
x = x + (Kx1 + 2 * Kx2 + 2 * Kx3 + Kx4) / 6.0;
vx = vx + (Lx1 + 2 * Lx2 + 2 * Lx3 + Lx4) / 6.0;
計算結果をファイルに保存しwgnuplotに食わせると回転運動。。。 ここまではOK牧場
次にPython Scipyにodeint()という微分方程式の数値計算をする関数があるので利用して同じ結果を得た。これは地道に数式を書かないで関数使うだけでして楽勝
そうしているうちに odeinet()とは別のsolve_ivp()が推奨だってchatGPT3.5先生が教えてくれたのでトライの試行のチャレンジィ!!
この一行
# 数値計算
solution = solve_ivp(orbit, t_span, y0, t_eval=t_eval)
orbitというのは計算式、そして時間範囲、初期値等を与えるで~す。そして実行したら
おお、蚊取り線香かーーーーーい!! というわけで対策をchatGPT3.5先生に訊くワタクシ。
solve_ivp
のオプション method
を使用して、エネルギー保存性が比較的高いソルバーを選択します。例: DOP853
(高精度のルンゲ=クッタ法)
solution = solve_ivp(orbit, t_span, y0, t_eval=t_eval, method='DOP853')
他にも精度を上げる対策等が出てきたが。chatGPT3.5先生は御親切であるぞのだぉ。
今日はsolve_ivp()が推奨で odeint()は初学者用だって知っただけでいいや。
solve_ivp()については解説記事を探して読むとしようかしらのだぉでありま~す。どや。
上の天体運動の数値解であれば自作でもodeint()でも何でもいいがなや(笑) なにを!
だって にんげんだもの byニセみつほ