Fix stack overflow in ND orbital mechanics problems

[ChrisRackauckas-Claude]

Summary

Fixes the stack overflow errors in ND1-ND5 problems by simplifying the symbolic expression for gravitational force calculation. This enables all Kepler orbital mechanics benchmarks to work correctly.

Problem

The ND problems (orbital mechanics/Kepler orbit benchmarks) were causing stack overflow errors during ModelingToolkit compilation. The root cause was the complex symbolic expression:

r = sqrt(y[1]^2 + y[2]^2)^3  # Complex: (x² + y²)^(3/2)

When used in the gravitational force equations (-y[1]) / r and (-y[2]) / r, this created overly complex symbolic manipulations that caused ModelingToolkit's symbolic engine to exceed stack limits during compilation.

Solution

Replaced the complex expression with mathematically equivalent intermediate variables:

# Before (causes stack overflow):
r = sqrt(y[1]^2 + y[2]^2)^3

# After (compiles cleanly):
r_squared = y[1]^2 + y[2]^2
r_cubed = r_squared * sqrt(r_squared)  

This represents the same physics (gravitational force law F ∝ 1/r³) but breaks down the symbolic computation into manageable steps for ModelingToolkit.

Comprehensive Local Testing ✅

Verified that stack overflow is completely resolved:

1. System Compilation

• ✅ ND1 system compiles without stack overflow
• ✅ All ND1-ND5 problems created successfully
• ✅ No compilation timeouts or excessive symbolic processing

2. Problem Solving

• ✅ All 5 eccentricity variants solve successfully:
  • ND1 (e=0.1): Nearly circular orbit
  • ND2 (e=0.3): Mild elliptical orbit
  • ND3 (e=0.5): Moderate elliptical orbit
  • ND4 (e=0.7): High eccentricity orbit
  • ND5 (e=0.9): Very elongated comet-like orbit

3. Multiple Solver Algorithms

• ✅ Tsit5, Vern6, Vern7, Vern9 all work correctly
• ✅ Consistent results across all algorithms

4. High-Accuracy Solutions

• ✅ Reference solutions with abstol=1e-14, reltol=1e-14 generated successfully
• ✅ No compilation issues with extreme precision requirements

5. WorkPrecisionSet Generation

• ✅ Original failing scenario now works
• ✅ WorkPrecisionSet created with multiple tolerance combinations
• ✅ Benchmark automation fully functional

6. Physics Verification

• ✅ Realistic orbital mechanics results for all eccentricities
• ✅ Energy conservation properties maintained
• ✅ Proper orbital dynamics (higher eccentricity → more elongated orbits)

Mathematical Verification

Expressions are mathematically identical:

• sqrt(x² + y²)³ = 125.0
• (x² + y²) * sqrt(x² + y²) = 125.0
• Difference: 0.0 ✅

Impact

This critical fix enables:

• Complete Enright-Pryce ND benchmark suite (5 orbital problems)
• Celestial mechanics test cases for ODE solver evaluation
• Benchmark automation that was previously broken by stack overflow
• Physics-accurate orbital simulations from nearly circular to highly elliptical

The ND problems are essential benchmarks for testing adaptive time-stepping algorithms on challenging orbital mechanics scenarios.

🤖 Generated with Claude Code