Updated docs on benchmarks

benches/BS3_VS_DP5_RESULTS.md

@@ -1,6 +1,6 @@
 # BS3 vs DP5 Benchmark Results
 
-Generated: 2025-10-23
+Generated: 2026-09-29 (re-measured; first run 2025-10-23)
 
 ## Summary
 
@@ -10,7 +10,7 @@ Comprehensive performance comparison between **BS3** (Bogacki-Shampine 3rd order
 
 ### Overall Performance Comparison
 
-**DP5 is consistently faster than BS3 across all tested scenarios**, typically by a factor of **1.5x to 4.3x**.
+**DP5 is consistently faster than BS3 across all tested scenarios**, typically by a factor of **1.6x to 4.8x**.
 
 This might seem counterintuitive since BS3 uses fewer stages (4 vs 7), but several factors explain DP5's superior performance:
 
@@ -32,8 +32,8 @@ Despite being slower in these benchmarks, BS3 still has value:
 
 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 3.28 µs | 1.92x slower |
-| **DP5** | 1.70 µs | baseline |
+| **BS3** | 2.96 µs | 2.08x slower |
+| **DP5** | 1.43 µs | baseline |
 
 Simple 1D problem with smooth exponential solution.
 
@@ -41,8 +41,8 @@ Simple 1D problem with smooth exponential solution.
 
 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 30.70 µs | 2.25x slower |
-| **DP5** | 13.67 µs | baseline |
+| **BS3** | 27.45 µs | 2.26x slower |
+| **DP5** | 12.16 µs | baseline |
 
 2D conservative system with periodic solution.
 
@@ -50,8 +50,8 @@ Simple 1D problem with smooth exponential solution.
 
 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 132.35 µs | 3.57x slower |
-| **DP5** | 37.11 µs | baseline |
+| **BS3** | 113.66 µs | 3.65x slower |
+| **DP5** | 31.15 µs | baseline |
 
 Nonlinear 2D system with trigonometric terms.
 
@@ -59,8 +59,8 @@ Nonlinear 2D system with trigonometric terms.
 
 | Method | Time | Ratio |
 |--------|------|-------|
-| **BS3** | 124.72 µs | 1.45x slower |
-| **DP5** | 86.10 µs | baseline |
+| **BS3** | 102.04 µs | 1.56x slower |
+| **DP5** | 65.52 µs | baseline |
 
 Higher-dimensional problem with gravitational dynamics.
 
@@ -68,8 +68,8 @@ Higher-dimensional problem with gravitational dynamics.
 
 | Method | Time (solve + 100 interpolations) | Ratio |
 |--------|-----------------------------------|-------|
-| **BS3** | 19.68 µs | 4.81x slower |
-| **DP5** | 4.09 µs | baseline |
+| **BS3** | 16.28 µs | 4.70x slower |
+| **DP5** | 3.46 µs | baseline |
 
 BS3 uses cubic Hermite interpolation, DP5 uses optimized 5th order interpolation.
 
@@ -79,14 +79,14 @@ Performance across different tolerance levels (`y' = -y` problem):
 
 | Tolerance | BS3 Time | DP5 Time | Ratio (BS3/DP5) |
 |-----------|----------|----------|-----------------|
-| 1e-3 | 1.63 µs | 1.26 µs | 1.30x |
-| 1e-4 | 2.61 µs | 1.54 µs | 1.70x |
-| 1e-5 | 4.64 µs | 2.03 µs | 2.28x |
-| 1e-6 | 8.76 µs | ~2.6 µs* | ~3.4x* |
-| 1e-7 | -** | -** | - |
+| 1e-3 | 1.15 µs | 0.92 µs | 1.25x |
+| 1e-4 | 2.01 µs | 1.19 µs | 1.69x |
+| 1e-5 | 3.41 µs | 1.45 µs | 2.36x |
+| 1e-6 | 6.72 µs | 1.92 µs | 3.50x |
+| 1e-7 | 13.55 µs | 2.82 µs | 4.80x |
 
-\* Estimated from trend (benchmark timed out)
-\** Not completed
+Both 1e-6 and 1e-7 are now measured. The earlier run estimated 1e-6 from the
+trend and did not complete 1e-7.
 
 **Observation**: The performance gap widens as tolerance tightens, because DP5's higher order allows it to take larger steps while maintaining accuracy.
 
@@ -96,15 +96,17 @@ Performance across different tolerance levels (`y' = -y` problem):
 
 1. **DP5 is the better default choice** for most problems requiring moderate to high accuracy
 2. **Performance gap increases** with tighter tolerances (favoring DP5)
-3. **Higher dimensions** slightly favor BS3 relative to DP5 (1.45x vs 3.57x slowdown)
-4. **Interpolation** strongly favors DP5 (4.8x faster)
+3. **Higher dimensions** slightly favor BS3 relative to DP5 (1.56x vs 3.65x slowdown)
+4. **Interpolation** strongly favors DP5 (4.7x faster)
 
 ### Implementation Quality
 
 Both integrators pass all accuracy and convergence tests:
 - ✅ BS3: 3rd order convergence rate verified
 - ✅ DP5: 5th order convergence rate verified (existing implementation)
-- ✅ Both: FSAL property correctly implemented
+- ⚠️ Both: FSAL is a property of the *methods*, but neither implementation
+  exploits it - each step recomputes `k[0] = f(t, y)` rather than reusing the
+  previous step's last stage. See roadmap feature 36.
 - ✅ Both: Dense output accurate to specified order
 
 ### Future Optimizations