使用 Matlab 的 ode45 进行基本的 2 体交互
Basic 2-body interaction using Matlab's ode45
我正在尝试为围绕着巨大物体的质量可忽略不计的物体建立基本引力模型。我遵循了 ODE 套件文档中提供的示例，但我得到的结果显然很荒谬。
这是我用 ode45 调用的函数：
function dy = rigid(t,y)
dy = zeros(4,1);
%Position
xx=y(1);
yy=y(2);
%Radius
r=(xx.^2+yy.^2).^0.5;
%Constants
M=10^30;
G=6.67*10^-11;
%dX/dt
dy(1)=y(3); %vx
dy(3)=-M.*G.*xx.*r.^-3; %ax
%dY/dt
dy(2)=y(4); %vy
dy(4)=-M.*G.*yy.*r.^-3; %ay
这里是求解器行：
%Init
x_0=1;
y_0=1;
vx_0=0;
vy_0=5;
%ODEs
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]);
%Vectors
x=Y(:,1);
y=Y(:,2);
%Plot
plot(x,y)
xlabel('x');
ylabel('y');
title('y=f(x)');
最后我得到了一个线性图。即使以初始速度，位置也不会在几步内发生变化。我唯一能想到的是我误解了设置 ODE 系统的方法。
这个问题我已经琢磨了一段时间了，在网上搜索了几次，确实没什么想法。
总结：哈密顿系统与普通数值积分器的积分存在问题，你的特殊初始条件加剧了这一点，以至于数值解与正确的解不相似。
您的实施本身没有任何问题。但是，您使用的初始条件并不是最好的。您使用的常数 G 和 M 以 SI 单位表示，这意味着坐标以 m 为单位，速度以 m/s 为单位，时间以 s 为单位。这些行
x_0=1;
y_0=1;
vx_0=0;
vy_0=5;
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]);
因此意味着您要求的轨道半径约为 1.4 米，轨道速度为 5 m/s，并且您希望该轨道运行 17 分钟。想象一下，实际上有一个物体距离 10^30 公斤的质量只有几米远！
所以让我们尝试要求更现实的东西，类似于Earths' orbit，并在1年内查看它：
x_0=149.513e9;
y_0=0;
vx_0=0;
vy_0=29.78e3;
[T,Y] = ode45(@rigid,[0 31.536e6],[x_0 y_0 vx_0 vy_0]);
结果是
看起来符合预期。
但是这里还有第二个问题。让我们看看这个 10 年 ([0 315.36e6]) 的轨道：
现在我们得到的不再是闭合轨道，而是螺旋！这是因为数值积分以有限的精度进行，对于这组方程，这导致（从物理上讲）loss of energy。可以使用 ode45 的参数提高精度，但最终问题仍然存在。
现在让我们回到您的原始参数并查看结果：
这个"orbit"是一条指向原点（太阳）的直线。这可能没问题，因为直线振荡是椭圆轨道的一种可能的特例。但是随着时间的推移查看坐标：
我们看到没有振荡。 "planet" 落入太阳并停留在那里。这里发生的事情与更大轨道的效果相同：不精确的积分导致能量损失。此外，数值积分卡住了：我们要求的周期为 1000 秒，但积分不会超过 1.6e-10 秒。
据我所知，ode45 和任何其他标准 Matlab 积分器都不足以解决这个问题。有专门为哈密顿系统设计的特殊数值积分算法，称为 symplectic integrators. There is a file exchange entry that provides different implementations. Also see this question and answers 更多指针。
我正在尝试为围绕着巨大物体的质量可忽略不计的物体建立基本引力模型。我遵循了 ODE 套件文档中提供的示例，但我得到的结果显然很荒谬。
这是我用 ode45 调用的函数：
function dy = rigid(t,y)
dy = zeros(4,1);
%Position
xx=y(1);
yy=y(2);
%Radius
r=(xx.^2+yy.^2).^0.5;
%Constants
M=10^30;
G=6.67*10^-11;
%dX/dt
dy(1)=y(3); %vx
dy(3)=-M.*G.*xx.*r.^-3; %ax
%dY/dt
dy(2)=y(4); %vy
dy(4)=-M.*G.*yy.*r.^-3; %ay
这里是求解器行：
%Init
x_0=1;
y_0=1;
vx_0=0;
vy_0=5;
%ODEs
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]);
%Vectors
x=Y(:,1);
y=Y(:,2);
%Plot
plot(x,y)
xlabel('x');
ylabel('y');
title('y=f(x)');
最后我得到了一个线性图。即使以初始速度，位置也不会在几步内发生变化。我唯一能想到的是我误解了设置 ODE 系统的方法。
这个问题我已经琢磨了一段时间了，在网上搜索了几次，确实没什么想法。
总结：哈密顿系统与普通数值积分器的积分存在问题，你的特殊初始条件加剧了这一点，以至于数值解与正确的解不相似。
您的实施本身没有任何问题。但是，您使用的初始条件并不是最好的。您使用的常数 G 和 M 以 SI 单位表示，这意味着坐标以 m 为单位，速度以 m/s 为单位，时间以 s 为单位。这些行
x_0=1;
y_0=1;
vx_0=0;
vy_0=5;
[T,Y] = ode45(@rigid,[0 1000],[x_0 y_0 vx_0 vy_0]);
因此意味着您要求的轨道半径约为 1.4 米，轨道速度为 5 m/s，并且您希望该轨道运行 17 分钟。想象一下，实际上有一个物体距离 10^30 公斤的质量只有几米远！
所以让我们尝试要求更现实的东西，类似于Earths' orbit，并在1年内查看它：
x_0=149.513e9;
y_0=0;
vx_0=0;
vy_0=29.78e3;
[T,Y] = ode45(@rigid,[0 31.536e6],[x_0 y_0 vx_0 vy_0]);
结果是
看起来符合预期。
但是这里还有第二个问题。让我们看看这个 10 年 ([0 315.36e6]) 的轨道：
现在我们得到的不再是闭合轨道，而是螺旋！这是因为数值积分以有限的精度进行，对于这组方程，这导致（从物理上讲）loss of energy。可以使用 ode45 的参数提高精度，但最终问题仍然存在。
现在让我们回到您的原始参数并查看结果：
这个"orbit"是一条指向原点（太阳）的直线。这可能没问题，因为直线振荡是椭圆轨道的一种可能的特例。但是随着时间的推移查看坐标：
我们看到没有振荡。 "planet" 落入太阳并停留在那里。这里发生的事情与更大轨道的效果相同：不精确的积分导致能量损失。此外，数值积分卡住了：我们要求的周期为 1000 秒，但积分不会超过 1.6e-10 秒。
据我所知，ode45 和任何其他标准 Matlab 积分器都不足以解决这个问题。有专门为哈密顿系统设计的特殊数值积分算法，称为 symplectic integrators. There is a file exchange entry that provides different implementations. Also see this question and answers 更多指针。