Solar System Mapping, ODE45 trouble

Question

The following code aims to map out the solar system by incorporating every significant body's effect on the others. Of course it should result in expected orbits. In the final embedded function GravityDE it cannot read the values from PlanetVec and, because of that, cannot produce the correct new results each time. We get the error

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

Any suggestions for how to solve this would be most welcome!

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

Thanks!!

Answer 1

OK, here we go. It's going to be a LONG answer, and it will probably be over-complete, but I think it'll be valuable also for future visitors.

The force on a celestial body of mass M due to another celestial body of mass m equals

F = -GMm / r²

where r is the distance between the two bodies. That is how we all learn Newton's equation in high school (I hope, anyway...). The basic equation above is somewhat flawed, for reasons I indicated in my comment above. Also, it is incomplete in the context of more-than-two celestial bodies.

First, GM is replaced by something that can actually be measured, μ -- a body's standard gravitational parameter. Second, calculations need to be easy to do in an arbitrary coordinate system. Third, the direction of the force is not included in the equation above. Fourth, the acceleration of a celestial body is usually what matters, not the force.

All this can be included by re-phrasing the equation as

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

where boldface letters indicate vectors, r is a position vector w.r.t. some arbitrary coordinate system, r_ij = r_i - r_j is the vector from body i to body j, and the double-dots are Newton's double-fluxions (equal to Leibniz' d²/dt²).

In words: the instantaneous directed acceleration of body i due to the gravitational effects of all other bodies j, is the summation of the gravitational parameters μ_j, divided by the squared distance between body i and j, and multiplied by the vector between the two bodies scaled to unity, over all other bodies j in the system (you see why we invented Mathematical notation? :)

This equation describes a system of j second-order differential equations, which cannot be solved analytically (in an easily computable, closed form anyway) and must therefore be solved numerically. Matlab's ode45 can do that, although its accuracy leaves something to be desired if you want to simulate the planet's orbits for several dozens of years or more (especially mercury's orbit is notoriously difficult to compute accurately numerically).

Anyway, you solve this with ode45 as follows. Define

y₀ = [r₁ṙ₁r₂ṙ₂ ... r_jṙ_j ]^T

which is the collection of initial state vectors of all j bodies. In code (independent of input sizes):

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

Note that you do not need to give all the state vectors individual names. This I also strongly discourage; how would you go about integrating the asteroid belt? Would you really give ~500.000 variable names to all the individual state vectors? No--use Matlab's vector/matrix nature to your advantage.

You can probably prevent the mess above by defining your xposition and yposition (and speeds) not separately, but together in a single vector.

The ode45 integrator works by computing

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

from an input y at each iteration. It is here that the modified Newton's equation above comes into play; to compute the double-fluxions r̈_j. In code:

% 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);

The way you were doing it was to compute the double-fluxion of the position of the planet, while the other planets were in their initial positions; you wanted to pass PlanetVec, the collection of initial statevectors (that I've called y₀) into GravityDE, and compute distances and accelerations of the planet w.r.t. that vector.

This is or course incorrect -- this will only move one planet, while keeping the other planets still. That's not how the Solar system works :p Now you might argue that it doesn't matter much, since the planets move so slowly, but that is true only for the outer planets; Mercury's influence on Venus' orbit for example is grossly miscalculated that way.

Now you know the general principles. DISCLAIMER: I've written this all from memory without much checking, so I might have missed a minus sign here and there. I think it's a good exercise for you to understand and check everything I did here.

Now, a complete, functional, copy-pastable summary:

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

Note that I have not even looked at why your original function failed to run. I bet that once you get this running, it doesn't really matter anymore.