Conservation-Law Correction and Adaptive-Precision Strategies for Long-Term ODE Integration: A Proof of Concept on the Sitnikov Problem

# Conservation-Law Correction and Adaptive-Precision Strategies for Long-Term ODE Integration: A Proof of Concept on the Sitnikov Problem --- **Repository**: FCE Engine Proof — Sitnikov Problem**Version**: 1.0**Date**: 2026**Keywords**: numerical integration, error correction, three-body problem, Poincaré section, energy conservation, adaptive precision, DOP853, Sitnikov problem --- ## Abstract Standard numerical integrators for ordinary differential equations accumulate error over long integration times, a fundamental limitation for orbital mechanics and other Hamiltonian systems. We present a post-processing correction framework — the Fractal Correction Engine (FCE) — that monitors an integrator's trajectory at Poincaré section crossings and intervenes to reduce accumulated error. We demonstrate the framework on the Sitnikov restricted three-body problem using two distinct strategies matched to the dynamics: **Strategy 1 (integrable regime, $e = 0$):** A conservation-law correction that exploits the exact constancy of the crossing velocity at $z = 0$. Any measured deviation from this constant is identified as pure integrator energy drift and corrected directly. This achieves a **437-fold reduction** in energy drift using a deliberately coarse DOP853 integrator ($\text{rtol} = 10^{-2}$). **Strategy 2 (chaotic regime, $e = 0.3$):** An adaptive-precision re-integration that replaces the coarse integrator with orbit-by-orbit integration at moderate precision ($\text{rtol} = 10^{-8}$), using Poincaré section crossings as natural restart points. This achieves a **46,248-fold reduction** in position error and **31,561-fold reduction** in energy error relative to a bare coarse integrator ($\text{rtol} = 10^{-1}$), while remaining **3.5 times cheaper** than a continuous ultra-fine reference integration ($\text{rtol} = 10^{-14}$). We are transparent that these are fundamentally different mechanisms: Strategy 1 is a genuine physics-based correction applied to a cheap integrator, while Strategy 2 is adaptive precision guided by orbital structure. Both are valuable, but they make different claims about what the FCE provides. --- ## 1. Introduction ### 1.1 The Problem: Long-Term Integrator Drift Numerical integration of Hamiltonian systems over long time spans is a central challenge in celestial mechanics, molecular dynamics, and dynamical astronomy. Even high-order integrators such as DOP853 (Dormand--Prince 8(5,3), an explicit Runge--Kutta method of order 8) accumulate truncation error that manifests as secular drift in conserved quantities. For an energy-conserving system integrated with a non-symplectic method at relative tolerance $\text{rtol}$, the energy error typically grows as: $$\Delta E(t) \sim \mathcal{O}(\text{rtol}) \cdot t$$ over long integration times $t$. This linear drift means that even small per-step errors compound into significant trajectory deviation. Symplectic integrators address this by preserving the symplectic structure of Hamiltonian phase space, bounding energy error to oscillations of size $\mathcal{O}(h^p)$ where $h$ is the step size and $p$ is the integrator order. However, symplectic integrators require the problem to be expressed in Hamiltonian form, use fixed step sizes (adaptive stepping breaks symplecticity), and are typically limited to low order ($p \leq 4$ in practice). ### 1.2 The FCE Concept The Fractal Correction Engine (FCE) takes a different approach: rather than modifying the integrator, it **monitors** the integrator's output at geometrically significant points and **corrects** detected errors. The key idea is: 1. **Monitor** the trajectory at Poincaré section crossings (geometrically natural reference points in phase space).2. **Detect** integrator error by comparing observed quantities against known invariants or high-precision recomputation.3. **Correct** the state at each crossing and restart the integrator from the corrected state. This framework is integrator-agnostic: it wraps any black-box ODE solver and improves its long-term accuracy. The cost of correction depends on the strategy used, ranging from negligible (conservation-law lookup) to moderate (local re-integration at higher precision). ### 1.3 Two Regimes, Two Strategies A central finding of this work is that the correction strategy must be matched to the dynamics: - **Integrable systems** possess exact conservation laws that provide a "free" error signal. The FCE can correct a cheap integrator toward the known invariant at essentially zero cost. This is a genuine correction — the cheap integrator does the work, and the FCE fixes the drift. - **Chaotic systems** lack persistent patterns that can be exploited for correction. Pattern-fitting approaches (Fourier decomposition, phase-based prediction) were tested and found to fail because the dynamics are genuinely irregular. The effective strategy is orbit-by-orbit re-integration at moderate precision, using Poincaré crossings as natural restart boundaries. This is more honestly described as *adaptive-precision integration guided by orbital structure* rather than a "correction" applied to a coarse integrator. We present both strategies transparently, with honest assessment of what each accomplishes and how. ### 1.4 Related Work The idea of correcting a numerical solution by estimating its local defect has a long history. Zadunaisky [6] introduced the *defect correction* principle: integrate the ODE at moderate precision, construct an interpolant of the numerical solution, substitute it back into the ODE to estimate the residual, and solve a correction equation for the error. This yields a higher-order solution at modest additional cost. Our Strategy 2 is conceptually adjacent — we re-integrate each orbit segment at higher precision rather than estimating the defect analytically, but the philosophy of local error estimation and correction is shared. The broader framework of *geometric numerical integration* [4] emphasizes preserving structural properties of the continuous system (symplecticity, energy conservation, time-reversal symmetry) in the discrete approximation. Symplectic integrators, the most successful family in this tradition, achieve bounded energy error over exponentially long times for Hamiltonian systems. Our Strategy 1 pursues a complementary goal: rather than building conservation into the integrator, we enforce it externally at Poincaré section crossings, allowing the use of any standard (non-symplectic) integrator. Poincaré sections have been used extensively in celestial mechanics for qualitative analysis of dynamics [5, 7], but less commonly as intervention points for numerical correction. The use of surface-of-section crossings as natural restart boundaries for orbit-by-orbit integration is related to the *Poincaré map* approach to studying stability, where the continuous flow is reduced to a discrete map at each section crossing. Our contribution is using this structure operationally — as a schedule for error correction or precision adaptation — rather than purely for analysis. --- ## 2. The Sitnikov Problem ### 2.1 Physical Setup The Sitnikov problem is a restricted three-body problem in which two equal-mass primaries ($m_1 = m_2 = 1/2$) orbit their common barycenter in the $xy$-plane, while a massless test particle moves along the perpendicular $z$-axis. The test particle feels the gravitational attraction of both primaries but does not affect their motion. This problem is a standard testbed for numerical methods in celestial mechanics because:- It reduces the three-body problem to a single degree of freedom ($z$).- It exhibits both integrable ($e = 0$) and chaotic ($e > 0$) dynamics depending on the eccentricity of the primary orbit.- The equation of motion is simple but captures genuine gravitational physics. ### 2.2 Equations of Motion The primaries orbit their barycenter with eccentricity $e$ and period $T = 2\pi$ (in normalized units). Their distance from the barycenter is: $$r(t) = \frac{1}{2}\bigl(1 - e\cos E(t)\bigr)$$ where the eccentric anomaly $E(t)$ satisfies Kepler's equation: $$E - e\sin E = M, \qquad M = t \bmod 2\pi$$ For circular orbits ($e = 0$), this simplifies to $r = 1/2$ (constant). The test particle's equation of motion along the $z$-axis is: $$\ddot{z} = -\frac{z}{\bigl(r(t)^2 + z^2\bigr)^{3/2}}$$ This is the $z$-component of the gravitational acceleration from both primaries, each located at distance $\sqrt{r^2 + z^2}$ from the test particle. ### 2.3 Energy The test particle's specific energy is: $$\mathcal{E}(z, \dot{z}, t) = \frac{1}{2}\dot{z}^2 - \frac{1}{\sqrt{r(t)^2 + z^2}}$$ The time derivative of this energy is: $$\frac{d\mathcal{E}}{dt} = \frac{r(t)\,\dot{r}(t)}{\bigl(r(t)^2 + z^2\bigr)^{3/2}}$$ **For $e = 0$:** Since $r(t) = 1/2$ is constant, $\dot{r} = 0$, giving $d\mathcal{E}/dt = 0$. The energy is an **exact invariant**. Any observed energy change in a numerical solution is purely integrator error. **For $e > 0$:** The primaries' distance varies in time ($\dot{r} \neq 0$), so the test particle energy genuinely oscillates. The energy is *not* conserved — it exchanges energy with the primary orbit through the time-dependent gravitational field. ### 2.4 Dynamics For $e = 0$, the Sitnikov problem is integrable. All bounded orbits are periodic, and the motion in the $(z, \dot{z})$ phase plane traces closed curves. For $e > 0$, the problem exhibits a transition to chaos. Orbits become irregular, with sensitive dependence on initial conditions. The Lyapunov exponent is positive for generic initial conditions at moderate eccentricity. This makes long-term trajectory prediction fundamentally difficult — nearby trajectories diverge exponentially. --- ## 3. Method: Poincaré Section Monitoring ### 3.1 The Poincaré Section We define the Poincaré section as the plane $z = 0$ in the $(z