<p>Here is what I did: </p>

<ul>
<li>Based on their masses, it is safest to initially consider Jupiter and Saturn as well as Uranus. It might also be fruitful to include the Earth in the analysis, to get relative positions, observation angles, etc. So, I will be considering: 

<ul>
<li>Sun</li>
<li>Earth</li>
<li>Jupiter</li>
<li>Saturn </li>
<li>Uranus</li>
<li>Neptune</li>
</ul></li>
<li>Get the <a href="https://en.wikipedia.org/wiki/Standard_gravitational_parameter" rel="noreferrer">standard gravitational parameters</a> (μ) for all of them</li>
<li>Get initial positions and velocities via <a href="http://ssd.jpl.nasa.gov/horizons.cgi" rel="noreferrer">JPL/HORIZONS</a> for all these planets. I had some data lying around from J2000.5, so I just used the state vectors from the 1st of January, 2000, at noon. </li>
<li>Write an N-body integrator with built-in MATLAB tools. Integrate this incomplete Solar system once without Neptune, and once with Neptune included. </li>
<li>Analyze and compare!</li>
</ul>

<p>So, here's my data, and N-body integrator: </p>

<pre><code>function [t, yout_noNeptune, yout_withNeptune] = discover_Neptune()

    % Time of integration (in years)
    tspan = [0 97] * 365.25 * 86400;

    % std. gravitational parameters [km/s²/kg]
    mus_noNeptune = [1.32712439940e11; % Sun
                     398600.4415       % Earth
                     1.26686534e8      % Jupiter
                     3.7931187e7       % Saturn
                     5.793939e6];      % Uranus

    mus_withNeptune = [mus_noNeptune
                       6.836529e6]; % Neptune

    % Initial positions [km] and velocities [km/s] on 2000/Jan/1, 00:00
    % These positions describe the barycenter of the associated system,
    % e.g., sJupiter equals the statevector of the Jovian system barycenter.
    % Coordinates are expressed in ICRF, Solar system barycenter
    sSun     = [0 0 0 0 0 0].';
    sEarth   = [-2.519628815461580E+07  1.449304809540383E+08 -6.175201582312584E+02,...
                -2.984033716426881E+01 -5.204660244783900E+00  6.043671763866776E-05].';
    sJupiter = [ 5.989286428194381E+08  4.390950273441353E+08 -1.523283183395675E+07,...
                -7.900977458946710E+00  1.116263478937066E+01  1.306377465321731E-01].';
    sSaturn  = [ 9.587405702749230E+08  9.825345942920649E+08 -5.522129405702555E+07,...
                -7.429660072417541E+00  6.738335806405299E+00  1.781138895399632E-01].';
    sUranus  = [ 2.158728913593440E+09 -2.054869688179662E+09 -3.562250313222718E+07,...
                 4.637622471852293E+00  4.627114800383241E+00 -4.290473194118749E-02].';
    sNeptune = [ 2.514787652167830E+09 -3.738894534538290E+09  1.904284739289832E+07,...
                 4.466005624145428E+00  3.075618250100339E+00 -1.666451179600835E-01].';

    y0_noNeptune   = [sSun; sEarth; sJupiter; sSaturn; sUranus];
    y0_withNeptune = [y0_noNeptune; sNeptune];

    % Integrate the partial Solar system 
    % once with Neptune, and once without
    options = odeset('AbsTol', 1e-8,...
                     'RelTol', 1e-10);

    [t, yout_noNeptune]   = ode113(@(t,y) odefcn(t,y,mus_noNeptune)  , tspan, y0_noNeptune  , options);
    [~, yout_withNeptune] = ode113(@(t,y) odefcn(t,y,mus_withNeptune),     t, y0_withNeptune, options);

end

% The differential equation 
%
%    dy/dt = d/dt [r₀ v₀ r₁ v₁ r₂ v₂ ... rₙ vₙ]    
%          = [v₀ a₀ v₁ a₁ v₂ a₂ ... vₙ aₙ]    
%
%  with 
%
%    aₓ = Σₘ -G·mₘ/|rₘ-rₓ|² · (rₘ-rₓ) / |rₘ-rₓ| 
%       = Σₘ -μₘ·(rₘ-rₓ)/|rₘ-rₓ|³  
%
function dydt = odefcn(~, y, mus)

    % Split up position and velocity
    rs = y([1:6:end; 2:6:end; 3:6:end]);
    vs = y([4:6:end; 5:6:end; 6:6:end]);

     % Number of celestial bodies
    N = size(rs,2);

    % Compute interplanetary distances to the power -3/2
    df  = bsxfun(@minus, permute(rs, [1 3 2]), rs);
    D32 = permute(sum(df.^2), [3 2 1]).^(-3/2);
    D32(1:N+1:end) = 0; % (remove infs)

    % Compute all accelerations     
    as = -bsxfun(@times, mus.', D32);              % (magnitudes)    
    as = bsxfun(@times, df, permute(as, [3 2 1])); % (directions)    
    as = reshape(sum(as,2), [],1);                 % (total)

    % Output derivatives of the state vectors
    dydt = y;
    dydt([1:6:end; 2:6:end; 3:6:end]) = vs;
    dydt([4:6:end; 5:6:end; 6:6:end]) = as;

end
</code></pre>

<p>Here is the driver script I used to get some nice plots out: </p>

<pre><code>clc
close all

% Get coordinates from N-body simulation
[t, yout_noNeptune, yout_withNeptune] = discover_Neptune();

% For plot titles etc.
bodies = {'Sun'
          'Earth'
          'Jupiter'
          'Saturn'
          'Uranus'
          'Neptune'};


% Extract positions
rs_noNeptune   = yout_noNeptune  (:, [1:6:end; 2:6:end; 3:6:end]);
rs_withNeptune = yout_withNeptune(:, [1:6:end; 2:6:end; 3:6:end]);



% Figure of the whole Solar sysetm, just to check
% whether everything went OK
figure, clf, hold on
for ii = 1:numel(bodies)
    plot3(rs_withNeptune(:,3*(ii-1)+1),...
          rs_withNeptune(:,3*(ii-1)+2),...
          rs_withNeptune(:,3*(ii-1)+3),...
          'color', rand(1,3));
end

axis equal
legend(bodies);
xlabel('X [km]');
ylabel('Y [km]');
title('Just the Solar system, nothing to see here');


% Compare positions of Uranus with and without Neptune
rs_Uranus_noNeptune   = rs_noNeptune  (:, 13:15);
rs_Uranus_withNeptune = rs_withNeptune(:, 13:15);

figure, clf, hold on

plot3(rs_Uranus_noNeptune(:,1),...
      rs_Uranus_noNeptune(:,2),...
      rs_Uranus_noNeptune(:,3),...
      'b.');

plot3(rs_Uranus_withNeptune(:,1),...
      rs_Uranus_withNeptune(:,2),...
      rs_Uranus_withNeptune(:,3),...
      'r.');

axis equal
xlabel('X [km]');
ylabel('Y [km]');
legend('Uranus, no Neptune',...
       'Uranus, with Neptune');


% Norm of the difference over time
figure, clf, hold on

rescaled_t = t/365.25/86400;

dx = sqrt(sum((rs_Uranus_noNeptune - rs_Uranus_withNeptune).^2,2));
plot(rescaled_t,dx);
xlabel('Time [years]');
ylabel('Absolute offset [km]');
title({'Euclidian distance between'
       'the two Uranuses'});


% Angles from Earth
figure, clf, hold on

rs_Earth_noNeptune   = rs_noNeptune  (:, 4:6);
rs_Earth_withNeptune = rs_withNeptune(:, 4:6);

v0 = rs_Uranus_noNeptune   - rs_Earth_noNeptune;
v1 = rs_Uranus_withNeptune - rs_Earth_withNeptune;

nv0 = sqrt(sum(v0.^2,2));
nv1 = sqrt(sum(v1.^2,2));

dPhi = 180/pi * 3600 * acos(min(1,max(0, sum(v0.*v1,2) ./ (nv0.*nv1) )));
plot(rescaled_t, dPhi);

xlabel('Time [years]');
ylabel('Separation [arcsec]')
title({'Angular separation between the two'
       'Uranuses when observed from Earth'});
</code></pre>

<p>which I'll describe here step by step. </p>

<p>First, a plot of the Solar system to check that the N-body integrator works as it should: </p>

<p><a href="https://i.sstatic.net/XcDzN.png" rel="noreferrer"><img src="https://i.sstatic.net/XcDzN.png" alt="the Solar system"></a></p>

<p>Nice! Next, I wanted to see the difference between the positions of Uranus with and without the influence of Neptune. So, I extracted just the positions of those two Uranuses, and plotted them: </p>

<p><a href="https://i.sstatic.net/erO5x.png" rel="noreferrer"><img src="https://i.sstatic.net/erO5x.png" alt="Two Uranuses, with and without Neptune"></a></p>

<p>...that's hardly useful. Even when zooming in greatly and rotating the heck out of it, this is just not a useful plot. So I looked at the evolution of the absolute Euclidian distance between the two Uranuses: </p>

<p><a href="https://i.sstatic.net/oL7G7.png" rel="noreferrer"><img src="https://i.sstatic.net/oL7G7.png" alt="Time evolution of the Euclidian distance between the two Uranuses"></a></p>

<p>That's starting to look more like it! Roughly 80 years after the start of our analysis, the two Uranuses are almost 6 million km apart! </p>

<p>Large as that may sound, in the grander scale of things this might drown in the noise when we take measurements here on Earth. Plus, it still does not tell the whole story, as we'll see in a moment. So next, let's look at the angular difference between the observation vectors, from Earth towards the two Uranuses to see how large that angle is, and if it can stand out above the observational error thresholds: </p>

<p><a href="https://i.sstatic.net/xwZlg.png" rel="noreferrer"><img src="https://i.sstatic.net/xwZlg.png" alt="Angular separation between the two Uranuses"></a></p>

<p>...whoa! Well over 300 arcseconds difference, plus all sorts of wobbly bobbley timey wimey rippling going on. That seems well within the observational capabilities of the time (although I can't find a reliable source on this so quickly; anyone?)</p>

<p>Just for good measure, I also produced that last plot leaving Jupiter and Saturn out of the picture. Although some perturbation theory had been developed in the 17<sup>th</sup> and 18<sup>th</sup> centuries, it was not very well developed and I doubt even Le Verrier took Jupiter into consideration (but again, I could be wrong; please correct me if you know more). </p>

<p>So, here is the last plot without Jupiter and Saturn: </p>

<p><a href="https://i.sstatic.net/GaLrT.png" rel="noreferrer"><img src="https://i.sstatic.net/GaLrT.png" alt="Angular separation between the two Uranuses, leaving Jupiter and Saturn out of the equation"></a></p>

<p>Although there are differences, they are minute, and most importantly irrelevant for discovering Neptune. </p>

CC BY-SA 3.0
2016-12-20T20:36:49.703Z
Brilliant answer!
zephyr
CC BY-SA 4.0
2020-09-24T23:47:29.117Z
@allthemikeysaretaken no, there's no source other than my code (which could also be wrong of course; I don't remember doing many validation tests on this to be honest). Those "squigglies" are the result of the Earth's orbit, and they _should_ therefore have a period of very close to 1 year. Eyeballing your first plot and my plot would indicate that your Earth is orbiting twice as slow as mine...Also, the picture you linked to uses different initial conditions (they put Uranus directly opposite Neptune), so we can't directly compare mine to theirs...
Rody Oldenhuis
CC BY-SA 4.0
2020-09-24T23:48:41.367Z
@allthemikeysaretaken Also: what exactly do you call "the" leap-frog integrator in Python? Verlet integration, though symplectic, may suffer from other sources of error that a variable-order multistep Adams-Bashforth-Moulton scheme (`ode113` in MATLAB) deals better with. It really depends on specifics whether one outperforms the other
Rody Oldenhuis
CC BY-SA 3.0
2016-12-20T16:14:36.600Z
2016-12-20T18:57:43.217Z
2016-12-20T18:57:43.217Z
Rody Oldenhuis
<p>I would like to demonstrate the existence of another planet (Neptune) by studying the discrepancy between the observation of the Uranus orbit and the mathematical prediction, this work was made from Le Verrier and I would like to understand his method.</p>

<p>I have read the Chapter 2, "The Discovery of Neptune (1845-1846)," in the biography Le Verrier -- Magnificent and Detestable Astronomer,
but it is too much in-depth and I did not understand very well his work.</p>

<p>I am studying the three-body problem (Sun, Uranus, Neptune) via Matlab and the two body problem (Sun, Uranus) taking the initial condition from here:</p>

<p><a href="http://nssdc.gsfc.nasa.gov/planetary/factsheet/uranusfact.html" rel="nofollow noreferrer">http://nssdc.gsfc.nasa.gov/planetary/factsheet/uranusfact.html</a></p>

<p>I have tried this method:
I put Uranus in the Perihelion with the Max. orbital velocity and I calculate the semi-major axis, and it is more accurate than the one obtained from putting Uranus and Neptune in the Perihelion with their respective Max. orbital velocity.</p>

<p>Here a cool pic made with Matlab:
<a href="https://i.sstatic.net/gxLML.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gxLML.jpg" alt="Here a cool pic"></a></p>

<p>Can anybody help me? what I have to do and with what data I have to compare my prediction? Even a simple link could be helpful.</p>

<orbit><fundamental-astronomy><n-body-simulations>
How to discover Neptune from the Uranus orbit (by computer simulation)
Sergio Piccione
# How to discover Neptune from the Uranus orbit (by computer simulation)

(Asked by: username_0 on 2016-12-19)

I would like to demonstrate the existence of another planet (Neptune) by studying the discrepancy between the observation of the Uranus orbit and the mathematical prediction, this work was made from Le Verrier and I would like to understand his method.

I have read the Chapter 2, "The Discovery of Neptune (1845-1846)," in the biography Le Verrier -- Magnificent and Detestable Astronomer,
but it is too much in-depth and I did not understand very well his work.

I am studying the three-body problem (Sun, Uranus, Neptune) via Matlab and the two body problem (Sun, Uranus) taking the initial condition from here:

<http://nssdc.gsfc.nasa.gov/planetary/factsheet/uranusfact.html>

I have tried this method:
I put Uranus in the Perihelion with the Max. orbital velocity and I calculate the semi-major axis, and it is more accurate than the one obtained from putting Uranus and Neptune in the Perihelion with their respective Max. orbital velocity.

Here a cool pic made with Matlab:
[![Here a cool pic](https://i.sstatic.net/gxLML.jpg)](https://i.sstatic.net/gxLML.jpg)

Can anybody help me? what I have to do and with what data I have to compare my prediction? Even a simple link could be helpful.

**Tags:** orbit, fundamental-astronomy, n-body-simulations

---

## 1 answers

### Answer by username_1 on 2016-12-20. Score: 12 | Accepted answer

Here is what I did: 

Based on their masses, it is safest to initially consider Jupiter and Saturn as well as Uranus. It might also be fruitful to include the Earth in the analysis, to get relative positions, observation angles, etc. So, I will be considering: 


Sun
Earth
Jupiter
Saturn 
Uranus
Neptune

Get the standard gravitational parameters (μ) for all of them
Get initial positions and velocities via JPL/HORIZONS for all these planets. I had some data lying around from J2000.5, so I just used the state vectors from the 1st of January, 2000, at noon. 
Write an N-body integrator with built-in MATLAB tools. Integrate this incomplete Solar system once without Neptune, and once with Neptune included. 
Analyze and compare!

So, here's my data, and N-body integrator: 
function [t, yout_noNeptune, yout_withNeptune] = discover_Neptune()

    % Time of integration (in years)
    tspan = [0 97] * 365.25 * 86400;

    % std. gravitational parameters [km/s²/kg]
    mus_noNeptune = [1.32712439940e11; % Sun
                     398600.4415       % Earth
                     1.26686534e8      % Jupiter
                     3.7931187e7       % Saturn
                     5.793939e6];      % Uranus

    mus_withNeptune = [mus_noNeptune
                       6.836529e6]; % Neptune

    % Initial positions [km] and velocities [km/s] on 2000/Jan/1, 00:00
    % These positions describe the barycenter of the associated system,
    % e.g., sJupiter equals the statevector of the Jovian system barycenter.
    % Coordinates are expressed in ICRF, Solar system barycenter
    sSun     = [0 0 0 0 0 0].';
    sEarth   = [-2.519628815461580E+07  1.449304809540383E+08 -6.175201582312584E+02,...
                -2.984033716426881E+01 -5.204660244783900E+00  6.043671763866776E-05].';
    sJupiter = [ 5.989286428194381E+08  4.390950273441353E+08 -1.523283183395675E+07,...
                -7.900977458946710E+00  1.116263478937066E+01  1.306377465321731E-01].';
    sSaturn  = [ 9.587405702749230E+08  9.825345942920649E+08 -5.522129405702555E+07,...
                -7.429660072417541E+00  6.738335806405299E+00  1.781138895399632E-01].';
    sUranus  = [ 2.158728913593440E+09 -2.054869688179662E+09 -3.562250313222718E+07,...
                 4.637622471852293E+00  4.627114800383241E+00 -4.290473194118749E-02].';
    sNeptune = [ 2.514787652167830E+09 -3.738894534538290E+09  1.904284739289832E+07,...
                 4.466005624145428E+00  3.075618250100339E+00 -1.666451179600835E-01].';

    y0_noNeptune   = [sSun; sEarth; sJupiter; sSaturn; sUranus];
    y0_withNeptune = [y0_noNeptune; sNeptune];

    % Integrate the partial Solar system 
    % once with Neptune, and once without
    options = odeset('AbsTol', 1e-8,...
                     'RelTol', 1e-10);

    [t, yout_noNeptune]   = ode113(@(t,y) odefcn(t,y,mus_noNeptune)  , tspan, y0_noNeptune  , options);
    [~, yout_withNeptune] = ode113(@(t,y) odefcn(t,y,mus_withNeptune),     t, y0_withNeptune, options);

end

% The differential equation 
%
%    dy/dt = d/dt [r₀ v₀ r₁ v₁ r₂ v₂ ... rₙ vₙ]    
%          = [v₀ a₀ v₁ a₁ v₂ a₂ ... vₙ aₙ]    
%
%  with 
%
%    aₓ = Σₘ -G·mₘ/|rₘ-rₓ|² · (rₘ-rₓ) / |rₘ-rₓ| 
%       = Σₘ -μₘ·(rₘ-rₓ)/|rₘ-rₓ|³  
%
function dydt = odefcn(~, y, mus)

    % Split up position and velocity
    rs = y([1:6:end; 2:6:end; 3:6:end]);
    vs = y([4:6:end; 5:6:end; 6:6:end]);

     % Number of celestial bodies
    N = size(rs,2);

    % Compute interplanetary distances to the power -3/2
    df  = bsxfun(@minus, permute(rs, [1 3 2]), rs);
    D32 = permute(sum(df.^2), [3 2 1]).^(-3/2);
    D32(1:N+1:end) = 0; % (remove infs)

    % Compute all accelerations     
    as = -bsxfun(@times, mus.', D32);              % (magnitudes)    
    as = bsxfun(@times, df, permute(as, [3 2 1])); % (directions)    
    as = reshape(sum(as,2), [],1);                 % (total)

    % Output derivatives of the state vectors
    dydt = y;
    dydt([1:6:end; 2:6:end; 3:6:end]) = vs;
    dydt([4:6:end; 5:6:end; 6:6:end]) = as;

end

Here is the driver script I used to get some nice plots out: 
clc
close all

% Get coordinates from N-body simulation
[t, yout_noNeptune, yout_withNeptune] = discover_Neptune();

% For plot titles etc.
bodies = {'Sun'
          'Earth'
          'Jupiter'
          'Saturn'
          'Uranus'
          'Neptune'};


% Extract positions
rs_noNeptune   = yout_noNeptune  (:, [1:6:end; 2:6:end; 3:6:end]);
rs_withNeptune = yout_withNeptune(:, [1:6:end; 2:6:end; 3:6:end]);



% Figure of the whole Solar sysetm, just to check
% whether everything went OK
figure, clf, hold on
for ii = 1:numel(bodies)
    plot3(rs_withNeptune(:,3*(ii-1)+1),...
          rs_withNeptune(:,3*(ii-1)+2),...
          rs_withNeptune(:,3*(ii-1)+3),...
          'color', rand(1,3));
end

axis equal
legend(bodies);
xlabel('X [km]');
ylabel('Y [km]');
title('Just the Solar system, nothing to see here');


% Compare positions of Uranus with and without Neptune
rs_Uranus_noNeptune   = rs_noNeptune  (:, 13:15);
rs_Uranus_withNeptune = rs_withNeptune(:, 13:15);

figure, clf, hold on

plot3(rs_Uranus_noNeptune(:,1),...
      rs_Uranus_noNeptune(:,2),...
      rs_Uranus_noNeptune(:,3),...
      'b.');

plot3(rs_Uranus_withNeptune(:,1),...
      rs_Uranus_withNeptune(:,2),...
      rs_Uranus_withNeptune(:,3),...
      'r.');

axis equal
xlabel('X [km]');
ylabel('Y [km]');
legend('Uranus, no Neptune',...
       'Uranus, with Neptune');


% Norm of the difference over time
figure, clf, hold on

rescaled_t = t/365.25/86400;

dx = sqrt(sum((rs_Uranus_noNeptune - rs_Uranus_withNeptune).^2,2));
plot(rescaled_t,dx);
xlabel('Time [years]');
ylabel('Absolute offset [km]');
title({'Euclidian distance between'
       'the two Uranuses'});


% Angles from Earth
figure, clf, hold on

rs_Earth_noNeptune   = rs_noNeptune  (:, 4:6);
rs_Earth_withNeptune = rs_withNeptune(:, 4:6);

v0 = rs_Uranus_noNeptune   - rs_Earth_noNeptune;
v1 = rs_Uranus_withNeptune - rs_Earth_withNeptune;

nv0 = sqrt(sum(v0.^2,2));
nv1 = sqrt(sum(v1.^2,2));

dPhi = 180/pi * 3600 * acos(min(1,max(0, sum(v0.*v1,2) ./ (nv0.*nv1) )));
plot(rescaled_t, dPhi);

xlabel('Time [years]');
ylabel('Separation [arcsec]')
title({'Angular separation between the two'
       'Uranuses when observed from Earth'});

which I'll describe here step by step. 
First, a plot of the Solar system to check that the N-body integrator works as it should: 

Nice! Next, I wanted to see the difference between the positions of Uranus with and without the influence of Neptune. So, I extracted just the positions of those two Uranuses, and plotted them: 

...that's hardly useful. Even when zooming in greatly and rotating the heck out of it, this is just not a useful plot. So I looked at the evolution of the absolute Euclidian distance between the two Uranuses: 

That's starting to look more like it! Roughly 80 years after the start of our analysis, the two Uranuses are almost 6 million km apart! 
Large as that may sound, in the grander scale of things this might drown in the noise when we take measurements here on Earth. Plus, it still does not tell the whole story, as we'll see in a moment. So next, let's look at the angular difference between the observation vectors, from Earth towards the two Uranuses to see how large that angle is, and if it can stand out above the observational error thresholds: 

...whoa! Well over 300 arcseconds difference, plus all sorts of wobbly bobbley timey wimey rippling going on. That seems well within the observational capabilities of the time (although I can't find a reliable source on this so quickly; anyone?)
Just for good measure, I also produced that last plot leaving Jupiter and Saturn out of the picture. Although some perturbation theory had been developed in the 17th and 18th centuries, it was not very well developed and I doubt even Le Verrier took Jupiter into consideration (but again, I could be wrong; please correct me if you know more). 
So, here is the last plot without Jupiter and Saturn: 

Although there are differences, they are minute, and most importantly irrelevant for discovering Neptune.

> **Comments:**
> **username_2**: Brilliant answer!  
> **username_1**: @allthemikeysaretaken no, there's no source other than my code (which could also be wrong of course; I don't remember doing many validation tests on this to be honest). Those "squigglies" are the result of the Earth's orbit, and they _should_ therefore have a period of very close to 1 year. Eyeballing your first plot and my plot would indicate that your Earth is orbiting twice as slow as mine...Also, the picture you linked to uses different initial conditions (they put Uranus directly opposite Neptune), so we can't directly compare mine to theirs...  
> **username_1**: @allthemikeysaretaken Also: what exactly do you call "the" leap-frog integrator in Python? Verlet integration, though symplectic, may suffer from other sources of error that a variable-order multistep Adams-Bashforth-Moulton scheme (`ode113` in MATLAB) deals better with. It really depends on specifics whether one outperforms the other  

---

