下ごしらえした万有引力の2体問題の運動方程式を Python で数値的に解き gnuplot でグラフにする

下ごしらえした万有引力の2体問題の運動方程式を Maxima で数値的に解く」の Python & gnuplot 版。 続きを読む

下ごしらえした万有引力の2体問題の運動方程式を Maxima で数値的に解く

無次元化された運動方程式

万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ」にまとめたように,系に特徴的な長さと時間スケールで無次元化した変数に対する運動方程式は以下のようになるのであった。

\begin{eqnarray}
\frac{d^2 X}{dT^2} = – \frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{d^2 Y}{dT^2} = – \frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}

初期条件

無次元化された変数に対する初期条件は $T = 0$ で

\begin{eqnarray}
X &=& 1 \\
Y &=& 0\\
V_x &=& \frac{dX}{dT} = 0 \\
V_y &=& \frac{dY}{dT} = \frac{v_{y0}}{v_0} = k \ \ ( 0 < k < \sqrt{2})
\end{eqnarray}

楕円の軌道要素(長半径,離心率,周期)

このとき,数値計算する前にもう,軌道は規格化された長半径が

$$ A \equiv \frac{a}{x_0} = \frac{1}{2 – k^2} $$

離心率が

$$e = |k^2 – 1|$$

の楕円になることがわかってしまう。また,規格化された周期が

$$P \equiv \frac{p}{\tau} = 2 \pi A^{\frac{3}{2}} = \frac{2 \pi}{\left(2 – k^2\right)^{\frac{3}{2}}}$$

となることもわかってしまっているのであった。

じゅあ,そこまでわかっているのなら,なぜわざわざ数値計算するのかというと,それは楕円軌道の解が時間 $t$ の陽関数として与えられているわけではないから。時刻 $t$ のとき,どこにいるかが解析的にわかっていないので,それを知りたいから数値計算することになる。

続きを読む

万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ

万有引力の2体問題(は結局1体問題に帰着するんだけど)の運動方程式を数値的に解く前の下ごしらえ。闇雲な初期条件からはじめるのではなく,そもそも数値計算しなくても楕円になることはわかっているのだから,ルンゲ・クッタ法なりを使って数値的に解く前に,それなりの下ごしらえをしておこう。

続きを読む

コメントへの対応

続きを読む

Maxima で領域の塗りつぶしの追記

Maxima で領域の塗りつぶしと回転体の表示 に書いた Maxima で領域を塗りつぶす例の追記。

当時は,塗りつぶした後に explicit(f(x), x, x0, x1) ができなくて,追加の線は parametric で描く などと嘯いていたが,以下のように filled_func = false とすれば塗りつぶし終了でその後は explicit(f(x), x, x0, x1) できます。
続きを読む

ケプラーの第1法則・第2法則のパラパラアニメ

Maxima でケプラー方程式を逐次近似的に解いて draw2d で描いてパラパラアニメにする」のバリエーション。学生の演習用。出来上がりの例として。

続きを読む

高さ h からの斜方投射の最大到達距離は角度45°のときではない

参考サイト

上記のページを参考に,地面から高さ \(h\) の地点から空気抵抗なしの斜方投射を行うと,水平方向の到達距離が最大となるのは打ち上げ角度(仰角)が \(45^{\circ}\) のときではなく,それより少し小さめになることをおさらいしておく。学生のコンピュータ演習用にと思ったが,けっこう手間取ったのでメモ。

なお,空気抵抗がある場合の斜方投射については,以下を参照。

続きを読む

gnuplot でケプラー方程式の近似解と数値解

Maxima でケプラー方程式の近似解と数値解」の gnuplot 版。 続きを読む

Python の SymPy でケプラー方程式の近似解と数値解

Maxima でケプラー方程式の近似解と数値解」の SymPy 版。 続きを読む

Maxima でケプラー方程式の近似解と数値解

続きを読む

Load more