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:


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

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.


2 Answers 2


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


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


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

close all

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

% For plot titles etc.
bodies = {'Sun'

% 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)
          'color', rand(1,3));

axis equal
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



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

the Solar system

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:

Two Uranuses, with and without Neptune

...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:

Time evolution of the 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:

Angular separation between the two Uranuses

...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:

Angular separation between the two Uranuses, leaving Jupiter and Saturn out of the equation

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

  • $\begingroup$ Brilliant answer! $\endgroup$
    – zephyr
    Commented Dec 20, 2016 at 20:36
  • $\begingroup$ Is there a source to confirm the accuracy of the last figure (separation in arcseconds of Uranus with and without Neptune, as viewed from respective Earth)? I ask because I tried to recreate this plot using the leap-frog integrator in python (and a time-step of 0.1 days); I got a similar shape but the scale of the angle is off. It may be a problem with my code, but the scale for the first 100 years in my simulation is comparable to the one found here (third image from the bottom). $\endgroup$
    – user33354
    Commented Sep 23, 2020 at 3:58
  • $\begingroup$ Viewing from the respective center-of-mass of the system (as opposed to viewing from the respective Earth frame) removes all the "squigglies" but the scale doesn't appear to change by doing this. $\endgroup$
    – user33354
    Commented Sep 23, 2020 at 4:02
  • 1
    $\begingroup$ @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... $\endgroup$ Commented Sep 24, 2020 at 23:47
  • $\begingroup$ @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 $\endgroup$ Commented Sep 24, 2020 at 23:48

If I understand correctly, you're modeling Uranus' orbit as an ellipse and want to compare it to Uranus' actual orbit as perturbed by Neptune? I don't have an answer, but Where can I find/visualize planets/stars/moons/etc positions? explains how to use SPICE, HORIZONS, and other tools to find Uranus' true position at a given time +-15000 years from now, including the best fit elliptical parameters (using HORIZONS "orbital elements" feature).

Of course, anything you do will be "circular" in some sense, since HORIZONS computed position of Uranus in the past already includes Neptune's perturbations.

If you could find tables of Uranus position predictions or something from the past, you might have something.

BTW, feel free to contact me (see profile for details) if this project extends beyond a stackexchange question.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .