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, palette, category_cycle, dual_render,
    SEQUENTIAL_CMAP, truth_line, annotate_stat, fig_size,
)
apply_style()

import logging
import numpy as np
import matplotlib.pyplot as plt

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

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.

By 1800, the synthetic benchmark already contains enough information to localize Neptune. The small 1810 deviation is a negligible optimizer-level wobble, not evidence that additional observations are harmful.

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)
    ax.legend()
    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=273.01 deg  error_1846= 0.00 deg
  cutoff=1810  n= 29  a=30.01 AU  lambda_1846=273.17 deg  error_1846= 0.16 deg
  cutoff=1820  n= 39  a=30.07 AU  lambda_1846=273.01 deg  error_1846= 0.00 deg
  cutoff=1830  n= 49  a=30.07 AU  lambda_1846=273.01 deg  error_1846= 0.00 deg
  cutoff=1840  n= 59  a=30.07 AU  lambda_1846=273.01 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=141.0 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.')