太阳系映射，ODE45 故障

Question

以下代码旨在通过结合每个重要天体对其他天体的影响来绘制太阳系。当然，它应该导致预期的轨道。在最终的嵌入函数GravityDE 中，它无法从PlanetVec 读取值，因此，每次都无法产生正确的新结果。我们得到错误

??? Undefined function 'GravityDE' for input arguments of type double. 

任何有关如何解决此问题的建议都将受到欢迎！

function Gravity1()

clear;
format long eng;
load('solar_system_data.mat');

StartTime = 0;
TimeStep = 24 * 3600 * 10;
EndTime = 24 * 3600 * 100;


TVec = StartTime:TimeStep:EndTime;
TimeStepMin = StartTime:2:TimeStep;

%Column Vectors for initial conditions
SunVec = [xposition(1), yposition(1), vx(1), vy(1),mass(1),1];
MercuryVec = [xposition(2), yposition(2), vx(2), vy(2),mass(2),2];
VenusVec = [xposition(3), yposition(3), vx(3), vy(3),mass(3),3];
EarthVec = [xposition(4), yposition(4),vx(4), vy(4),mass(4),4];
MoonVec = [xposition(10), yposition(10), vx(10), vy(10),mass(10),10];
MarsVec = [xposition(5), yposition(5), vx(5), vy(5),mass(5),5];
JupiterVec = [xposition(6), yposition(6), vx(6), vy(6),mass(6),6];
SaturnVec = [xposition(7), yposition(7), vx(7), vy(7),mass(7),7];
UranusVec = [xposition(8), yposition(8), vx(8), vy(8),mass(8),8];
NeptuneVec = [xposition(9), yposition(9), vx(9), vy(9),mass(9),9];
PlanetVec=[SunVec(1),SunVec(2),SunVec(3),SunVec(4),SunVec(5),SunVec(6);MercuryVec(1),       MercuryVec(2), MercuryVec(3), MercuryVec(4),MercuryVec(5),MercuryVec(6);VenusVec(1), VenusVec(2), VenusVec(3), VenusVec(4),VenusVec(5),VenusVec(6);EarthVec(1), EarthVec(2), EarthVec(3), EarthVec(4),EarthVec(5),EarthVec(6);MoonVec(1),MoonVec(2),MoonVec(3),MoonVec(4),MoonVec(5),MoonVec(6);MarsVec(1), MarsVec(2), MarsVec(3), MarsVec(4),MarsVec(5),MarsVec(6);JupiterVec(1), JupiterVec(2), JupiterVec(3), JupiterVec(4),JupiterVec(5),JupiterVec(6);SaturnVec(1), SaturnVec(2), SaturnVec(3), SaturnVec(4),SaturnVec(5),SaturnVec(6);UranusVec(1), UranusVec(2),UranusVec(3), UranusVec(4),UranusVec(5),UranusVec(6);NeptuneVec(1),       NeptuneVec(2), NeptuneVec(3), NeptuneVec(4),NeptuneVec(5),NeptuneVec(6)];
n=0;
while n<EndTime;
%Built in solver
[TimeVec, SunMat] = ode45(@GravityDE, TimeStepMin, SunVec);
[TimeVec, MercuryMat] = ode45(@GravityDE, TimeStepMin, MercuryVec);
[TimeVec, VenusMat] = ode45(@GravityDE, TimeStepMin, VenusVec);
[TimeVec, EarthMat] = ode45(@GravityDE, TimeStepMin, EarthVec);
[TimeVec, MoonMat] = ode45(@GravityDE, TimeStepMin, MoonVec);
[TimeVec,  MarsMat] = ode45(@GravityDE, TimeStepMin, MarsVec);
[TimeVec, JupiterMat] = ode45(@GravityDE, TimeStepMin, JupiterVec);
[TimeVec, SaturnMat] = ode45(@GravityDE, TimeStepMin, SaturnVec);
[TimeVec,  UranusMat] = ode45(@GravityDE, TimeStepMin, UranusVec);
[TimeVec, NeptuneMat] = ode45(@GravityDE, TimeStepMin, NeptuneVec);

SunXVec = SunMat (end,1);
SunYVec = SunMat (end,2);
SunVXVec = SunMat(end,3);
SunVYVec = SunMat(end,4);

MercuryXVec = MercuryMat (end,1);
MercuryYVec = MercuryMat (end,2);
MercuryVXVec = MercuryMat(end,3);
MercuryVYVec = MercuryMat(end,4);

VenusXVec = VenusMat (end,1);
VenusYVec = VenusMat (end,2);
VenusVXVec = VenusMat(end,3);
VenusVYVec = VenusMat(end,4);

EarthXVec = EarthMat (end,1);
EarthYVec = EarthMat (end,2);
EarthVXVec = EarthMat(end,3);
EarthVYVec = EarthMat(end,4);

MoonXVec = MoonMat (end,1);
MoonYVec = MoonMat (end,2);
MoonVXVec = MoonMat(end,3);
MoonVYVec =MoonMat(end,4);

MarsXVec = MarsMat (end,1);
MarsYVec = MarsMat (end,2);
MarsVXVec = MarsMat(end,3);
MarsVYVec = MarsMat(end,4);

JupiterXVec = JupiterMat (end,1);
JupiterYVec = JupiterMat (end,2);
JupiterVXVec = JupiterMat(end,3);
JupiterVYVec =JupiterMat(end,4);

SaturnXVec = SaturnMat (end,1);
SaturnYVec = SaturnMat (end,2);
SaturnVXVec = SaturnMat(end,3);
SaturnVYVec =SaturnMat(end,4);

UranusXVec = UranusMat (end,1);
UranusYVec = UranusMat (end,2);
UranusVXVec = UranusMat(end,3);
UranusVYVec =UranusMat(end,4);

NeptuneXVec = NeptuneMat (end,1);
NeptuneYVec = NeptuneMat (end,2);
NeptuneVXVec = NeptuneMat(end,3);
NeptuneVYVec =NeptuneMat(end,4);

SunVec=[SunXVec,SunYVec,SunVXVec,SunVYVec,mass(1),1];
MercuryVec = [MercuryXVec, MercuryYVec, MercuryVXVec, MercuryVYVec,mass(2),2];
VenusVec = [VenusXVec, VenusYVec, VenusVXVec, VenusVYVec,mass(3),3];
EarthVec = [EarthXVec, EarthYVec, EarthVXVec, EarthVYVec,mass(4),4];
MoonVec = [MoonXVec,MoonYVec,MoonVXVec,MoonVYVec,mass(10),10];
MarsVec = [MarsXVec, MarsYVec, MarsVXVec, MarsVYVec,mass(5),5];
JupiterVec = [JupiterXVec, JupiterYVec, JupiterVXVec, JupiterVYVec,mass(6),6];
SaturnVec = [SaturnXVec, SaturnYVec, SaturnVXVec, SaturnVYVec,mass(7),7];
UranusVec = [UranusXVec, UranusYVec,UranusVXVec, UranusVYVec,mass(8),8];
NeptuneVec = [NeptuneXVec, NeptuneYVec, NeptuneVXVec, NeptuneVYVec,mass(9),9];
PlanetVec=[SunVec(1),SunVec(2),SunVec(3),SunVec(4),SunVec(5),SunVec(6);MercuryVec(1),                MercuryVec(2), MercuryVec(3), MercuryVec(4),MercuryVec(5),MercuryVec(6);VenusVec(1),  VenusVec(2), VenusVec(3), VenusVec(4),VenusVec(5),VenusVec(6);EarthVec(1), EarthVec(2),   EarthVec(3),   EarthVec(4),EarthVec(5),EarthVec(6);MoonVec(1),MoonVec(2),MoonVec(3),MoonVec(4),MoonVec(5),MoonVec(6);MarsVec(1), MarsVec(2), MarsVec(3), MarsVec(4),MarsVec(5),MarsVec(6);JupiterVec(1),  JupiterVec(2), JupiterVec(3), JupiterVec(4),JupiterVec(5),JupiterVec(6);SaturnVec(1),  SaturnVec(2), SaturnVec(3), SaturnVec(4),SaturnVec(5),SaturnVec(6);UranusVec(1),  UranusVec(2),UranusVec(3), UranusVec(4),UranusVec(15),UranusVec(6);NeptuneVec(1),  NeptuneVec(2), NeptuneVec(3), NeptuneVec(4),NeptuneVec(5),NeptuneVec(6)];

plot (SunXVec,SunYVec,'.','Color','yellow');
hold on;
plot (MercuryXVec,MercuryYVec,'.','Color','green');
hold on;
plot (VenusXVec,VenusYVec,'.','Color','blue');
hold on;
plot (EarthXVec,EarthYVec, '.','Color', 'red');
hold on;
plot (MoonXVec,MoonYVec, '.','Color','black');
hold on;
plot (MarsXVec,MarsYVec, '.','Color','black');
hold on;
plot (JupiterXVec,JupiterYVec,'.','Color','green');
hold on;
plot (SaturnXVec,SaturnYVec, '.','Color','blue');
hold on;
plot (UranusXVec,UranusYVec, '.','Color','red');
hold on;
plot (NeptuneXVec,NeptuneYVec, '.','Color','blue');
hold on;
n=n+TimeStep;
end

function dYVec = GravityDE (TimeStep, YVec,PlanetVec)
load('solar_system_data.mat');
GravConst = 6.67259e-11;
Xi = YVec(1);
Yi = YVec(2);
VXi = YVec(3);
VYi = YVec(4);
Massi=YVec(5);
BodyName=YVec(6);

AccXtotal=0;
AccYtotal=0;
j=1;
while j<=10

Massj=PlanetVec(j,5);
Yj=PlanetVec(j,2);
Xj=PlanetVec(j,1);



RangeSq = (Xi-Xj).^2 + (Yi-Yj).^2;
if RangeSq==0
    AccMag=0;
    Theta = atan2(Yi-Yj,Xi-Xj);
    AccX = -AccMag .* cos (Theta);
    AccY = -AccMag .* sin (Theta);
    j=j+1;
    AccXtotal=AccXtotal+AccX;
    AccYtotal=AccYtotal+AccY;
else
    Theta = atan2(Yi-Yj,Xi-Xj);
    AccMag = (GravConst .* Massj ./ RangeSq);
    AccX = -AccMag .* cos (Theta);
    AccY = -AccMag .* sin (Theta);
    j=j+1;
    AccXtotal=AccXtotal+AccX;
    AccYtotal=AccYtotal+AccY;
    VXi=VXi+AccXtotal.*TimeStep;
    VYi=VYi+AccYtotal.*TimeStep;
end

dYVec = [VXi; VYi; AccXtotal; AccYtotal;Massi;BodyName];

end

谢谢！！

Answer 1

好的，我们开始吧。这将是一个 LONG 答案，并且可能会过于完整，但我认为它对未来的访问者也很有价值。

另一个天体m对一个天体M的作用力等于

F = -GMm / r²

其中 r 是两个物体之间的距离。这就是我们在高中学习牛顿方程的方式（我希望，无论如何......）。由于我在上面的评论中指出的原因，上面的基本等式有些缺陷。此外，在两个以上天体的情况下，它是不完整的。

首先，GM 被替换为可以实际测量的东西，μ——身体的standard gravitational parameter。其次，计算需要在任意坐标系中易于进行。第三，力的方向不包括在上面的等式中。第四，天体的加速度通常是重要的，而不是力。

所有这些都可以通过将等式重新表述为

r̈_i = Σ_j≠iμ_j r_ij / |r_ij|³

其中粗体字母表示向量，r 是位置向量 w.r.t。一些任意坐标系，r_ij = r_i - r _j 是体 i 到体 j 的向量，双点是牛顿双通量（等于莱布尼茨的 d²/dt²）。

换句话说：由于所有其他物体j的引力作用，物体i的瞬时定向加速度，是引力参数μ_j，除以主体 i 和 j 之间的平方距离，再乘以缩放为单位的两个主体之间的向量，超过系统中的所有其他物体j（你明白我们为什么发明数学符号了吗？:)

这个方程描述了一个由 j 个二阶微分方程组成的系统，它不能解析求解（无论如何以易于计算的封闭形式），因此必须用数值​​求解。 Matlab 的ode45 可以做到这一点，尽管如果您想模拟数十年或更长时间的行星轨道（尤其是众所周知的水星轨道很难用数值精确计算），它的准确性仍有待改进。

无论如何，您可以使用ode45 解决此问题，如下所示。定义

y₀ = [r₁ṙ_{1r2ṙ2 ... rj ṙj ]^T}

是所有j体的初始状态向量的集合。在代码中（与输入大小无关）：

y0 = [xposition(:) yposition(:) vx(:) vy(:)].';  
y0 = y0(:);

请注意，您确实不需要需要为所有状态向量单独命名。这我也强烈反对；你将如何整合小行星带？你真的会给所有单独的状态向量提供大约 500.000 个变量名称吗？不——使用 Matlab 的向量/矩阵特性对您有利。

您可以通过定义 xposition 和 yposition（和速度）而不是单独定义，而是一起定义在一个向量中来防止上述混乱。

ode45 积分器通过计算工作

ẏ = [ṙ₁r̈₁ṙ₂r̈₂ ... ṙ_jr̈ _j ]^T

来自每次迭代的输入 y。正是在这里，上面的修改后的牛顿方程开始发挥作用。计算双通量r̈_j。在代码中：

% collect data     
r = [y(1:4:end) y(2:4:end)];  % X/Y positions
V = [y(3:4:end) y(4:4:end)];  % Vx/Vy speeds

% initialize output
ydotdot = zeros(size(y));
ydotdot(1:4:end) = V(:,1);  % we already know the first half;
ydotdot(2:4:end) = V(:,2);  % it's simply equal to the speeds

% Compute all accelerations. 
% This is where ALL the computational burden is -- when optimizing 
% for speed, this is where to start!
sz   = size(r,1);
accx = zeros(sz);
accy = zeros(sz);
for ii = 1:sz
    ri = r(ii,:);
    for jj = ii+1:sz
        rij = ri - r(jj,:);
        sc  = (rij*rij.')^(-3/2);
        accx(jj,ii) = rij(1) * sc;
        accy(jj,ii) = rij(2) * sc;
    end
end

accx = bsxfun(@times, -mu(:), accx-accx.');
accy = bsxfun(@times, -mu(:), accy-accy.');

% insert accelerations
ydotdot(3:4:end) = sum(accx);
ydotdot(4:4:end) = sum(accy);

你这样做的方式是计算行星位置的双通量，而其他行星处于它们的初始位置；你想将PlanetVec，初始状态向量的集合（我称之为y₀）传递到GravityDE，并计算行星w.r.t的距离和加速度。那个向量。

这当然是不正确的——这只会移动一个行星，而使其他行星保持静止。太阳系不是这样工作的：p 现在你可能会争辩说这并不重要，因为行星移动得如此缓慢，但这仅适用于外行星；例如，水星对金星轨道的影响就被严重误算了。

现在您知道了一般原则。 免责声明：我都是凭记忆写的，没有做太多检查，所以我可能在这里和那里漏掉了一个减号。我认为这是一个很好的练习，让您了解和检查我在这里所做的一切。

现在，一个完整的、实用的、可复制粘贴的摘要：

function Gravity1

    % NOTE: 'clear' has no meaning at the start of a function; a function
    % has its own variable space, meaning it is empty to begin with.

    % NOTE: this assumes you put have all the planetary mu's inside this datafile
    load('solar_system_data.mat');

    % NOTE: ode45 chooses its own time steps; its an adaptive method.
    % Passing it custom time steps is hopelessly inefficient.
    t0   = 0;
    tend = 24 * 3600 * 100;

    % re-format initial vectors
    y0 = [xposition(:) yposition(:) vx(:) vy(:)].';  
    y0 = y0(:);

    % perform integration 
    [t y] = ode45(@d2ydt2, [t0 tend], y0);

    % and do plot
    h = figure; hold on % NOTE: only a single 'hold on' is needed to turn it on :) 
    plot(y(:,1),y(:,2), 'y.');    % Sun
    plot(y(:,1),y(:,2), 'g.');    % Mercury
    plot(y(:,1),y(:,2), 'b.');    % Venus
    plot(y(:,1),y(:,2), 'r.');    % Earth
    plot(y(:,1),y(:,2), 'k.');    % Mars
    plot(y(:,1),y(:,2), 'g.');    % Jupiter
    plot(y(:,1),y(:,2), 'b.');    % Saturn
    plot(y(:,1),y(:,2), 'r.');    % Uranus
    plot(y(:,1),y(:,2), 'b.');    % Neptune


    % It's easiest to put the differential equation in a nested function    
    function ydotdot = d2ydt2(~,y)

        % rename data     
        r = [y(1:4:end) y(2:4:end)];  % X/Y positions
        V = [y(3:4:end) y(4:4:end)];  % Vx/Vy speeds

        % initialize output
        ydotdot = zeros(size(y));
        ydotdot(1:4:end) = V(:,1);  % we already know the first half;
        ydotdot(2:4:end) = V(:,2);  % it's simply equal to the speeds

        % Compute all accelerations. 
        % This is where ALL the computational burden is -- when optimizing 
        % for speed, this is where to start!
        sz   = size(r,1);
        accx = zeros(sz);
        accy = zeros(sz);
        for ii = 1:sz
            ri = r(ii,:);
            for jj = ii+1:sz
                rij = ri - r(jj,:);
                sc  = (rij*rij.')^(-3/2);
                accx(jj,ii) = rij(1) * sc;
                accy(jj,ii) = rij(2) * sc;
            end
        end

        accx = bsxfun(@times, -mu(:), accx-accx.');
        accy = bsxfun(@times, -mu(:), accy-accy.');

        % insert accelerations
        ydotdot(3:4:end) = sum(accx);
        ydotdot(4:4:end) = sum(accy);

    end

end

请注意，我什至没有查看您的原始函数为何无法运行。我敢打赌，一旦你开始运行它，它就不再重要了。

Answer 2

没有solar_system_data.mat，很难调试您的代码，但问题是您没有将PlanetVec 传递给您的函数

如果我没记错的话，你应该试试

[TimeVec, SunMat] = ode45(@(t,y)GravityDE(t,y,PlanetVec), TimeStepMin, SunVec);

请阅读parameterizing functions