matlab - Matlab、ハミルトニアン システム、オイラーによるプロットとエネルギー計算 (後方、前方、中点、および半陰的)

Question

これが二体ケプラー問題です。

ハミルトニアンは

これは「楕円」であるはずです。初期値が 1 つの軸上に速度成分のみを持ち、開始座標が同じ軸上にある場合、どのようにして楕円になることができますか?

後方オイラーを使用して数値解をどのようにプロットしますか?

Answer 1

2 つの 2 次方程式を 4 つの 1 次方程式に変換し、ode23(逆オイラー法ではなく、別の数値法を使用して) を使用して解きます。

odefun = @(t, y) [y(3);
                  y(4);
                  -y(1) / (y(1)^2 + y(2)^2)^(3/2);
                  -y(2) / (y(1)^2 + y(2)^2)^(3/2)]
p0 = [1 0];    % initial position
v0 = [0 1];    % initial velocity
[t y] = ode23(odefun, [0 20], [p0 v0])
comet(y(:,1), y(:,2))

トレースされたパスが楕円であることがわかります (実際には私の初速度の円です)。初期条件で遊ぶことができます（v0 = [0 1.1]楕円が表示されます）。

Answer 2

システムが実際に何を意味するのかを誤解しているかもしれませんが、私が見ているのは、互いの重力の影響下にある 2 つのボディを持つシステムです。ボディ 1 の開始位置はベクトルではなくスカラーのみであるため、X 軸上にのみ配置できます。この場合、座標は x = 1-a (0 < a < 1) です。ボディ 2 の開始位置は x=0 です。現在、重力勾配は常に同じ軸に平行であるため、x 軸から遠ざかる力のベクトルはありません。もちろん、開始速度もスカラーであるため、x 軸に平行な軌道以外の場所を指す初期速度はありません。

「これは楕円を形成する」と言うとき、私は間違った座標系を見ているのかもしれません。次のように順オイラープロットを行いました。

clear all, close all, clc

stl= 0.0005; % steplength
nos = 50000; % number of steps

q1 = zeros(1,nos);
q2 = zeros(1,nos);
q1v = zeros(1,nos);
q2v = zeros(1,nos);
q1a = zeros(1,nos);
q2a = zeros(1,nos);

% Starting positions
a=0.5;
q1(1) = 1-a;
q2(1) = 0;
q1v(1) = 0;
q2v(1) = sqrt((1+a)/(1-a));

% P is inertia, and has the same vector components as the body velocity

for i=1:nos-1    

[q1a(i),q2a(i)] = u1FUNC(q1(i),q2(i)); 

q1v(i+1) = q1v(i) + q1a(i).*stl; 
q2v(i+1) = q2v(i) + q2a(i).*stl;

q1(i+1) = q1(i) + q1v(i+1)*stl;
q2(i+1) = q2(i) + q2v(i+1)*stl;  

end

plot(q1)
hold on
plot(q2,'r')

u1FUNC

function [a,b] = func(c,d)

% [a,b] = [q1'',q2'']
% [c,d] = [q1,q2]

a = -c/(((c^2)+(d^2))^(3/2));
b = -d/(((c^2)+(d^2))^(3/2));

end

ここで、このコードは q1 と q2 が q1 と q2 を y 軸に、時間を x 軸にプロットしたプロットを提供します。動かない参照系。しかし、どうすれば楕円のように見えるのでしょうか?