2025年01月22日

一日一項目 天体運動計算で四苦八苦 なんですとぉーーー

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牧場
image

次にPython Scipyにodeint()という微分方程式の数値計算をする関数があるので利用して同じ結果を得た。これは地道に数式を書かないで関数使うだけでして楽勝Open-mouthed smile

そうしているうちに odeinet()とは別のsolve_ivp()が推奨だってchatGPT3.5先生が教えてくれたのでトライの試行のチャレンジィ!! 
この一行
# 数値計算
solution = solve_ivp(orbit, t_span, y0, t_eval=t_eval)
orbitというのは計算式、そして時間範囲、初期値等を与えるで~す。そして実行したら
image

おい、ズレたぞ!!  時間範囲を拡大してみると
image

  おお、蚊取り線香かーーーーーい!!    というわけで対策をchatGPT3.5先生に訊くワタクシ。

solve_ivp のオプション method を使用して、エネルギー保存性が比較的高いソルバーを選択します。例: DOP853(高精度のルンゲ=クッタ法)
solution = solve_ivp(orbit, t_span, y0, t_eval=t_eval, method='DOP853')

これで実行したら あららら つながったぞでありま~す!! 
image

他にも精度を上げる対策等が出てきたが。chatGPT3.5先生は御親切であるぞのだぉ。

今日はsolve_ivp()が推奨で odeint()は初学者用だって知っただけでいいや。

solve_ivp()については解説記事を探して読むとしようかしらのだぉでありま~す。どや。

上の天体運動の数値解であれば自作でもodeint()でも何でもいいがなや(笑) なにを!

  だって にんげんだもの byニセみつほ

posted by toinohni at 08:24| 東京 ☀| Comment(0) | 物理科学雑学 | このブログの読者になる | 更新情報をチェックする
この記事へのコメント
コメントを書く
コチラをクリックしてください