Sensitivity and Robustness Analysis

Three experiments probing the robustness of the 4-parameter Neptune solver:

#	Question	Method
1	How badly does observational noise degrade the prediction?	Inject Gaussian noise at σ = 0, 1, 2, 5 arcsec; 10 trials each
2	How many years of data are needed?	Truncate at 1800, 1810, 1820, 1830, 1840
3	Does Bode’s Law bias hurt?	Vary the a prior: modern (25–35 AU), Bode (35–45 AU), agnostic (20–50 AU)

All experiments use 5-body synthetic observations as the benchmark.

Setup

Code

from IPython.display import display

from discoverneptune.plot_style import (
    apply_style,
    dual_render,
    fig_size,
    palette,
    truth_line,
)

apply_style()

import logging

import matplotlib.pyplot as plt
import numpy as np

from discoverneptune.plotting import (
    plot_bode_comparison,
    plot_noise_sensitivity,
    plot_residuals,
)
from discoverneptune.sensitivity import (
    run_bode_experiment,
    run_noise_experiment,
    run_truncation_experiment,
)
from discoverneptune.simulation import (
    MASS_NEPTUNE_APPROX,
    NeptuneCandidate,
    build_simulation_from_vectors,
    simulate_uranus_longitudes_from_vectors,
)
from discoverneptune.solver import solve
from discoverneptune.synthetic import build_synthetic_problem
from discoverneptune.validation import longitude_error_deg

logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s  %(levelname)-8s  %(name)s  %(message)s',
    datefmt='%H:%M:%S',
)
logger = logging.getLogger('sensitivity_nb')

%matplotlib inline

Build synthetic benchmark

5-body synthetic observations for Jupiter, Saturn, Uranus, and a known Neptune, spanning 1781-1846 (66 annual observations). These serve as the ground truth against which all experiments are evaluated.

Code

TRUE_NEPTUNE = NeptuneCandidate(
    mass_solar=MASS_NEPTUNE_APPROX,
    a=30.07,
    e=0.009,
    l_rad=np.radians(132.0),
)

logger.info("Building synthetic benchmark problem ...")
problem = build_synthetic_problem(TRUE_NEPTUNE)
jd_times    = problem["jd_times"]
truth_longs = problem["truth_longs"]
start_jd    = problem["start_jd"]
jupiter_sv  = problem["jupiter_sv"]
saturn_sv   = problem["saturn_sv"]
uranus_sv   = problem["uranus_sv"]

logger.info(
    "Loaded %d synthetic observations spanning %.1f years (JD %.1f - %.1f)",
    len(jd_times),
    (jd_times[-1] - jd_times[0]) / 365.25,
    jd_times[0],
    jd_times[-1],
)

Baseline optimizer run

A full-budget solve establishes the reference point: how well does the optimizer recover Neptune with clean synthetic data and the full 66-year baseline?

Code

NEPTUNE_TRUE_A = 30.07  # AU

logger.info('Running baseline solve ...')
baseline_result = solve(
    jd_times=jd_times,
    truth_longs=truth_longs,
    start_jd=start_jd,
    jupiter_sv=jupiter_sv,
    saturn_sv=saturn_sv,
    uranus_sv=uranus_sv,
    seed=42,
)

n = baseline_result.neptune
rms_arcsec = float(np.degrees(np.sqrt(baseline_result.sse / len(truth_longs))) * 3600.0)
print(f'Baseline:  a = {n.a:.3f} AU  (true: {NEPTUNE_TRUE_A} AU)')
print(f'           l = {np.degrees(n.l_rad):.2f} deg')
print(f'           RMS = {rms_arcsec:.2f} arcsec')
print(f'           converged = {baseline_result.converged}')

Baseline:  a = 30.070 AU  (true: 30.07 AU)
           l = 132.00 deg
           RMS = 0.00 arcsec
           converged = True

Observation noise tolerance

19th-century transit instruments achieved ~1 arcsec precision on a good night. How much noise can the solver tolerate before the prediction degrades? We inject Gaussian noise at σ = 0, 1, 2, and 5 arcsec and run 10 independent trials at each level.

Code

logger.info('Running noise experiment (4 levels x 10 trials) ...')
noise_trials = run_noise_experiment(
    jd_times, truth_longs, start_jd,
    jupiter_sv, saturn_sv, uranus_sv,
    sigma_levels=(0.0, 1.0, 2.0, 5.0),
    n_trials=10,
)
logger.info('Noise experiment done: %d trials', len(noise_trials))

def _build_noise(noise_trials, *, theme="light"):
    return plot_noise_sensitivity(noise_trials, theme=theme)

display(dual_render(_build_noise, noise_trials, alt="Neptune semi-major axis error vs. observation noise"))

Data sparsity: how many years of observation?

On this exact synthetic benchmark, RMS is not the informative quantity: once the solver has enough baseline to identify Neptune, several truncation windows drive the fit essentially to zero. The more interesting question is the 1846 localization problem: whether each truncated fit still localizes Neptune’s sky longitude well enough to point a telescope in the right part of the sky.

On the current regenerated surface, every displayed cutoff localizes the synthetic Neptune within the table’s shown precision. This result is specific to the closed-world benchmark and optimizer settings; it is not a claim about the information content of Le Verrier’s historical observations.

Code

logger.info('Running truncation experiment ...')
trunc_trials = run_truncation_experiment(
    jd_times, truth_longs, start_jd,
    jupiter_sv, saturn_sv, uranus_sv,
    cutoff_years=(1800, 1810, 1820, 1830, 1840),
)

def neptune_longitude_deg_at_final_epoch(neptune: NeptuneCandidate) -> float:
    sim = build_simulation_from_vectors(jupiter_sv, saturn_sv, uranus_sv, neptune)
    t_yr = (jd_times[-1] - start_jd) / 365.25
    sim.integrate(t_yr)
    sun = sim.particles[0]
    nep = sim.particles['Neptune']
    return float(np.degrees(np.arctan2(nep.y - sun.y, nep.x - sun.x)) % 360.0)

true_lon_1846_deg = neptune_longitude_deg_at_final_epoch(TRUE_NEPTUNE)
cutoff_years = []
lon_1846_deg = []
error_1846_deg = []
n_observations = []

for t in trunc_trials:
    fit_lon_deg = neptune_longitude_deg_at_final_epoch(t.result.neptune)
    err_deg = longitude_error_deg(np.radians(fit_lon_deg), np.radians(true_lon_1846_deg))
    cutoff_years.append(t.cutoff_year)
    lon_1846_deg.append(fit_lon_deg)
    error_1846_deg.append(err_deg)
    n_observations.append(t.n_observations)
    print(
        f'  cutoff={t.cutoff_year}  n={t.n_observations:3d}  '
        f'a={t.result.neptune.a:.2f} AU  '
        f'lambda_1846={fit_lon_deg:6.2f} deg  '
        f'error_1846={err_deg:5.2f} deg'
    )

def _build_truncation(cutoff_years, error_1846_deg, n_observations, *, theme="light"):
    apply_style(theme)
    pal = palette(theme)
    fig, ax = plt.subplots(figsize=fig_size("wide"))
    ax.plot(cutoff_years, error_1846_deg, 'o-', color=pal.after, linewidth=1.8)
    truth_line(ax, 5.0, label='5° band', theme=theme)
    ax.invert_xaxis()
    ax.set_xlabel('Last observation year included')
    ax.set_ylabel('1846 localization error (deg)')
    ax.set_title('How early can the truncated record still localize Neptune?')
    for year, err, n_obs in zip(cutoff_years, error_1846_deg, n_observations, strict=True):
        ax.annotate(f'n={n_obs}', (year, err), xytext=(0, 7), textcoords='offset points',
                    ha='center', fontsize=8)
    return fig

display(dual_render(_build_truncation, cutoff_years, error_1846_deg, n_observations,
            alt="1846 Neptune localization error vs. observation cutoff year"))

  cutoff=1800  n= 19  a=30.07 AU  lambda_1846=272.57 deg  error_1846= 0.00 deg
  cutoff=1810  n= 29  a=30.07 AU  lambda_1846=272.57 deg  error_1846= 0.00 deg
  cutoff=1820  n= 39  a=30.07 AU  lambda_1846=272.57 deg  error_1846= 0.00 deg
  cutoff=1830  n= 49  a=30.07 AU  lambda_1846=272.57 deg  error_1846= 0.00 deg
  cutoff=1840  n= 59  a=30.07 AU  lambda_1846=272.57 deg  error_1846= 0.00 deg

Bode’s Law prior bias

Bode’s Law predicted Neptune at ~38.8 AU (actual: 30.07 AU). If Le Verrier had restricted his search to a Bode-consistent range, would the solver still converge? We compare three bound choices: modern (25–35 AU), Bode (35–45 AU), and agnostic (20–50 AU).

Code

logger.info("Running Bode's Law prior experiment ...")
bode_trials = run_bode_experiment(
    jd_times, truth_longs, start_jd,
    jupiter_sv, saturn_sv, uranus_sv,
    priors=('modern', 'bode', 'agnostic'),
)
for t in bode_trials:
    print(f'  prior={t.prior_name:8s}  a in {t.a_bounds}  '
          f'found a={t.result.neptune.a:.2f} AU  '
          f'l={np.degrees(t.result.neptune.l_rad):.1f} deg')

def _build_bode(bode_trials, *, theme="light"):
    return plot_bode_comparison(bode_trials, theme=theme)

display(dual_render(_build_bode, bode_trials, alt="Effect of Bode's Law prior on Neptune recovery"))

  prior=modern    a in (25.0, 35.0)  found a=30.07 AU  l=132.0 deg
  prior=bode      a in (35.0, 45.0)  found a=35.00 AU  l=140.9 deg
  prior=agnostic  a in (20.0, 50.0)  found a=30.07 AU  l=132.0 deg

Residuals: before and after Neptune

A direct visualization of the problem Neptune solves. The main panel uses a shared y-axis so the scale contrast is visible at a glance; the zoomed inset shows the post-fit residual structure after the best-fit Neptune is added.

Code

neptune_absent = NeptuneCandidate(
    mass_solar=1e-12,   # negligible mass -- effectively no Neptune
    a=30.07,
    e=0.009,
    l_rad=baseline_result.neptune.l_rad,
)

sim_baseline = simulate_uranus_longitudes_from_vectors(
    jd_times, start_jd, jupiter_sv, saturn_sv, uranus_sv, neptune_absent
)
sim_best = simulate_uranus_longitudes_from_vectors(
    jd_times, start_jd, jupiter_sv, saturn_sv, uranus_sv, baseline_result.neptune
)

def _build_residuals(jd_times, truth_longs, sim_baseline, sim_best, *, theme="light"):
    return plot_residuals(jd_times, truth_longs, sim_baseline, sim_best, theme=theme)

display(dual_render(_build_residuals, jd_times, truth_longs, sim_baseline, sim_best,
            alt="Uranus residuals before and after Neptune"))

logger.info('Sensitivity analysis complete.')