matlab - Matlab，汉密尔顿系统，用 Euler 绘图和能量计算（后向、前向、中点和半隐式）

Question

这就是二体开普勒问题。

哈密​​顿量是

这应该是一个“椭圆”。当初始值仅在一个轴上具有速度分量并且起始坐标在同一轴上时，以所有神圣的名义，它怎么可能是一个椭圆？

你如何使用反向欧拉绘制数值解？

Answer 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

也许我误解了系统的实际含义，但我看到的是一个只有在彼此引力影响下才有两个物体的系统。body 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 在 y 轴上绘制 q1 和 q2，在 x 轴上绘制时间，这对我来说似乎是合乎逻辑的，如果 q1 和 q2 描述了有问题的身体与 origo 的距离一个不动的参考系统。但是如何让它看起来像一个椭圆呢？