diff --git a/Cargo.toml b/Cargo.toml index 84d1fe2..afa9433 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,16 +9,20 @@ repository = "https://git.rcjohnstone.com/connor/ordinary-diffeq" license = "MIT" exclude = ["roadmap/"] -# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html +[features] +default = ["autodiff"] +autodiff = ["dep:num-dual"] [dependencies] serde = { version = "1.0", features = ["derive"] } -nalgebra = { version = "0.34", features = ["serde-serialize"] } +nalgebra = { version = "0.35", features = ["serde-serialize"] } num-traits = "0.2.19" +num-dual = { version = "0.15", optional = true } roots = "0.0.8" [dev-dependencies] approx = "0.5" +num-dual = "0.15" criterion = "0.7.0" [[bench]] diff --git a/VERN7_BENCHMARK_REPORT.md b/VERN7_BENCHMARK_REPORT.md index 69b2811..ad26b98 100644 --- a/VERN7_BENCHMARK_REPORT.md +++ b/VERN7_BENCHMARK_REPORT.md @@ -1,17 +1,17 @@ # Vern7 Performance Benchmark Report -**Date**: 2025-10-24 -**Test System**: Linux 6.17.4-arch2-1 +**Date**: 2026-09-29 +**Test System**: Linux 7.2.6-arch2-1 **Optimization Level**: Release build with full optimizations ## Executive Summary Vern7 demonstrates **substantial performance advantages** over lower-order methods (BS3 and DP5) at tight tolerances (1e-8 to 1e-12), achieving: -- **2.7x faster** than DP5 at 1e-10 tolerance (exponential problem) -- **3.8x faster** than DP5 in harmonic oscillator -- **8.8x faster** than DP5 for orbital mechanics -- **51x faster** than BS3 in harmonic oscillator -- **1.65x faster** than DP5 for interpolation workloads +- **2.9x faster** than DP5 at 1e-10 tolerance (exponential problem) +- **4.4x faster** than DP5 in harmonic oscillator +- **1.1x faster** than DP5 for orbital mechanics +- **83x faster** than BS3 in harmonic oscillator +- **2.3x faster** than DP5 for interpolation workloads These results confirm Vern7's design goal: **maximum efficiency for high-accuracy requirements**. @@ -23,14 +23,14 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac | Method | Time (μs) | Relative Speed | Speedup vs BS3 | |--------|-----------|----------------|----------------| -| **Vern7** | **3.81** | **1.00x** (baseline) | **51.8x** | -| DP5 | 10.43 | 2.74x slower | 18.9x | -| BS3 | 197.37 | 51.8x slower | 1.0x | +| **Vern7** | **2.96** | **1.00x** (baseline) | **53.2x** | +| DP5 | 8.46 | 2.86x slower | 18.6x | +| BS3 | 157.35 | 53.2x slower | 1.0x | **Analysis**: -- Vern7 is **2.7x faster** than DP5 and **51x faster** than BS3 +- Vern7 is **2.9x faster** than DP5 and **53x faster** than BS3 - BS3's 3rd-order method requires many tiny steps to maintain 1e-10 accuracy -- DP5's 5th-order is better but still requires ~2.7x more work than Vern7 +- DP5's 5th-order is better but still requires ~2.9x more work than Vern7 - Vern7's 7th-order allows much larger step sizes while maintaining accuracy --- @@ -41,15 +41,15 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac | Method | Time (μs) | Relative Speed | Speedup vs BS3 | |--------|-----------|----------------|----------------| -| **Vern7** | **26.89** | **1.00x** (baseline) | **55.1x** | -| DP5 | 102.74 | 3.82x slower | 14.4x | -| BS3 | 1,481.4 | 55.1x slower | 1.0x | +| **Vern7** | **19.33** | **1.00x** (baseline) | **83.2x** | +| DP5 | 84.19 | 4.36x slower | 19.1x | +| BS3 | 1,608.7 | 83.2x slower | 1.0x | **Analysis**: -- Vern7 is **3.8x faster** than DP5 and **55x faster** than BS3 +- Vern7 is **4.4x faster** than DP5 and **83x faster** than BS3 - Smooth periodic problems like harmonic oscillators are ideal for high-order methods -- BS3 requires ~1.5ms due to tiny steps needed for tight tolerance -- DP5 needs ~103μs, still significantly more than Vern7's 27μs +- BS3 requires ~1.6ms due to tiny steps needed for tight tolerance +- DP5 needs ~84μs, still significantly more than Vern7's 19μs - Higher dimensionality (2D vs 1D) amplifies the advantage of larger steps --- @@ -60,15 +60,18 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac | Method | Time (μs) | Relative Speed | Speedup | |--------|-----------|----------------|---------| -| **Vern7** | **98.75** | **1.00x** (baseline) | **8.77x** | -| DP5 | 865.79 | 8.77x slower | 1.0x | +| **Vern7** | **83.88** | **1.00x** (baseline) | **1.13x** | +| DP5 | 94.68 | 1.13x slower | 1.0x | **Analysis**: -- Vern7 is **8.8x faster** than DP5 for this challenging 6D problem +- Vern7 is **1.1x faster** than DP5 on this problem +- Step counts: Vern7 184, DP5 440. Vern7 takes 2.4x fewer steps, but each step + costs 16 stages against DP5's 7, so the advantage very nearly cancels - Orbital mechanics requires tight tolerances to maintain energy conservation - BS3 was too slow to include in the benchmark at this tolerance -- 6D problem with long integration time shows Vern7's scalability -- This represents realistic astrodynamics/orbital mechanics workloads +- Vern7's advantage here is in step count, which is what matters once the + right-hand side is expensive - a real gravity model rather than a two-body + point mass --- @@ -78,17 +81,20 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac | Method | Time (μs) | Relative Speed | Notes | |--------|-----------|----------------|-------| -| **Vern7** | **11.05** | **1.00x** (baseline) | Lazy extra stages | -| DP5 | 18.27 | 1.65x slower | Standard dense output | +| **Vern7** | **7.28** | **1.00x** (baseline) | Lazy extra stages | +| DP5 | 16.73 | 2.30x slower | Standard dense output | **Analysis**: -- Vern7 with lazy computation is **1.65x faster** than DP5 +- Vern7 with lazy computation is **2.3x faster** than DP5 - First interpolation triggers lazy computation of 6 extra stages (k11-k16) - Subsequent interpolations reuse cached extra stages (~10ns RefCell overhead) - Despite computing extra stages, Vern7 is still faster overall due to: 1. Fewer total integration steps (larger step sizes) 2. Higher accuracy interpolation (7th order vs 5th order) -- Lazy computation adds minimal overhead (~6μs for 6 stages, amortized over 100 interpolations) +- Lazy computation adds minimal overhead, amortized over 100 interpolations +- Each benchmark iteration builds a fresh `Solution`, so the lazy stages are + computed once per iteration rather than served from the `RefCell` cache - the + cost of interpolating into a solution you have just produced --- @@ -100,11 +106,11 @@ These results confirm Vern7's design goal: **maximum efficiency for high-accurac | Tolerance | DP5 (μs) | Vern7 (μs) | Speedup | Winner | |-----------|----------|------------|---------|--------| -| 1e-6 | 2.63 | 2.05 | 1.28x | Vern7 | -| 1e-7 | 3.71 | 2.74 | 1.35x | Vern7 | -| 1e-8 | 5.43 | 3.12 | 1.74x | Vern7 | -| 1e-9 | 7.97 | 3.86 | 2.06x | **Vern7** | -| 1e-10 | 11.33 | 5.33 | 2.13x | **Vern7** | +| 1e-6 | 1.94 | 1.51 | 1.29x | Vern7 | +| 1e-7 | 2.82 | 1.93 | 1.46x | Vern7 | +| 1e-8 | 4.20 | 2.29 | 1.83x | Vern7 | +| 1e-9 | 6.25 | 2.97 | 2.11x | **Vern7** | +| 1e-10 | 9.52 | 3.89 | 2.45x | **Vern7** | ### Performance Scaling Chart (Conceptual) @@ -128,7 +134,7 @@ Time (μs) **Analysis**: - At **moderate tolerances (1e-6)**: Vern7 is 1.3x faster -- At **tight tolerances (1e-10)**: Vern7 is 2.1x faster +- At **tight tolerances (1e-10)**: Vern7 is 2.5x faster - **Crossover point**: Vern7 becomes increasingly advantageous as tolerance tightens - DP5's time scales roughly quadratically with tolerance - Vern7's time scales more slowly (higher order = larger steps) @@ -157,22 +163,23 @@ Time (μs) ### Lazy Computation Impact The lazy computation of extra stages (k11-k16) provides: -- **Minimal overhead**: ~6μs to compute 6 extra stages +- **Minimal overhead**: the 6 extra stages are computed once per step that is interpolated into - **Cache efficiency**: Extra stages computed once per interval, reused for multiple interpolations - **Memory efficiency**: Only computed when interpolation is requested -- **Performance**: Despite extra computation, still 1.65x faster than DP5 for interpolation workloads +- **Performance**: Despite extra computation, still 2.3x faster than DP5 for interpolation workloads ### Step Size Comparison -Estimated step sizes at 1e-10 tolerance for exponential problem: +Measured at 1e-10 tolerance on the exponential problem (`y' = y`, t = 0 to 4): | Method | Avg Step Size | Steps Required | Function Evals | |--------|---------------|----------------|----------------| -| BS3 | ~0.002 | ~2000 | ~8000 | -| DP5 | ~0.01 | ~400 | ~2400 | -| **Vern7** | ~0.05 | **~80** | **~800** | +| BS3 | 0.0014 | 2,916 | 11,664 | +| DP5 | 0.0354 | 113 | 791 | +| **Vern7** | 0.1379 | **29** | **290** | -**Vern7 requires ~3x fewer function evaluations than DP5.** +**Vern7 requires 2.7x fewer function evaluations than DP5**, taking 3.9x fewer +steps at 10 stages each against DP5's 7. --- @@ -185,7 +192,7 @@ Our Rust implementation achieves performance comparable to Julia's highly-optimi | Step computation | Highly optimized, FSAL | Optimized, no FSAL | | Lazy interpolation | ✓ | ✓ | | Stage caching | RefCell-based | RefCell-based (~10ns) | -| Memory allocation | Minimal | Minimal | +| Memory allocation | Cached in the integrator | Two heap allocations per step | | Relative speed | Baseline | ~Comparable | **Note**: Direct comparison difficult due to different hardware and problems, but algorithmic approach is identical. @@ -213,8 +220,8 @@ Our Rust implementation achieves performance comparable to Julia's highly-optimi Vern7 successfully achieves its design goal of being the **most efficient method for high-accuracy non-stiff problems**. The implementation with lazy computation of extra stages provides: -- ✅ **2-9x speedup** over DP5 at tight tolerances -- ✅ **50x+ speedup** over BS3 at tight tolerances +- ✅ **1.1-4.4x speedup** over DP5 at tight tolerances +- ✅ **53-83x speedup** over BS3 at tight tolerances - ✅ **Efficient lazy interpolation** with minimal overhead - ✅ **Full 7th-order accuracy** for both steps and interpolation - ✅ **Memory-efficient caching** with RefCell @@ -225,7 +232,7 @@ The results validate the effort invested in implementing the complex 16-stage in ## Appendix: Benchmark Configuration -**Hardware**: Not specified (Linux system) +**Hardware**: Linux 7.2.6-arch2-1 **Compiler**: rustc (release mode, full optimizations) **Measurement Tool**: Criterion.rs v0.7.0 **Sample Size**: 100 samples per benchmark diff --git a/benches/BS3_VS_DP5_RESULTS.md b/benches/BS3_VS_DP5_RESULTS.md index 57fb055..940b37a 100644 --- a/benches/BS3_VS_DP5_RESULTS.md +++ b/benches/BS3_VS_DP5_RESULTS.md @@ -1,6 +1,6 @@ # BS3 vs DP5 Benchmark Results -Generated: 2025-10-23 +Generated: 2026-09-29 ## 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.2x to 5.1x**. 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.99 µs | 2.21x slower | +| **DP5** | 1.35 µ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.37 µs | 2.76x slower | +| **DP5** | 9.92 µ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.44 µs | 3.81x slower | +| **DP5** | 29.74 µ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** | 98.02 µs | 5.07x slower | +| **DP5** | 19.32 µ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.34 µs | 4.67x slower | +| **DP5** | 3.50 µs | baseline | BS3 uses cubic Hermite interpolation, DP5 uses optimized 5th order interpolation. @@ -79,14 +79,11 @@ 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 | -** | -** | - | - -\* Estimated from trend (benchmark timed out) -\** Not completed +| 1e-3 | 1.10 µs | 0.91 µs | 1.21x | +| 1e-4 | 1.87 µs | 1.19 µs | 1.57x | +| 1e-5 | 3.46 µs | 1.44 µs | 2.40x | +| 1e-6 | 6.77 µs | 1.95 µs | 3.46x | +| 1e-7 | 13.73 µs | 2.79 µs | 4.91x | **Observation**: The performance gap widens as tolerance tightens, because DP5's higher order allows it to take larger steps while maintaining accuracy. @@ -96,15 +93,18 @@ 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 favor DP5 most strongly** (5.07x on the 6D orbit against + 2.21x on 1D exponential decay) +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 diff --git a/benches/bs3_vs_dp5.rs b/benches/bs3_vs_dp5.rs index b59e5cc..ae75310 100644 --- a/benches/bs3_vs_dp5.rs +++ b/benches/bs3_vs_dp5.rs @@ -27,9 +27,8 @@ fn bench_exponential_decay(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, params); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, bs3, controller).solve(); - }); + let mut problem = Problem::new(ode, bs3, controller); + black_box(problem.solve()); }); }); @@ -37,9 +36,8 @@ fn bench_exponential_decay(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, params); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -66,9 +64,8 @@ fn bench_harmonic_oscillator(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 20.0, y0, ()); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, bs3, controller).solve(); - }); + let mut problem = Problem::new(ode, bs3, controller); + black_box(problem.solve()); }); }); @@ -76,9 +73,8 @@ fn bench_harmonic_oscillator(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 20.0, y0, ()); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -109,9 +105,8 @@ fn bench_pendulum(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, params); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, bs3, controller).solve(); - }); + let mut problem = Problem::new(ode, bs3, controller); + black_box(problem.solve()); }); }); @@ -119,9 +114,8 @@ fn bench_pendulum(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, params); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -160,9 +154,8 @@ fn bench_orbit_6d(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, bs3, controller).solve(); - }); + let mut problem = Problem::new(ode, bs3, controller); + black_box(problem.solve()); }); }); @@ -170,9 +163,8 @@ fn bench_orbit_6d(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -199,12 +191,10 @@ fn bench_interpolation(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 5.0, y0, ()); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - let mut problem = Problem::new(ode, bs3, controller); - let solution = problem.solve(); - // Interpolate at 100 points - let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect(); - }); + let mut problem = Problem::new(ode, bs3, controller); + let solution = problem.solve(); + // Interpolate at 100 points + black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::>()); }); }); @@ -213,12 +203,10 @@ fn bench_interpolation(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 5.0, y0, ()); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - let mut problem = Problem::new(ode, dp45, controller); - let solution = problem.solve(); - // Interpolate at 100 points - let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect(); - }); + let mut problem = Problem::new(ode, dp45, controller); + let solution = problem.solve(); + // Interpolate at 100 points + black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::>()); }); }); @@ -245,9 +233,8 @@ fn bench_tolerance_scaling(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, ()); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, bs3, controller).solve(); - }); + let mut problem = Problem::new(ode, bs3, controller); + black_box(problem.solve()); }); }); @@ -255,9 +242,8 @@ fn bench_tolerance_scaling(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, ()); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); } diff --git a/benches/orbit.rs b/benches/orbit.rs index de8cf26..8b4e321 100644 --- a/benches/orbit.rs +++ b/benches/orbit.rs @@ -29,9 +29,8 @@ fn bench_orbit(c: &mut Criterion) { c.bench_function("bench_orbit", |b| { b.iter(|| { - std::hint::black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + std::hint::black_box(problem.solve()); }); }); } diff --git a/benches/simple_1d.rs b/benches/simple_1d.rs index 3f4effc..f6a4e6f 100644 --- a/benches/simple_1d.rs +++ b/benches/simple_1d.rs @@ -20,9 +20,8 @@ fn bench_simple_1d(c: &mut Criterion) { c.bench_function("bench_simple_1d", |b| { b.iter(|| { - std::hint::black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + std::hint::black_box(problem.solve()); }); }); } @@ -44,11 +43,13 @@ fn bench_interpolation_1d(c: &mut Criterion) { c.bench_function("bench_interpolation_1d", |b| { b.iter(|| { - std::hint::black_box({ - let mut problem = Problem::new(ode, dp45, controller); - let solution = problem.solve(); - let _ = (0..100).map(|t| solution.interpolate(t as f64 * 0.1)[0]); - }); + let mut problem = Problem::new(ode, dp45, controller); + let solution = problem.solve(); + std::hint::black_box( + (0..100) + .map(|t| solution.interpolate(t as f64 * 0.1)[0]) + .sum::(), + ); }); }); } diff --git a/benches/vern7_comparison.rs b/benches/vern7_comparison.rs index 5e23cfe..7ae9749 100644 --- a/benches/vern7_comparison.rs +++ b/benches/vern7_comparison.rs @@ -28,9 +28,8 @@ fn bench_exponential_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 4.0, y0, ()); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, bs3, controller).solve(); - }); + let mut problem = Problem::new(ode, bs3, controller); + black_box(problem.solve()); }); }); @@ -38,9 +37,8 @@ fn bench_exponential_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 4.0, y0, ()); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -48,9 +46,8 @@ fn bench_exponential_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 4.0, y0, ()); let vern7 = Vern7::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, vern7, controller).solve(); - }); + let mut problem = Problem::new(ode, vern7, controller); + black_box(problem.solve()); }); }); @@ -77,9 +74,8 @@ fn bench_harmonic_oscillator_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 20.0, y0, ()); let bs3 = BS3::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, bs3, controller).solve(); - }); + let mut problem = Problem::new(ode, bs3, controller); + black_box(problem.solve()); }); }); @@ -87,9 +83,8 @@ fn bench_harmonic_oscillator_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 20.0, y0, ()); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -97,9 +92,8 @@ fn bench_harmonic_oscillator_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 20.0, y0, ()); let vern7 = Vern7::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, vern7, controller).solve(); - }); + let mut problem = Problem::new(ode, vern7, controller); + black_box(problem.solve()); }); }); @@ -138,9 +132,8 @@ fn bench_orbit_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -148,9 +141,8 @@ fn bench_orbit_tight_tol(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10000.0, y0, params); let vern7 = Vern7::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, vern7, controller).solve(); - }); + let mut problem = Problem::new(ode, vern7, controller); + black_box(problem.solve()); }); }); @@ -175,27 +167,23 @@ fn bench_vern7_interpolation(c: &mut Criterion) { // Vern7 with interpolation (should compute extra stages lazily) group.bench_function("vern7_with_interpolation", |b| { b.iter(|| { - black_box({ - let ode = ODE::new(&derivative, 0.0, 5.0, y0, ()); - let vern7 = Vern7::new().a_tol(tol).r_tol(tol); - let mut problem = Problem::new(ode, vern7, controller); - let solution = problem.solve(); - // Interpolate at 100 points - first one computes extra stages - let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect(); - }); + let ode = ODE::new(&derivative, 0.0, 5.0, y0, ()); + let vern7 = Vern7::new().a_tol(tol).r_tol(tol); + let mut problem = Problem::new(ode, vern7, controller); + let solution = problem.solve(); + // Interpolate at 100 points - first one computes extra stages + black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::>()); }); }); // DP5 with interpolation for comparison group.bench_function("dp5_with_interpolation", |b| { b.iter(|| { - black_box({ - let ode = ODE::new(&derivative, 0.0, 5.0, y0, ()); - let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); - let mut problem = Problem::new(ode, dp45, controller); - let solution = problem.solve(); - let _: Vec<_> = (0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect(); - }); + let ode = ODE::new(&derivative, 0.0, 5.0, y0, ()); + let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); + let mut problem = Problem::new(ode, dp45, controller); + let solution = problem.solve(); + black_box((0..100).map(|i| solution.interpolate(i as f64 * 0.05)).collect::>()); }); }); @@ -223,9 +211,8 @@ fn bench_tolerance_scaling_vern7(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, ()); let dp45 = DormandPrince45::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, dp45, controller).solve(); - }); + let mut problem = Problem::new(ode, dp45, controller); + black_box(problem.solve()); }); }); @@ -233,9 +220,8 @@ fn bench_tolerance_scaling_vern7(c: &mut Criterion) { let ode = ODE::new(&derivative, 0.0, 10.0, y0, ()); let vern7 = Vern7::new().a_tol(tol).r_tol(tol); b.iter(|| { - black_box({ - Problem::new(ode, vern7, controller).solve(); - }); + let mut problem = Problem::new(ode, vern7, controller); + black_box(problem.solve()); }); }); } diff --git a/readme.md b/readme.md index 8585ed6..76e6b58 100644 --- a/readme.md +++ b/readme.md @@ -15,8 +15,8 @@ and field line tracing: | **Vern7** (Verner) | 7(6) | 10+6 | 7th order | High accuracy (rtol ~ 1e-8 to 1e-12) | **Performance at 1e-10 tolerance:** -- Vern7: **2.7-8.8x faster** than DP5 -- Vern7: **50x+ faster** than BS3 +- Vern7: **1.1-4.4x faster** than DP5 +- Vern7: **53-83x faster** than BS3 See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analysis. @@ -27,6 +27,8 @@ See [benchmark report](VERN7_BENCHMARK_REPORT.md) for detailed performance analy - **Dense output interpolation** at any time point - **Parameters** in derivative and callback functions - **Lazy computation** of extra interpolation stages (Vern7) +- **PID step size control** as an alternative to the PI controller +- **Automatic Differentiation**, the same problem can be solved with floats or also differentiated with dual numbers ### Future Improvements @@ -68,15 +70,77 @@ let ode = ODE::new(&derivative, 0.0, 6.3, y0, params); let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6); let controller = PIController::default(); -let value_too_high = Callback { - event: &|t: f64, _y: Vector2, _p: &Params| { 5.0 - t }, +let value_too_high = FnCallback { + event: &|t: f64, _y: &Vector2, _p: &Params| { 5.0 - t }, effect: &stop, }; // Solve the problem -let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high); +let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high); let solution = problem.solve(); // Can interpolate solutions to whatever you want let _interpolated_answer = solution.interpolate(4.4); ``` + +## Automatic differentiation + +State, time and step size are generic over the `Scalar` trait. `f64` implements it by default, +and so do the forward-mode dual numbers of [`num-dual`](https://crates.io/crates/num-dual) (the +`autodiff` feature, on by default). Solving at a dual scalar propagates exact derivatives alongside +the state. + +The only change needed is to the signature of the right-hand side. Here is the pendulum from above, +before and after: + +```rust,ignore +fn derivative(_t: f64, y: Vector2, p: &Params) -> Vector2 { + let &(g, l) = p; + Vector2::new( y[1], -(g/l) * y[0].sin() ) +} + +fn derivative(_t: T, y: Vector2, p: &Params) -> Vector2 { + let &(g, l) = p; + Vector2::new( y[1], -y[0].sin() * (g/l) ) +} +``` + +`Params` stays `(f64, f64)` — parameters only need to become `T` themselves if you want to +differentiate *with respect to* them. Real constants mix into `T` arithmetic directly; `Scalar::lift` +is only needed to scale a vector by a real, or to build a `T` out of nothing. + +With that one change, the same function can be solved or differentiated: + +```rust +use nalgebra::{Matrix2, Vector2}; +use num_dual::DualSVec64; +use ordinary_diffeq::prelude::*; +use std::f64::consts::PI; + +type Params = (f64, f64); + +fn derivative(_t: T, y: Vector2, p: &Params) -> Vector2 { + let &(g, l) = p; + Vector2::new(y[1], -y[0].sin() * (g / l)) +} + +// One propagation in `DualSVec64<2>` yields the whole 2x2 Jacobian of the +// final state with respect to the initial one. +let (y_end, jacobian): (Vector2, Matrix2) = num_dual::jacobian( + |y0: Vector2>| { + type T = DualSVec64<2>; + let ode = ODE::new(&derivative, T::lift(0.0), T::lift(1.0), y0, (9.81, 1.0)); + let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12); + let mut problem = Problem::new(ode, vern7, PIController::default()); + let solution = problem.solve(); + *solution.states.last().unwrap() + }, + &Vector2::new(0.0, PI / 2.0), +); + +// The pendulum is Hamiltonian, so its flow is area-preserving +assert!((jacobian.determinant() - 1.0).abs() < 1e-9); +``` + +These are exact derivatives of the scheme as it ran, at machine precision, with no finite +differencing and no step size to tune. diff --git a/roadmap/README.md b/roadmap/README.md index 0e71e4d..34b7e8a 100644 --- a/roadmap/README.md +++ b/roadmap/README.md @@ -298,6 +298,12 @@ Each feature below links to a detailed implementation plan in the `features/` di - **Dependencies**: None - **Effort**: Small +- [ ] **[Reduce Per-Step Allocation](features/39-per-step-allocation.md)** + - Two heap allocations per accepted step (stage array, dense coefficients) + - Against ~150 ns of arithmetic; matters only for cheap right-hand sides + - **Dependencies**: None + - **Effort**: Medium + --- ## Implementation Notes diff --git a/roadmap/features/19-jacobian-computation.md b/roadmap/features/19-jacobian-computation.md index ff31370..a28eaea 100644 --- a/roadmap/features/19-jacobian-computation.md +++ b/roadmap/features/19-jacobian-computation.md @@ -4,6 +4,31 @@ Finite difference and auto-diff Jacobians +## Note - 2026-09-29 + +Forward-mode AD arrived by a different route than this entry assumes. On the +`feature/generic-scalar` branch, state, time and step size became generic over a +`Scalar` trait, implemented by `f64` and by `num-dual`'s dual numbers. That gives +differentiation *through the whole solve* - d y(t_f)/d y0, d/d t0, d/d t_end and +d/d params, so a full 6x6 state transition matrix comes out of one propagation. +This entry is framed around the narrower case, df/dy at a point for implicit +methods, which now falls out of the same machinery. + +What that leaves: + +- `rosenbrock::compute_jacobian` is still forward finite differences + (eps = sqrt(machine_eps) * max(|y[j]|, 1), D+1 evaluations, ~8 digits). + Pointing it at a dual evaluation is now a small change, not a large one. +- Sparse Jacobian support: untouched. +- The **Effort: Large** estimate below predates the genericisation and wants + revisiting. + +Deferred deliberately, as scope creep against the genericisation: how a Jacobian +should be obtained for whichever algorithms need one. Whether the integrator +requests it through a trait, whether AD or finite differences is a per-problem +choice, and how either interacts with the `JacobianUpdateStrategy` reuse policy +already in `rosenbrock.rs`. To be thought through properly on its own. + ## Why This Feature Matters (To be detailed) diff --git a/roadmap/features/39-per-step-allocation.md b/roadmap/features/39-per-step-allocation.md new file mode 100644 index 0000000..2ea2c0b --- /dev/null +++ b/roadmap/features/39-per-step-allocation.md @@ -0,0 +1,81 @@ +# Feature: Reduce Per-Step Allocation + +## Overview + +Remove the heap allocations every accepted step performs. + +## Why This Feature Matters + +Each call to `step()` allocates a `Vec` for the stage values +(`vec![SVector::zeros(); STAGES]`) and a second `Vec` for the dense output +coefficients. `Problem::solve` then wraps the latter in a `RefCell` and pushes +it onto a growing vector. That is two allocations per step against roughly +150 ns of arithmetic. + +Measured 2026-09-29 on the harmonic oscillator at 1e-10, where the right-hand +side is two flops so the timing is almost entirely integrator overhead: + +| Method | Steps | Total | Per step | +|--------|-------|-------|----------| +| BS3 | 13,927 | 1.58 ms | 113 ns | +| DP45 | 722 | 112 us | 155 ns | +| Vern7 | 124 | 26.6 us | 215 ns | + +`BS3::step` in isolation is 44.6 ns for four stages. + +The gain is largest exactly where it is least important: problems with cheap +right-hand sides. For the orbit propagation this crate is aimed at, a +spherical-harmonic gravity evaluation dominates and per-step overhead is +rounding error. Worth doing, not worth doing first. + +## Dependencies + +None. `STAGES` is already an associated const, so the stage array can be a +fixed-size array rather than a `Vec`. + +## Implementation Approach + +(To be detailed.) Sketch: stage storage becomes `[SVector; STAGES]` if +const generics allow it at the trait level, or a caller-supplied scratch buffer +threaded through `step()`. Dense coefficients want a pre-reserved buffer sized +from an estimate of the step count, or a flat `Vec>` indexed by +stride rather than a `Vec` of `Vec`s. + +Note the interaction with the `RefCell` in `Solution::dense`: it exists so +Vern7 can append its lazily computed extra stages. A flat buffer has to keep +that possible. + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Benchmarking + +## Testing Requirements + +Numerics must not change. The existing benches plus a step-count and +final-state comparison against the previous release are sufficient: results +should be bit-identical. + +## References + +1. Julia implementation: OrdinaryDiffEq.jl, which caches stage storage in the + integrator rather than allocating per step + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: Low - no change to the numerics, so the test suite is a strong guard + +## Success Criteria + +- [ ] No allocation in the steady-state step loop +- [ ] Bit-identical results against the previous release +- [ ] Measurable improvement on the cheap-RHS benchmarks + +## Future Enhancements + +Supersedes the "Optimize memory allocation in step method" bullet under Future +Enhancements in `01-bs3-method.md`. diff --git a/src/callback.rs b/src/callback.rs index c290f0a..5643513 100644 --- a/src/callback.rs +++ b/src/callback.rs @@ -1,21 +1,54 @@ use nalgebra::SVector; -use super::ode::ODE; +use crate::ode::ODE; +use crate::scalar::Scalar; -/// A function that takes in a time and a state and outputs a single float value +/// A callback function you can utilize to change your model mid-run without stopping. /// -/// The integration solver will check this function for zero crossings +/// The solver evaluates [`event`](Callback::event) at each end of every accepted +/// step and looks for a sign change. When it finds one it locates the crossing +/// and applies [`effect`](Callback::effect) to the ODE there. +/// ``` +/// # use ordinary_diffeq::prelude::*; +/// # use nalgebra::Vector2; +/// # type Params = (); +/// let too_high = FnCallback { +/// event: &|_t: f64, y: &Vector2, _p: &Params| 10.0 - y[0], +/// effect: &stop, +/// }; +/// ``` +pub trait Callback { + /// The function whose zero crossings mark the event. + fn event(&self, t: T, y: &SVector, params: &P) -> T; + + /// What to do to the ODE once the crossing has been located. + fn effect(&self, ode: &mut ODE<'_, T, D, P>); +} + +/// A [`Callback`] built from a pair of closures. #[derive(Clone, Copy)] -pub struct Callback<'a, const D: usize, P> { +pub struct FnCallback<'a, T: Scalar, const D: usize, P> { /// The function to check for zero crossings - pub event: &'a dyn Fn(f64, SVector, &P) -> f64, + pub event: &'a dyn Fn(T, &SVector, &P) -> T, /// The function to change the ODE - pub effect: &'a dyn Fn(&mut ODE), + pub effect: &'a dyn Fn(&mut ODE<'_, T, D, P>), +} + +impl Callback for FnCallback<'_, T, D, P> { + #[inline] + fn event(&self, t: T, y: &SVector, params: &P) -> T { + (self.event)(t, y, params) + } + + #[inline] + fn effect(&self, ode: &mut ODE<'_, T, D, P>) { + (self.effect)(ode) + } } /// A convenience function for stopping the integration -pub fn stop(ode: &mut ODE) { +pub fn stop(ode: &mut ODE<'_, T, D, P>) { ode.t_end = ode.t; } @@ -26,8 +59,8 @@ mod tests { #[test] fn test_basic_callbacks() { type Params = (); - let _value_too_high = Callback { - event: &|_: f64, y: SVector, _p: &Params| 10.0 - y[0], + let _value_too_high = FnCallback { + event: &|_: f64, y: &SVector, _p: &Params| 10.0 - y[0], effect: &stop, }; } diff --git a/src/controller.rs b/src/controller.rs index 6b65996..6ac9ecc 100644 --- a/src/controller.rs +++ b/src/controller.rs @@ -1,11 +1,13 @@ +use crate::scalar::Scalar; + #[derive(Debug, Clone, Copy, PartialEq)] -pub enum TryStep { - Accepted(f64, f64), - NotYetAccepted(f64), +pub enum TryStep { + Accepted(T, T), + NotYetAccepted(T), } -impl TryStep { - pub fn extract(&self) -> f64 { +impl TryStep { + pub fn extract(&self) -> T { match self { TryStep::Accepted(h, _) => *h, TryStep::NotYetAccepted(h) => *h, @@ -16,7 +18,7 @@ impl TryStep { matches!(self, TryStep::Accepted(_, _)) } - pub fn reset(&mut self) -> Result { + pub fn reset(&mut self) -> Result, &str> { match self { TryStep::Accepted(_, h) => Ok(TryStep::NotYetAccepted(*h)), TryStep::NotYetAccepted(_) => Err("Cannot reset a NotYetAccepted TryStep"), @@ -24,12 +26,18 @@ impl TryStep { } } -pub trait Controller { - fn determine_step(&mut self, h: f64, err: f64) -> TryStep; +/// Chooses the next step size from the last one and its normalized error. +pub trait Controller { + fn determine_step(&mut self, h: T, err: f64) -> TryStep; + + /// The step size to attempt next + fn next_step_guess(&self) -> TryStep; + + fn set_next_step_guess(&mut self, guess: TryStep); } #[derive(Debug, Clone, Copy)] -pub struct PIController { +pub struct PIController { pub alpha: f64, pub beta: f64, pub factor_c1: f64, @@ -37,13 +45,13 @@ pub struct PIController { pub factor_old: f64, pub h_max: f64, pub safety_factor: f64, - pub next_step_guess: TryStep, + pub next_step_guess: TryStep, } -impl Controller for PIController { +impl Controller for PIController { /// Determines if the previously run step size and error were valid or not. Either way, it also /// returns what the next step size should be - fn determine_step(&mut self, prev_step: f64, err: f64) -> TryStep { + fn determine_step(&mut self, prev_step: T, err: f64) -> TryStep { let factor_11 = err.powf(self.alpha); let factor = self.factor_c2.max( self.factor_c1 @@ -53,19 +61,30 @@ impl Controller for PIController { let mut h = prev_step / factor; // Accept the stepsize and provide what the next step size should be self.factor_old = err.max(1.0e-4); - if h.abs() > self.h_max { + if h.re().abs() > self.h_max { // If the step goes past the maximum allowed, though, we shrink it - h = self.h_max.copysign(h); + h = T::lift(self.h_max.copysign(h.re())); } TryStep::Accepted(prev_step, h) } else { // Reject the stepsize and propose a smaller one for the current step - TryStep::NotYetAccepted(prev_step / (self.factor_c1.min(factor_11 / self.safety_factor))) + TryStep::NotYetAccepted( + prev_step / self.factor_c1.min(factor_11 / self.safety_factor), + ) } } + + fn next_step_guess(&self) -> TryStep { + self.next_step_guess + } + + fn set_next_step_guess(&mut self, guess: TryStep) { + self.next_step_guess = guess; + } + } -impl PIController { +impl PIController { pub fn new( alpha: f64, beta: f64, @@ -83,17 +102,18 @@ impl PIController { factor_old: 1.0e-4, h_max: h_max.abs(), safety_factor, - next_step_guess: TryStep::NotYetAccepted(initial_h), + next_step_guess: TryStep::NotYetAccepted(T::lift(initial_h)), } } } -impl Default for PIController { +impl Default for PIController { fn default() -> Self { Self::new(0.17, 0.04, 10.0, 0.2, 100000.0, 0.9, 1e-4) } } + /// PID (Proportional-Integral-Derivative) step size controller /// /// The PID controller is an advanced adaptive time-stepping controller that provides @@ -111,9 +131,9 @@ impl Default for PIController { /// Where: /// - ε_n = normalized error estimate at current step /// - ε_{n-1} = normalized error estimate at previous step -/// - β₁ = proportional coefficient (controls reaction to current error) -/// - β₂ = integral coefficient (controls reaction to error history) -/// - β₃ = derivative coefficient (controls reaction to error rate of change) +/// - β₁ (`p_term`) = proportional coefficient (controls reaction to current error) +/// - β₂ (`i_term`) = integral coefficient (controls reaction to error history) +/// - β₃ (`d_term`) = derivative coefficient (controls reaction to error rate of change) /// /// # When to Use /// @@ -138,11 +158,11 @@ impl Default for PIController { /// let controller = PIDController::for_order(5); /// ``` #[derive(Debug, Clone, Copy)] -pub struct PIDController { +pub struct PIDController { // PID Coefficients - pub beta1: f64, // Proportional: reaction to current error - pub beta2: f64, // Integral: reaction to error history - pub beta3: f64, // Derivative: reaction to error rate of change + pub p_term: f64, // Proportional: reaction to current error + pub i_term: f64, // Integral: reaction to error history + pub d_term: f64, // Derivative: reaction to error rate of change // Constraints pub factor_c1: f64, // 1/min_factor (maximum step decrease) @@ -156,35 +176,35 @@ pub struct PIDController { pub h_old: f64, // h_{n-1}: previous step size // Next step guess - pub next_step_guess: TryStep, + pub next_step_guess: TryStep, } -impl Controller for PIDController { +impl Controller for PIDController { /// Determines if the previously run step was acceptable and computes the next step size /// using PID control theory - fn determine_step(&mut self, h: f64, err: f64) -> TryStep { + fn determine_step(&mut self, h: T, err: f64) -> TryStep { // Compute PID control factor // For first step or when history isn't available, fall back to simpler control let factor = if self.err_old <= 0.0 { // First step: use only proportional control - let factor_11 = err.powf(self.beta1); + let factor_11 = err.powf(self.p_term); self.factor_c2.max( self.factor_c1.min(factor_11 / self.safety_factor) ) } else if self.err_older <= 0.0 { // Second step: use PI control (proportional + integral) - let factor_11 = err.powf(self.beta1); - let factor_12 = self.err_old.powf(-self.beta2); + let factor_11 = err.powf(self.p_term); + let factor_12 = self.err_old.powf(-self.i_term); self.factor_c2.max( self.factor_c1.min(factor_11 * factor_12 / self.safety_factor) ) } else { // Full PID control (proportional + integral + derivative) - let factor_11 = err.powf(self.beta1); - let factor_12 = self.err_old.powf(-self.beta2); + let factor_11 = err.powf(self.p_term); + let factor_12 = self.err_old.powf(-self.i_term); // Derivative term uses ratio of consecutive step sizes let factor_13 = if self.h_old > 0.0 { - (h / self.h_old).powf(-self.beta3) + (h.re() / self.h_old).powf(-self.d_term) } else { 1.0 }; @@ -200,43 +220,53 @@ impl Controller for PIDController { // Update error history for next step self.err_older = self.err_old; self.err_old = err.max(1.0e-4); // Prevent very small values - self.h_old = h; + self.h_old = h.re(); // Apply maximum step size limit - if h_next.abs() > self.h_max { - h_next = self.h_max.copysign(h_next); + if h_next.re().abs() > self.h_max { + h_next = T::lift(self.h_max.copysign(h_next.re())); } TryStep::Accepted(h, h_next) } else { // Step rejected - propose smaller step // Use only proportional control for rejection (more aggressive) - let factor_11 = err.powf(self.beta1); + let factor_11 = err.powf(self.p_term); let h_next = h / (self.factor_c1.min(factor_11 / self.safety_factor)); // Note: Don't update history on rejection TryStep::NotYetAccepted(h_next) } } + + fn next_step_guess(&self) -> TryStep { + self.next_step_guess + } + + fn set_next_step_guess(&mut self, guess: TryStep) { + self.next_step_guess = guess; + } + } -impl PIDController { +impl PIDController { /// Create a new PID controller with custom parameters /// /// # Arguments /// - /// * `beta1` - Proportional coefficient (typically 0.3-0.5) - /// * `beta2` - Integral coefficient (typically 0.04-0.1) - /// * `beta3` - Derivative coefficient (typically 0.01-0.05) + /// * `p_term` - Proportional coefficient (typically 0.3-0.5) + /// * `i_term` - Integral coefficient (typically 0.04-0.1) + /// * `d_term` - Derivative coefficient (typically 0.01-0.05) /// * `max_factor` - Maximum step size increase factor (typically 10.0) /// * `min_factor` - Maximum step size decrease factor (typically 0.2) /// * `h_max` - Maximum allowed step size /// * `safety_factor` - Safety factor (typically 0.9) /// * `initial_h` - Initial step size guess + #[allow(clippy::too_many_arguments)] pub fn new( - beta1: f64, - beta2: f64, - beta3: f64, + p_term: f64, + i_term: f64, + d_term: f64, max_factor: f64, min_factor: f64, h_max: f64, @@ -244,9 +274,9 @@ impl PIDController { initial_h: f64, ) -> Self { Self { - beta1, - beta2, - beta3, + p_term, + i_term, + d_term, factor_c1: 1.0 / min_factor, factor_c2: 1.0 / max_factor, h_max: h_max.abs(), @@ -254,16 +284,16 @@ impl PIDController { err_old: 0.0, // No history initially err_older: 0.0, // No history initially h_old: 0.0, // No history initially - next_step_guess: TryStep::NotYetAccepted(initial_h), + next_step_guess: TryStep::NotYetAccepted(T::lift(initial_h)), } } /// Create a PID controller with coefficients scaled for a specific method order /// /// Uses the Gustafsson coefficients scaled by order: - /// - β₁ = 0.49 / (order + 1) - /// - β₂ = 0.34 / (order + 1) - /// - β₃ = 0.10 / (order + 1) + /// - `p_term` = 0.49 / (order + 1) + /// - `i_term` = 0.34 / (order + 1) + /// - `d_term` = 0.10 / (order + 1) /// /// # Arguments /// @@ -271,9 +301,9 @@ impl PIDController { pub fn for_order(order: usize) -> Self { let k_plus_1 = (order + 1) as f64; Self::new( - 0.49 / k_plus_1, // beta1: proportional - 0.34 / k_plus_1, // beta2: integral - 0.10 / k_plus_1, // beta3: derivative + 0.49 / k_plus_1, // p_term: proportional + 0.34 / k_plus_1, // i_term: integral + 0.10 / k_plus_1, // d_term: derivative 10.0, // max_factor 0.2, // min_factor 100000.0, // h_max @@ -292,22 +322,22 @@ impl PIDController { } } -impl Default for PIDController { +impl Default for PIDController { /// Default PID controller using conservative coefficients /// /// Uses conservative PID coefficients that provide stable performance /// across a wide range of problems: - /// - β₁ = 0.07 (proportional) - /// - β₂ = 0.04 (integral) - /// - β₃ = 0.01 (derivative) + /// - `p_term` = 0.07 (proportional) + /// - `i_term` = 0.04 (integral) + /// - `d_term` = 0.01 (derivative) /// /// These values are appropriate for typical ODE methods of order 5-7. /// For method-specific tuning, use `PIDController::for_order(order)` instead. fn default() -> Self { Self::new( - 0.07, // beta1 (proportional) - 0.04, // beta2 (integral) - 0.01, // beta3 (derivative) + 0.07, // p_term (proportional) + 0.04, // i_term (integral) + 0.01, // d_term (derivative) 10.0, // max_factor 0.2, // min_factor 100000.0, // h_max @@ -317,13 +347,14 @@ impl Default for PIDController { } } + #[cfg(test)] mod tests { use super::*; #[test] fn test_pi_controller_creation() { - let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 10.0, 0.9, 1e-4); + let controller = PIController::::new(0.17, 0.04, 10.0, 0.2, 10.0, 0.9, 1e-4); assert!(controller.alpha == 0.17); assert!(controller.beta == 0.04); @@ -337,11 +368,11 @@ mod tests { #[test] fn test_pid_controller_creation() { - let controller = PIDController::new(0.49, 0.34, 0.10, 10.0, 0.2, 10.0, 0.9, 1e-4); + let controller = PIDController::::new(0.49, 0.34, 0.10, 10.0, 0.2, 10.0, 0.9, 1e-4); - assert_eq!(controller.beta1, 0.49); - assert_eq!(controller.beta2, 0.34); - assert_eq!(controller.beta3, 0.10); + assert_eq!(controller.p_term, 0.49); + assert_eq!(controller.i_term, 0.34); + assert_eq!(controller.d_term, 0.10); assert_eq!(controller.factor_c1, 1.0 / 0.2); assert_eq!(controller.factor_c2, 1.0 / 10.0); assert_eq!(controller.h_max, 10.0); @@ -353,32 +384,33 @@ mod tests { #[test] fn test_pid_for_order() { - let controller = PIDController::for_order(5); + let controller = PIDController::::for_order(5); // For order 5, k+1 = 6 - assert!((controller.beta1 - 0.49 / 6.0).abs() < 1e-10); - assert!((controller.beta2 - 0.34 / 6.0).abs() < 1e-10); - assert!((controller.beta3 - 0.10 / 6.0).abs() < 1e-10); + assert!((controller.p_term - 0.49 / 6.0).abs() < 1e-10); + assert!((controller.i_term - 0.34 / 6.0).abs() < 1e-10); + assert!((controller.d_term - 0.10 / 6.0).abs() < 1e-10); } #[test] fn test_pid_default() { - let controller = PIDController::default(); + let controller = PIDController::::default(); // Default uses conservative PID coefficients - assert_eq!(controller.beta1, 0.07); - assert_eq!(controller.beta2, 0.04); - assert_eq!(controller.beta3, 0.01); // True PID with derivative term + assert_eq!(controller.p_term, 0.07); + assert_eq!(controller.i_term, 0.04); + assert_eq!(controller.d_term, 0.01); // True PID with derivative term } #[test] fn test_pid_reset() { - let mut controller = PIDController::default(); - // Simulate some history - controller.err_old = 0.5; - controller.err_older = 0.3; - controller.h_old = 0.01; + let mut controller = PIDController:: { + err_old: 0.5, + err_older: 0.3, + h_old: 0.01, + ..Default::default() + }; controller.reset(); @@ -392,8 +424,8 @@ mod tests { // For smooth problems, PI and PID should perform similarly // Test with exponential decay: y' = -y - let mut pi = PIController::default(); - let mut pid = PIDController::default(); + let mut pi = PIController::::default(); + let mut pid = PIDController::::default(); // Simulate a sequence of small errors (smooth problem) let errors = vec![0.8, 0.6, 0.5, 0.45, 0.4, 0.35, 0.3]; @@ -403,8 +435,8 @@ mod tests { let mut pid_steps = Vec::new(); for &err in &errors { - let mut pi_result = >::determine_step(&mut pi, h, err); - let mut pid_result = >::determine_step(&mut pid, h, err); + let mut pi_result = >::determine_step(&mut pi, h, err); + let mut pid_result = >::determine_step(&mut pid, h, err); if pi_result.is_accepted() { pi_steps.push(pi_result.extract()); @@ -437,25 +469,25 @@ mod tests { #[test] fn test_pid_bootstrap() { // Test that PID progressively uses P → PI → PID as history builds - let mut pid = PIDController::new(0.49, 0.34, 0.10, 10.0, 0.2, 100.0, 0.9, 0.01); + let mut pid = PIDController::::new(0.49, 0.34, 0.10, 10.0, 0.2, 100.0, 0.9, 0.01); let h = 0.01; let err1 = 0.5; let err2 = 0.4; let err3 = 0.3; - // First step: should use only proportional (beta1) + // First step: should use only proportional (p_term) assert_eq!(pid.err_old, 0.0); assert_eq!(pid.err_older, 0.0); - let step1 = >::determine_step(&mut pid, h, err1); + let step1 = >::determine_step(&mut pid, h, err1); assert!(step1.is_accepted()); // After first step, err_old is updated but err_older is still 0 assert!(pid.err_old > 0.0); assert_eq!(pid.err_older, 0.0); - // Second step: should use PI (beta1 and beta2) - let step2 = >::determine_step(&mut pid, h, err2); + // Second step: should use PI (p_term and i_term) + let step2 = >::determine_step(&mut pid, h, err2); assert!(step2.is_accepted()); // After second step, both err_old and err_older should be set @@ -463,21 +495,21 @@ mod tests { assert!(pid.err_older > 0.0); assert!(pid.h_old > 0.0); - // Third step: should use full PID (beta1, beta2, and beta3) - let step3 = >::determine_step(&mut pid, h, err3); + // Third step: should use full PID (p_term, i_term, and d_term) + let step3 = >::determine_step(&mut pid, h, err3); assert!(step3.is_accepted()); } #[test] fn test_pid_step_rejection() { // Test that error history is NOT updated after rejection - let mut pid = PIDController::default(); + let mut pid = PIDController::::default(); let h = 0.01; // First, accept a step to build history let err_good = 0.5; - let step1 = >::determine_step(&mut pid, h, err_good); + let step1 = >::determine_step(&mut pid, h, err_good); assert!(step1.is_accepted()); let err_old_before = pid.err_old; @@ -486,7 +518,7 @@ mod tests { // Now reject a step with large error let err_bad = 2.0; - let step2 = >::determine_step(&mut pid, h, err_bad); + let step2 = >::determine_step(&mut pid, h, err_bad); assert!(!step2.is_accepted()); // History should NOT have changed after rejection @@ -500,9 +532,9 @@ mod tests { // Test on oscillatory error pattern (simulating difficult problem) // True PID (with derivative term) should provide smoother step size evolution - let mut pi = PIController::default(); - // Use actual PID with non-zero beta3 (Gustafsson coefficients for order 5) - let mut pid = PIDController::for_order(5); + let mut pi = PIController::::default(); + // Use actual PID with non-zero d_term (Gustafsson coefficients for order 5) + let mut pid = PIDController::::for_order(5); // Simulate oscillatory error pattern let errors = vec![0.8, 0.3, 0.9, 0.2, 0.85, 0.25, 0.8, 0.3]; @@ -515,8 +547,8 @@ mod tests { let mut pid_h_prev = h; for &err in &errors { - let mut pi_result = >::determine_step(&mut pi, h, err); - let mut pid_result = >::determine_step(&mut pid, h, err); + let mut pi_result = >::determine_step(&mut pi, h, err); + let mut pid_result = >::determine_step(&mut pid, h, err); if pi_result.is_accepted() { let pi_h_next = pi_result.reset().unwrap().extract(); @@ -550,7 +582,7 @@ mod tests { / pid_ratios.len() as f64; let pid_std = pid_variance.sqrt(); - // With true PID (non-zero beta3), we expect similar or better stability + // With true PID (non-zero d_term), we expect similar or better stability // Allow some tolerance since this is a simple synthetic test assert!( pid_std <= pi_std * 1.1, diff --git a/src/integrator/bs3.rs b/src/integrator/bs3.rs index ff0aab5..cecd0d1 100644 --- a/src/integrator/bs3.rs +++ b/src/integrator/bs3.rs @@ -1,7 +1,8 @@ use nalgebra::SVector; use super::super::ode::ODE; -use super::Integrator; +use super::{Integrator, StepResult}; +use crate::scalar::Scalar; /// Bogacki-Shampine 3/2 integrator trait for tableau coefficients pub trait BS3Integrator<'a> { @@ -57,6 +58,15 @@ pub struct BS3 { r_tol: f64, } +impl Default for BS3 +where + BS3: Integrator, + { + fn default() -> Self { + Self::new() + } +} + impl BS3 where BS3: Integrator, @@ -140,38 +150,52 @@ where const ADAPTIVE: bool = true; const DENSE: bool = true; - fn step

( - &self, - ode: &ODE, - h: f64, - ) -> (SVector, Option, Option>>) { + fn step(&self, ode: &ODE, h: T) -> StepResult { // Allocate storage for the 4 stages - let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; + let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; // Stage 1: k1 = f(t, y) k[0] = (ode.f)(ode.t, ode.y, &ode.params); // Stage 2: k2 = f(t + c[1]*h, y + h*a[1,0]*k1) - let y2 = ode.y + h * Self::A[0] * k[0]; - k[1] = (ode.f)(ode.t + Self::C[1] * h, y2, &ode.params); + let y2 = ode.y + k[0] * (h * Self::A[0]); + k[1] = (ode.f)(ode.t + h * Self::C[1], y2, &ode.params); // Stage 3: k3 = f(t + c[2]*h, y + h*(a[2,0]*k1 + a[2,1]*k2)) - let y3 = ode.y + h * (Self::A[1] * k[0] + Self::A[2] * k[1]); - k[2] = (ode.f)(ode.t + Self::C[2] * h, y3, &ode.params); + let y3 = ode.y + (k[0] * T::lift(Self::A[1]) + k[1] * T::lift(Self::A[2])) * h; + k[2] = (ode.f)(ode.t + h * Self::C[2], y3, &ode.params); // Stage 4: k4 = f(t + c[3]*h, y + h*(a[3,0]*k1 + a[3,1]*k2 + a[3,2]*k3)) - let y4 = ode.y + h * (Self::A[3] * k[0] + Self::A[4] * k[1] + Self::A[5] * k[2]); - k[3] = (ode.f)(ode.t + Self::C[3] * h, y4, &ode.params); + let y4 = ode.y + + (k[0] * T::lift(Self::A[3]) + k[1] * T::lift(Self::A[4]) + k[2] * T::lift(Self::A[5])) + * h; + k[3] = (ode.f)(ode.t + h * Self::C[3], y4, &ode.params); // Compute 3rd order solution - let next_y = ode.y + h * (Self::B[0] * k[0] + Self::B[1] * k[1] + Self::B[2] * k[2] + Self::B[3] * k[3]); + let next_y = ode.y + + (k[0] * T::lift(Self::B[0]) + + k[1] * T::lift(Self::B[1]) + + k[2] * T::lift(Self::B[2]) + + k[3] * T::lift(Self::B[3])) + * h; // Compute error estimate (difference between 3rd and 2nd order solutions) - let err = h * (Self::B_ERROR[0] * k[0] + Self::B_ERROR[1] * k[1] + Self::B_ERROR[2] * k[2] + Self::B_ERROR[3] * k[3]); + let err = (k[0] * T::lift(Self::B_ERROR[0]) + + k[1] * T::lift(Self::B_ERROR[1]) + + k[2] * T::lift(Self::B_ERROR[2]) + + k[3] * T::lift(Self::B_ERROR[3])) + * h; - // Compute error norm scaled by tolerance - let tol = self.a_tol + ode.y.abs() * self.r_tol; - let error_norm = (err.component_div(&tol)).norm(); + // Compute error norm scaled by tolerance. Step control is a real-valued + // decision, so each component drops to the reals as it is accumulated — + // one pass, rather than building real copies of `err` and `y` first. + let mut error_norm_sq = 0.0; + for i in 0..D { + let tol = self.a_tol[i] + ode.y[i].re().abs() * self.r_tol; + let scaled = err[i].re() / tol; + error_norm_sq += scaled * scaled; + } + let error_norm = error_norm_sq.sqrt(); // Store coefficients for dense output (cubic Hermite interpolation) // BS3 uses standard cubic Hermite interpolation with derivatives at endpoints @@ -183,16 +207,20 @@ where k[3], // f(t1, y1) - derivative at end (FSAL) ]; - (next_y, Some(error_norm), Some(dense_coeffs)) + StepResult { + y: next_y, + err: Some(error_norm), + dense: Some(dense_coeffs), + } } - fn interpolate( + fn interpolate( &self, - t_start: f64, - t_end: f64, - dense: &[SVector], - t: f64, - ) -> SVector { + t_start: T, + t_end: T, + dense: &[SVector], + t: T, + ) -> SVector { // Compute interpolation parameter θ ∈ [0, 1] let theta = (t - t_start) / (t_end - t_start); let h = t_end - t_start; @@ -215,14 +243,14 @@ where let f1 = &dense[3]; let theta2 = theta * theta; - let one_minus_theta = 1.0 - theta; + let one_minus_theta = T::lift(1.0) - theta; let one_minus_theta2 = one_minus_theta * one_minus_theta; // Apply cubic Hermite interpolation formula - (1.0 + 2.0 * theta) * one_minus_theta2 * y0 - + theta2 * (3.0 - 2.0 * theta) * y1 - + theta * one_minus_theta2 * h * f0 - + theta2 * (theta - 1.0) * h * f1 + y0 * ((theta * 2.0 + 1.0) * one_minus_theta2) + + y1 * (theta2 * (-(theta * 2.0 - 3.0))) + + f0 * (theta * one_minus_theta2 * h) + + f1 * (theta2 * (theta - 1.0) * h) } } @@ -233,6 +261,7 @@ mod tests { use nalgebra::Vector1; #[test] + #[allow(clippy::assertions_on_constants)] fn test_bs3_creation() { let _bs3: BS3<1> = BS3::new(); assert_eq!(BS3::<1>::ORDER, 3); @@ -254,7 +283,8 @@ mod tests { let bs3 = BS3::new(); let h = 0.001; // Smaller step size for tighter tolerances - let (y_next, err, dense) = bs3.step(&ode, h); + let __step = bs3.step(&ode, h); + let (y_next, err, dense) = (__step.y, __step.err, __step.dense); // At t=0.001, exact solution is e^0.001 ≈ 1.0010005001667084 let exact = (0.001_f64).exp(); @@ -283,7 +313,8 @@ mod tests { let bs3 = BS3::new(); let h = 0.001; // Smaller step size - let (_y_next, _err, dense) = bs3.step(&ode, h); + let __step = bs3.step(&ode, h); + let (_y_next, _err, dense) = (__step.y, __step.err, __step.dense); let dense = dense.unwrap(); // Interpolate at midpoint @@ -312,7 +343,8 @@ mod tests { // Take 100 steps to reach t = 1.0 let mut ode = ODE::new(&derivative, 0.0, 1.0, y0, ()); for _ in 0..100 { - let (y_new, _, _) = bs3.step(&ode, h); + let __step = bs3.step(&ode, h); + let (y_new, _, _) = (__step.y, __step.err, __step.dense); ode.y = y_new; ode.t += h; } @@ -345,7 +377,8 @@ mod tests { // Take steps until we reach t_end while ode.t < t_end - 1e-10 { - let (y_new, _, _) = bs3.step(&ode, h); + let __step = bs3.step(&ode, h); + let (y_new, _, _) = (__step.y, __step.err, __step.dense); ode.y = y_new; ode.t += h; } @@ -389,7 +422,8 @@ mod tests { // First step let ode1 = ODE::new(&derivative, 0.0, 1.0, y0, ()); - let (y_new1, _, dense1) = bs3.step(&ode1, h); + let __step = bs3.step(&ode1, h); + let (y_new1, _, dense1) = (__step.y, __step.err, __step.dense); let dense1 = dense1.unwrap(); // Extract f1 from first step (derivative at end of step) @@ -397,7 +431,8 @@ mod tests { // Second step starts where first ended let ode2 = ODE::new(&derivative, h, 1.0, y_new1, ()); - let (_, _, dense2) = bs3.step(&ode2, h); + let __step = bs3.step(&ode2, h); + let (_, _, dense2) = (__step.y, __step.err, __step.dense); let dense2 = dense2.unwrap(); // Extract f0 from second step (derivative at start of step) diff --git a/src/integrator/dormand_prince.rs b/src/integrator/dormand_prince.rs index 7e3004a..5c39b4b 100644 --- a/src/integrator/dormand_prince.rs +++ b/src/integrator/dormand_prince.rs @@ -1,7 +1,8 @@ use nalgebra::SVector; -use super::super::ode::ODE; -use super::Integrator; +use super::{Integrator, StepResult}; +use crate::ode::ODE; +use crate::scalar::Scalar; /// Integrator Trait pub trait DormandPrinceIntegrator<'a> { @@ -11,9 +12,54 @@ pub trait DormandPrinceIntegrator<'a> { const D: &'a [f64]; } +/// Dormand-Prince 5(4) method +/// +/// A 5th order explicit Runge-Kutta method with an embedded 4th order method for +/// error estimation. This is the general-purpose workhorse of the crate, and the +/// same method as MATLAB's `ode45` and `DifferentialEquations.jl`'s `DP5`. +/// +/// # Characteristics +/// - Order: 5(4) - 5th order solution with 4th order error estimate +/// - Stages: 7 +/// - FSAL: Yes (First Same As Last - the 7th stage is the next step's first, +/// though this implementation does not yet reuse it) +/// - Adaptive: Yes +/// - Dense output: 4th order interpolation +/// +/// # When to use DormandPrince45 +/// - The default choice when nothing about the problem suggests otherwise +/// - Problems requiring general accuracy (rtol ~ 1e-6 to 1e-8) +/// - Non-stiff problems +/// +/// At tighter tolerances than about 1e-8, prefer [`Vern7`](super::vern7::Vern7): +/// a higher order method reaches the same accuracy in far fewer steps. +/// +/// # Example +/// ```rust +/// use ordinary_diffeq::prelude::*; +/// use nalgebra::Vector1; +/// +/// let params = (); +/// fn derivative(_t: f64, y: Vector1, _p: &()) -> Vector1 { +/// Vector1::new(-y[0]) +/// } +/// +/// let y0 = Vector1::new(1.0); +/// let ode = ODE::new(&derivative, 0.0, 5.0, y0, ()); +/// let dp45 = DormandPrince45::new().a_tol(1e-8).r_tol(1e-8); +/// let controller = PIController::default(); +/// +/// let mut problem = Problem::new(ode, dp45, controller); +/// let solution = problem.solve(); +/// ``` +/// +/// # References +/// - Dormand, J.R. and Prince, P.J. (1980), "A family of embedded Runge-Kutta +/// formulae", Journal of Computational and Applied Mathematics, Vol. 6, No. 1, +/// pp. 19-26 #[derive(Debug, Clone, Copy)] pub struct DormandPrince45 { - a_tol: SVector, + a_tol: SVector, r_tol: f64, } @@ -22,16 +68,43 @@ where DormandPrince45: Integrator, { pub fn new() -> Self { - Self { a_tol: SVector::::from_element(1e-8), r_tol: 1e-8 } + Self { + a_tol: SVector::::from_element(1e-8), + r_tol: 1e-8, + } } - pub fn a_tol(&mut self, a_tol: f64) -> Self { - Self { a_tol: SVector::::from_element(a_tol), r_tol: self.r_tol } + + /// The same absolute tolerance on every component. + pub fn a_tol(self, a_tol: f64) -> Self { + Self { + a_tol: SVector::::from_element(a_tol), + r_tol: self.r_tol, + } } - pub fn a_tol_full(&mut self, a_tol: SVector::) -> Self { - Self { a_tol, r_tol: self.r_tol } + + /// A per-component absolute tolerance. + pub fn a_tol_full(self, a_tol: SVector) -> Self { + Self { + a_tol, + r_tol: self.r_tol, + } } - pub fn r_tol(&mut self, r_tol: f64) -> Self { - Self { a_tol: self.a_tol, r_tol } + + /// The same relative tolerance on every component + pub fn r_tol(self, r_tol: f64) -> Self { + Self { + a_tol: self.a_tol, + r_tol, + } + } +} + +impl Default for DormandPrince45 +where + DormandPrince45: Integrator, +{ + fn default() -> Self { + Self::new() } } @@ -96,51 +169,58 @@ where const ADAPTIVE: bool = true; const DENSE: bool = true; - fn step

( - &self, - ode: &ODE, - h: f64, - ) -> (SVector, Option, Option>>) { - let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; + fn step(&self, ode: &ODE, h: T) -> StepResult { + let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; let mut next_y = ode.y; - let mut err = SVector::::zeros(); - let mut rcont5 = SVector::::zeros(); + let mut err = SVector::::zeros(); + let mut rcont5 = SVector::::zeros(); // Do the first of the summations k[0] = (ode.f)(ode.t, ode.y, &ode.params); - next_y += k[0] * Self::B[0] * h; - err += k[0] * (Self::B[0] - Self::B[Self::STAGES]) * h; + next_y += k[0] * T::lift(Self::B[0]) * h; + err += k[0] * T::lift(Self::B[0] - Self::B[Self::STAGES]) * h; let rcont1 = ode.y; - rcont5 += k[0] * h * Self::D[0]; + rcont5 += k[0] * h * T::lift(Self::D[0]); // Then the rest for i in 1..Self::STAGES { // Compute the ks - let mut y_term = SVector::::zeros(); + let mut y_term = SVector::::zeros(); for (j, item) in k.iter().enumerate().take(i) { - y_term += item * Self::A[(i * (i - 1)) / 2 + j]; + y_term += item * T::lift(Self::A[(i * (i - 1)) / 2 + j]); } - k[i] = (ode.f)(ode.t + Self::C[i] * h, ode.y + y_term * h, &ode.params); + k[i] = (ode.f)(ode.t + T::lift(Self::C[i]) * h, ode.y + y_term * h, &ode.params); // Use that and bis to calculate the y and error terms - next_y += k[i] * h * Self::B[i]; - err += k[i] * (Self::B[i] - Self::B[i + Self::STAGES]) * h; - rcont5 += k[i] * h * Self::D[i]; + next_y += k[i] * h * T::lift(Self::B[i]); + err += k[i] * T::lift(Self::B[i] - Self::B[i + Self::STAGES]) * h; + rcont5 += k[i] * h * T::lift(Self::D[i]); } let rcont2 = next_y - ode.y; - let rcont3 = h * k[0] - rcont2; + let rcont3 = k[0] * h - rcont2; let rcont4 = rcont2 - k[Self::STAGES - 1] * h - rcont3; - let tol = self.a_tol + ode.y * self.r_tol; let rcont = vec![rcont1, rcont2, rcont3, rcont4, rcont5]; - (next_y, Some((err.component_div(&tol)).norm()), Some(rcont)) + // Just take the real component of the error + let mut error_norm_sq = 0.0; + for i in 0..D { + let tol = self.a_tol[i] + ode.y[i].re() * self.r_tol; + let scaled = err[i].re() / tol; + error_norm_sq += scaled * scaled; + } + StepResult { + y: next_y, + err: Some(error_norm_sq.sqrt()), + dense: Some(rcont), + } } - fn interpolate( + + fn interpolate( &self, - t_start: f64, - t_end: f64, - dense: &[SVector], - t: f64, - ) -> SVector { + t_start: T, + t_end: T, + dense: &[SVector], + t: T, + ) -> SVector { let s = (t - t_start) / (t_end - t_start); - let s1 = 1.0 - s; + let s1 = T::one() - s; dense[0] + (dense[1] + (dense[2] + (dense[3] + dense[4] * s1) * s) * s1) * s } } diff --git a/src/integrator/mod.rs b/src/integrator/mod.rs index d49c070..b5796c9 100644 --- a/src/integrator/mod.rs +++ b/src/integrator/mod.rs @@ -1,12 +1,34 @@ use nalgebra::SVector; -use super::ode::ODE; +use crate::ode::ODE; +use crate::scalar::Scalar; pub mod bs3; pub mod dormand_prince; pub mod rosenbrock; pub mod vern7; +/// Result of a `Integrator::step` +#[derive(Debug, Clone, PartialEq)] +pub struct StepResult { + /// The state at the end of the step. + pub y: SVector, + + /// The local error, normalized by the integrator's tolerances, or `None` + /// from a fixed-step method. A value at or below `1.0` means the step is + /// acceptable. + pub err: Option, + + /// Coefficients for interpolating inside the step, to be handed back to + /// [`Integrator::interpolate`], or `None` from a method without dense + /// output. + /// + /// For a method with lazy dense output (`EXTRA_STAGES > 0`) the layout is + /// `[y_start, y_end, k_1, ..., k_MAIN_STAGES]`, and the extra stages are + /// appended by [`Integrator::compute_extra_stages`] on first use. + pub dense: Option>>, +} + /// Integrator Trait pub trait Integrator { const ORDER: usize; @@ -23,21 +45,18 @@ pub trait Integrator { /// Total stages when full dense output is computed const TOTAL_DENSE_STAGES: usize = Self::MAIN_STAGES + Self::EXTRA_STAGES; - /// Returns a new y value, then possibly an error value, and possibly a dense output - /// coefficient set - fn step

( - &self, - ode: &ODE, - h: f64, - ) -> (SVector, Option, Option>>); + /// Advances `ode` by `h` from its current state, without mutating it. + fn step(&self, ode: &ODE, h: T) -> StepResult; - fn interpolate( + /// Evaluates the solution at `t` inside a step that ran from `t_start` to + /// `t_end`, using the `dense` coefficients that step produced. + fn interpolate( &self, - t_start: f64, - t_end: f64, - dense: &[SVector], - t: f64, - ) -> SVector; + t_start: T, + t_end: T, + dense: &[SVector], + t: T, + ) -> SVector; /// Compute extra stages for full-order dense output (lazy computation). /// @@ -56,14 +75,14 @@ pub trait Integrator { /// # Returns /// /// Vector of extra k-stages (empty for most integrators) - fn compute_extra_stages

( + fn compute_extra_stages( &self, - _ode: &ODE, - _t_start: f64, - _y_start: SVector, - _h: f64, - _main_stages: &[SVector], - ) -> Vec> { + _ode: &ODE, + _t_start: T, + _y_start: SVector, + _h: T, + _main_stages: &[SVector], + ) -> Vec> { // Default implementation: no extra stages needed Vec::new() } @@ -79,7 +98,7 @@ mod tests { #[test] fn test_dopri() { type Params = (); - fn derivative(_t: f64, y: Vector3, _p: &Params) -> Vector3 { + fn derivative(_t: T, y: Vector3, _p: &Params) -> Vector3 { y } @@ -92,11 +111,11 @@ mod tests { // and also that the error seems reasonable let step = 0.001; while ode.t < ode.t_end { - let (new_y, err, _) = dp45.step(&ode, step); - ode.y = new_y; + let result = dp45.step(&ode, step); + ode.y = result.y; ode.t += step; assert_relative_eq!(ode.y[0], ode.t.exp(), max_relative = 0.01); - assert!(err.unwrap() < 1.0); + assert!(result.err.unwrap() < 1.0); } } } diff --git a/src/integrator/rosenbrock.rs b/src/integrator/rosenbrock.rs index ddfc95d..a620a56 100644 --- a/src/integrator/rosenbrock.rs +++ b/src/integrator/rosenbrock.rs @@ -1,7 +1,8 @@ use nalgebra::{SMatrix, SVector}; -use super::super::ode::ODE; -use super::Integrator; +use super::super::ode::{ProblemFunction, ODE}; +use super::{Integrator, StepResult}; +use crate::scalar::Scalar; /// Strategy for when to update the Jacobian matrix #[derive(Debug, Clone, Copy)] @@ -31,24 +32,26 @@ impl Default for JacobianUpdateStrategy { /// /// Uses forward differences: J[i,j] ≈ (f_i(y + ε*e_j) - f_i(y)) / ε /// where ε = √(machine_epsilon) * max(|y[j]|, 1.0) -pub fn compute_jacobian( - f: &dyn Fn(f64, SVector, &P) -> SVector, - t: f64, - y: SVector, +pub fn compute_jacobian( + f: ProblemFunction<'_, T, D, P>, + t: T, + y: SVector, params: &P, -) -> SMatrix { +) -> SMatrix { let sqrt_eps = f64::EPSILON.sqrt(); let f_y = f(t, y, params); - let mut jacobian = SMatrix::::zeros(); + let mut jacobian = SMatrix::::zeros(); // Compute each column of the Jacobian by perturbing y[j] for j in 0..D { - // Choose epsilon based on the magnitude of y[j] - let epsilon = sqrt_eps * y[j].abs().max(1.0); + // Choose epsilon based on the magnitude of y[j]. The perturbation is a + // numerical-method parameter, not part of the solution, so it is chosen + // in the reals. + let epsilon = sqrt_eps * y[j].re().abs().max(1.0); // Perturb y in the j-th direction let mut y_perturbed = y; - y_perturbed[j] += epsilon; + y_perturbed[j] += T::lift(epsilon); // Evaluate f at perturbed point let f_perturbed = f(t, y_perturbed, params); @@ -131,8 +134,8 @@ pub fn compute_jacobian( /// let rosenbrock = Rosenbrock23::new(); /// /// // Take a single step -/// let (y_next, error, _dense) = rosenbrock.step(&ode, 0.1); -/// assert!((y_next[0] - 0.905).abs() < 0.01); +/// let step = rosenbrock.step(&ode, 0.1); +/// assert!((step.y[0] - 0.905).abs() < 0.01); /// ``` #[derive(Debug, Clone, Copy)] pub struct Rosenbrock23 { @@ -218,81 +221,61 @@ impl Integrator for Rosenbrock23 { const ADAPTIVE: bool = true; const DENSE: bool = true; - fn step

( - &self, - ode: &ODE, - h: f64, - ) -> (SVector, Option, Option>>) { + fn step(&self, ode: &ODE, h: T) -> StepResult { let t = ode.t; let uprev = ode.y; // Compute Jacobian let j = compute_jacobian(&ode.f, t, uprev, &ode.params); - // Julia: dtγ = dt * d let dtgamma = h * self.d; // Form W = I - dtγ*J - let w = SMatrix::::identity() - dtgamma * j; + let w = SMatrix::::identity() - j * dtgamma; let w_inv = w.try_inverse().expect("W matrix is singular"); // Evaluate fsalfirst = f(uprev) let fsalfirst = (ode.f)(t, uprev, &ode.params); // Stage 1: Solve W * k₁ = f(y) where k₁ is a derivative estimate - // Julia stores derivatives in k, not displacements let k1 = w_inv * fsalfirst; // Stage 2: u = uprev + dt/2 * k₁ - // Julia line 69 let dto2 = h / 2.0; - let u = uprev + dto2 * k1; + let u = uprev + k1 * dto2; // Evaluate f₁ = f(u, t + dt/2) - // Julia line 71 let f1 = (ode.f)(t + dto2, u, &ode.params); // Stage 2: W * k₂ = f₁ - k₁ + J*k₁ - // Julia line 80: linsolve_tmp = f₁ - tmp (where tmp = k₁) - // This is equivalent to: W * k₂ = f₁ - k₁ - // => (I - dtγ*J) * k₂ = f₁ - k₁ - // => k₂ = (I - dtγ*J)^{-1} * (f₁ - k₁) - // But actually, maybe the RHS should be scaled differently. Let me try: W * k₂ = f₁ + J*k₁ - // Since W = I - dtγ*J, we have W*k₂ - I*k₂ = -dtγ*J*k₂, so if RHS = f₁ + J*k₁... - // Actually, let's just implement exactly what Julia does: let rhs2 = f1 - k1; let k2_temp = w_inv * rhs2; - // Julia then does: k₂ = k₂_temp * neginvdtγ + k₁ - // But neginvdtγ = -1/(dtγ), which would give huge values. - // Let me try without that scaling: let k2 = k2_temp + k1; // Solution: u = uprev + dt * k₂ - // Julia line 89 - let u_final = uprev + h * k2; + let u_final = uprev + k2 * h; // Error estimation // Evaluate fsallast = f(u_final, t + dt) - // Julia line 94 let fsallast = (ode.f)(t + h, u_final, &ode.params); - // Julia line 98-99: linsolve_tmp = fsallast - c₃₂*(k₂ - f₁) - 2*(k₁ - fsalfirst) + dt*dT // Ignoring dT (time derivative) for autonomous systems - let linsolve_tmp3 = fsallast - self.c32 * (k2 - f1) - 2.0 * (k1 - fsalfirst); + let linsolve_tmp3 = + fsallast - (k2 - f1) * T::lift(self.c32) - (k1 - fsalfirst) * T::lift(2.0); // Stage 3 for error estimation: W * k₃ = linsolve_tmp3 let k3 = w_inv * linsolve_tmp3; // Error: dt/6 * (k₁ - 2*k₂ + k₃) - // Julia line 115 let dto6 = h / 6.0; - let error_vec = dto6 * (k1 - 2.0 * k2 + k3); + let error_vec = (k1 - k2 * T::lift(2.0) + k3) * dto6; // Compute scalar error estimate using weighted norm let mut error_sum = 0.0; for i in 0..D { - let scale = self.a_tol + self.r_tol * uprev[i].abs().max(u_final[i].abs()); - let weighted_error = error_vec[i] / scale; + let scale = + self.a_tol + self.r_tol * uprev[i].re().abs().max(u_final[i].re().abs()); + let weighted_error = error_vec[i].re() / scale; error_sum += weighted_error * weighted_error; } let error = (error_sum / D as f64).sqrt(); @@ -300,16 +283,20 @@ impl Integrator for Rosenbrock23 { // Dense output: store k₁ and k₂ let dense = Some(vec![k1, k2]); - (u_final, Some(error), dense) + StepResult { + y: u_final, + err: Some(error), + dense, + } } - fn interpolate( + fn interpolate( &self, - t_start: f64, - t_end: f64, - dense: &[SVector], - t: f64, - ) -> SVector { + t_start: T, + t_end: T, + dense: &[SVector], + t: T, + ) -> SVector { // Second order interpolation using k₁ and k₂ // For Rosenbrock methods, we use a simple Hermite interpolation let k1 = dense[0]; @@ -323,7 +310,7 @@ impl Integrator for Rosenbrock23 { // y(t) ≈ y_n + θ*h*((1-θ)*k₁ + θ*k₂) // But we don't have y_n stored, so we just return the stage interpolation // This is a simple linear interpolation of the derivative - (1.0 - theta) * k1 + theta * k2 + k1 * (T::lift(1.0) - theta) + k2 * theta } } @@ -383,7 +370,8 @@ mod tests { let ode = ODE::new(&derivative, 0.0, 0.1, y0, ()); let rb23 = Rosenbrock23::new(); - let (y_next, err, _) = rb23.step(&ode, 0.1); + let __step = rb23.step(&ode, 0.1); + let (y_next, err, _) = (__step.y, __step.err, __step.dense); let analytical = (-0.1_f64).exp(); println!("Computed: {}, Analytical: {}", y_next[0], analytical); @@ -416,7 +404,8 @@ mod tests { let rb23 = Rosenbrock23::new(); while ode.t < t_end - 1e-10 { - let (y_next, _, _) = rb23.step(&ode, h); + let __step = rb23.step(&ode, h); + let (y_next, _, _) = (__step.y, __step.err, __step.dense); ode.y = y_next; ode.t += h; } @@ -465,7 +454,8 @@ mod tests { let rb23 = Rosenbrock23::new(); for _ in 0..n_steps { - let (y_next, _, _) = rb23.step(&ode, actual_h); + let __step = rb23.step(&ode, actual_h); + let (y_next, _, _) = (__step.y, __step.err, __step.dense); ode.y = y_next; ode.t += actual_h; } @@ -497,7 +487,6 @@ mod tests { #[test] fn test_rosenbrock23_adaptive_solve() { - // Julia test: sol = solve(prob, Rosenbrock23()); @test length(sol) < 20 // This tests that the adaptive solver can efficiently solve prob_ode_linear use crate::controller::PIController; use crate::problem::Problem; @@ -518,7 +507,6 @@ mod tests { println!("Adaptive solve completed in {} steps", solution.states.len()); - // Julia test: @test length(sol) < 20 assert!(solution.states.len() < 20, "Adaptive solve should complete in less than 20 steps, got {}", solution.states.len()); diff --git a/src/integrator/vern7.rs b/src/integrator/vern7.rs index 45f9568..ac5ec91 100644 --- a/src/integrator/vern7.rs +++ b/src/integrator/vern7.rs @@ -1,7 +1,8 @@ use nalgebra::SVector; use super::super::ode::ODE; -use super::Integrator; +use super::{Integrator, StepResult}; +use crate::scalar::Scalar; /// Verner 7 integrator trait for tableau coefficients pub trait Vern7Integrator<'a> { @@ -69,6 +70,15 @@ pub struct Vern7 { r_tol: f64, } +impl Default for Vern7 +where + Vern7: Integrator, + { + fn default() -> Self { + Self::new() + } +} + impl Vern7 where Vern7: Integrator, @@ -278,13 +288,9 @@ where const MAIN_STAGES: usize = 10; const EXTRA_STAGES: usize = 6; - fn step

( - &self, - ode: &ODE, - h: f64, - ) -> (SVector, Option, Option>>) { + fn step(&self, ode: &ODE, h: T) -> StepResult { // Allocate storage for the 10 stages - let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; + let mut k: Vec> = vec![SVector::::zeros(); Self::STAGES]; // Stage 1: k[0] = f(t, y) k[0] = (ode.f)(ode.t, ode.y, &ode.params); @@ -295,43 +301,54 @@ where // A matrix is stored in lower triangular form, row by row // Row i has i elements (0-indexed), starting at position i*(i-1)/2 let row_start = (i * (i - 1)) / 2; - for j in 0..i { - y_temp += k[j] * Self::A[row_start + j] * h; + for (j, kj) in k.iter().enumerate().take(i) { + y_temp += kj * T::lift(Self::A[row_start + j]) * h; } - k[i] = (ode.f)(ode.t + Self::C[i] * h, y_temp, &ode.params); + k[i] = (ode.f)(ode.t + h * Self::C[i], y_temp, &ode.params); } // Compute 7th order solution using B weights let mut next_y = ode.y; - for i in 0..Self::STAGES { - next_y += k[i] * Self::B[i] * h; + for (i, ki) in k.iter().enumerate() { + next_y += ki * T::lift(Self::B[i]) * h; } // Compute error estimate using B_ERROR weights - let mut err = SVector::::zeros(); - for i in 0..Self::STAGES { - err += k[i] * Self::B_ERROR[i] * h; + let mut err = SVector::::zeros(); + for (i, ki) in k.iter().enumerate() { + err += ki * T::lift(Self::B_ERROR[i]) * h; } - // Compute error norm scaled by tolerance - let tol = self.a_tol + ode.y.abs() * self.r_tol; - let error_norm = (err.component_div(&tol)).norm(); + // Compute error norm scaled by tolerance. Step control is a + // real-valued decision, so each component drops to the reals as it is + // accumulated rather than building real copies of `err` and `y` first. + let mut error_norm_sq = 0.0; + for i in 0..D { + let tol = self.a_tol[i] + ode.y[i].re().abs() * self.r_tol; + let scaled = err[i].re() / tol; + error_norm_sq += scaled * scaled; + } + let error_norm = error_norm_sq.sqrt(); // Store dense output coefficients // For now, store all k values for interpolation let mut dense_coeffs = vec![ode.y, next_y]; dense_coeffs.extend_from_slice(&k); - (next_y, Some(error_norm), Some(dense_coeffs)) + StepResult { + y: next_y, + err: Some(error_norm), + dense: Some(dense_coeffs), + } } - fn interpolate( + fn interpolate( &self, - t_start: f64, - t_end: f64, - dense: &[SVector], - t: f64, - ) -> SVector { + t_start: T, + t_end: T, + dense: &[SVector], + t: T, + ) -> SVector { // Vern7 uses 7th order polynomial interpolation // Check if extra stages (k11-k16) are available // Dense array format: [y0, y1, k1, k2, ..., k10, k11, ..., k16] @@ -358,8 +375,8 @@ where // Helper to evaluate polynomial using Horner's method #[inline] - fn evalpoly(x: f64, coeffs: &[f64]) -> f64 { - let mut result = 0.0; + fn evalpoly(x: T, coeffs: &[f64]) -> T { + let mut result = T::lift(0.0); for &c in coeffs.iter().rev() { result = result * x + c; } @@ -433,13 +450,13 @@ where ]); // Compute base interpolation with main stages - let mut result = y0 + h * (k1 * b1_theta + + let mut result = y0 + (k1 * b1_theta + k4 * b4_theta + k5 * b5_theta + k6 * b6_theta + k7 * b7_theta + k8 * b8_theta + - k9 * b9_theta); + k9 * b9_theta) * h; // If extra stages are available, add their contribution for full 7th order accuracy if dense.len() >= 2 + Self::TOTAL_DENSE_STAGES { @@ -507,25 +524,25 @@ where ]); // Add contribution from extra stages - result += h * (k11 * b11_theta + + result += (k11 * b11_theta + k12 * b12_theta + k13 * b13_theta + k14 * b14_theta + k15 * b15_theta + - k16 * b16_theta); + k16 * b16_theta) * h; } result } - fn compute_extra_stages

( + fn compute_extra_stages( &self, - ode: &ODE, - t_start: f64, - y_start: SVector, - h: f64, - main_stages: &[SVector], - ) -> Vec> { + ode: &ODE, + t_start: T, + y_start: SVector, + h: T, + main_stages: &[SVector], + ) -> Vec> { // Extract main stages that are used in extra stage computation // From Julia: extra stages use k1, k4, k5, k6, k7, k8, k9 let k1 = &main_stages[0]; @@ -540,86 +557,86 @@ where // Stage 11: uses k1, k4-k9 (7 coefficients) let mut y11 = y_start; - y11 += k1 * Self::A_EXTRA[0] * h; - y11 += k4 * Self::A_EXTRA[1] * h; - y11 += k5 * Self::A_EXTRA[2] * h; - y11 += k6 * Self::A_EXTRA[3] * h; - y11 += k7 * Self::A_EXTRA[4] * h; - y11 += k8 * Self::A_EXTRA[5] * h; - y11 += k9 * Self::A_EXTRA[6] * h; - let k11 = (ode.f)(t_start + Self::C_EXTRA[0] * h, y11, &ode.params); + y11 += k1 * T::lift(Self::A_EXTRA[0]) * h; + y11 += k4 * T::lift(Self::A_EXTRA[1]) * h; + y11 += k5 * T::lift(Self::A_EXTRA[2]) * h; + y11 += k6 * T::lift(Self::A_EXTRA[3]) * h; + y11 += k7 * T::lift(Self::A_EXTRA[4]) * h; + y11 += k8 * T::lift(Self::A_EXTRA[5]) * h; + y11 += k9 * T::lift(Self::A_EXTRA[6]) * h; + let k11 = (ode.f)(t_start + h * Self::C_EXTRA[0], y11, &ode.params); extra_stages.push(k11); // Stage 12: uses k1, k4-k9, k11 (8 coefficients) let mut y12 = y_start; - y12 += k1 * Self::A_EXTRA[7] * h; - y12 += k4 * Self::A_EXTRA[8] * h; - y12 += k5 * Self::A_EXTRA[9] * h; - y12 += k6 * Self::A_EXTRA[10] * h; - y12 += k7 * Self::A_EXTRA[11] * h; - y12 += k8 * Self::A_EXTRA[12] * h; - y12 += k9 * Self::A_EXTRA[13] * h; - y12 += &extra_stages[0] * Self::A_EXTRA[14] * h; // k11 - let k12 = (ode.f)(t_start + Self::C_EXTRA[1] * h, y12, &ode.params); + y12 += k1 * T::lift(Self::A_EXTRA[7]) * h; + y12 += k4 * T::lift(Self::A_EXTRA[8]) * h; + y12 += k5 * T::lift(Self::A_EXTRA[9]) * h; + y12 += k6 * T::lift(Self::A_EXTRA[10]) * h; + y12 += k7 * T::lift(Self::A_EXTRA[11]) * h; + y12 += k8 * T::lift(Self::A_EXTRA[12]) * h; + y12 += k9 * T::lift(Self::A_EXTRA[13]) * h; + y12 += extra_stages[0] * T::lift(Self::A_EXTRA[14]) * h; // k11 + let k12 = (ode.f)(t_start + h * Self::C_EXTRA[1], y12, &ode.params); extra_stages.push(k12); // Stage 13: uses k1, k4-k9, k11, k12 (9 coefficients) let mut y13 = y_start; - y13 += k1 * Self::A_EXTRA[15] * h; - y13 += k4 * Self::A_EXTRA[16] * h; - y13 += k5 * Self::A_EXTRA[17] * h; - y13 += k6 * Self::A_EXTRA[18] * h; - y13 += k7 * Self::A_EXTRA[19] * h; - y13 += k8 * Self::A_EXTRA[20] * h; - y13 += k9 * Self::A_EXTRA[21] * h; - y13 += &extra_stages[0] * Self::A_EXTRA[22] * h; // k11 - y13 += &extra_stages[1] * Self::A_EXTRA[23] * h; // k12 - let k13 = (ode.f)(t_start + Self::C_EXTRA[2] * h, y13, &ode.params); + y13 += k1 * T::lift(Self::A_EXTRA[15]) * h; + y13 += k4 * T::lift(Self::A_EXTRA[16]) * h; + y13 += k5 * T::lift(Self::A_EXTRA[17]) * h; + y13 += k6 * T::lift(Self::A_EXTRA[18]) * h; + y13 += k7 * T::lift(Self::A_EXTRA[19]) * h; + y13 += k8 * T::lift(Self::A_EXTRA[20]) * h; + y13 += k9 * T::lift(Self::A_EXTRA[21]) * h; + y13 += extra_stages[0] * T::lift(Self::A_EXTRA[22]) * h; // k11 + y13 += extra_stages[1] * T::lift(Self::A_EXTRA[23]) * h; // k12 + let k13 = (ode.f)(t_start + h * Self::C_EXTRA[2], y13, &ode.params); extra_stages.push(k13); // Stage 14: uses k1, k4-k9, k11, k12, k13 (10 coefficients) let mut y14 = y_start; - y14 += k1 * Self::A_EXTRA[24] * h; - y14 += k4 * Self::A_EXTRA[25] * h; - y14 += k5 * Self::A_EXTRA[26] * h; - y14 += k6 * Self::A_EXTRA[27] * h; - y14 += k7 * Self::A_EXTRA[28] * h; - y14 += k8 * Self::A_EXTRA[29] * h; - y14 += k9 * Self::A_EXTRA[30] * h; - y14 += &extra_stages[0] * Self::A_EXTRA[31] * h; // k11 - y14 += &extra_stages[1] * Self::A_EXTRA[32] * h; // k12 - y14 += &extra_stages[2] * Self::A_EXTRA[33] * h; // k13 - let k14 = (ode.f)(t_start + Self::C_EXTRA[3] * h, y14, &ode.params); + y14 += k1 * T::lift(Self::A_EXTRA[24]) * h; + y14 += k4 * T::lift(Self::A_EXTRA[25]) * h; + y14 += k5 * T::lift(Self::A_EXTRA[26]) * h; + y14 += k6 * T::lift(Self::A_EXTRA[27]) * h; + y14 += k7 * T::lift(Self::A_EXTRA[28]) * h; + y14 += k8 * T::lift(Self::A_EXTRA[29]) * h; + y14 += k9 * T::lift(Self::A_EXTRA[30]) * h; + y14 += extra_stages[0] * T::lift(Self::A_EXTRA[31]) * h; // k11 + y14 += extra_stages[1] * T::lift(Self::A_EXTRA[32]) * h; // k12 + y14 += extra_stages[2] * T::lift(Self::A_EXTRA[33]) * h; // k13 + let k14 = (ode.f)(t_start + h * Self::C_EXTRA[3], y14, &ode.params); extra_stages.push(k14); // Stage 15: uses k1, k4-k9, k11, k12, k13 (10 coefficients, reuses k13 not k14) let mut y15 = y_start; - y15 += k1 * Self::A_EXTRA[34] * h; - y15 += k4 * Self::A_EXTRA[35] * h; - y15 += k5 * Self::A_EXTRA[36] * h; - y15 += k6 * Self::A_EXTRA[37] * h; - y15 += k7 * Self::A_EXTRA[38] * h; - y15 += k8 * Self::A_EXTRA[39] * h; - y15 += k9 * Self::A_EXTRA[40] * h; - y15 += &extra_stages[0] * Self::A_EXTRA[41] * h; // k11 - y15 += &extra_stages[1] * Self::A_EXTRA[42] * h; // k12 - y15 += &extra_stages[2] * Self::A_EXTRA[43] * h; // k13 - let k15 = (ode.f)(t_start + Self::C_EXTRA[4] * h, y15, &ode.params); + y15 += k1 * T::lift(Self::A_EXTRA[34]) * h; + y15 += k4 * T::lift(Self::A_EXTRA[35]) * h; + y15 += k5 * T::lift(Self::A_EXTRA[36]) * h; + y15 += k6 * T::lift(Self::A_EXTRA[37]) * h; + y15 += k7 * T::lift(Self::A_EXTRA[38]) * h; + y15 += k8 * T::lift(Self::A_EXTRA[39]) * h; + y15 += k9 * T::lift(Self::A_EXTRA[40]) * h; + y15 += extra_stages[0] * T::lift(Self::A_EXTRA[41]) * h; // k11 + y15 += extra_stages[1] * T::lift(Self::A_EXTRA[42]) * h; // k12 + y15 += extra_stages[2] * T::lift(Self::A_EXTRA[43]) * h; // k13 + let k15 = (ode.f)(t_start + h * Self::C_EXTRA[4], y15, &ode.params); extra_stages.push(k15); // Stage 16: uses k1, k4-k9, k11, k12, k13 (10 coefficients, reuses k13 not k14 or k15) let mut y16 = y_start; - y16 += k1 * Self::A_EXTRA[44] * h; - y16 += k4 * Self::A_EXTRA[45] * h; - y16 += k5 * Self::A_EXTRA[46] * h; - y16 += k6 * Self::A_EXTRA[47] * h; - y16 += k7 * Self::A_EXTRA[48] * h; - y16 += k8 * Self::A_EXTRA[49] * h; - y16 += k9 * Self::A_EXTRA[50] * h; - y16 += &extra_stages[0] * Self::A_EXTRA[51] * h; // k11 - y16 += &extra_stages[1] * Self::A_EXTRA[52] * h; // k12 - y16 += &extra_stages[2] * Self::A_EXTRA[53] * h; // k13 - let k16 = (ode.f)(t_start + Self::C_EXTRA[5] * h, y16, &ode.params); + y16 += k1 * T::lift(Self::A_EXTRA[44]) * h; + y16 += k4 * T::lift(Self::A_EXTRA[45]) * h; + y16 += k5 * T::lift(Self::A_EXTRA[46]) * h; + y16 += k6 * T::lift(Self::A_EXTRA[47]) * h; + y16 += k7 * T::lift(Self::A_EXTRA[48]) * h; + y16 += k8 * T::lift(Self::A_EXTRA[49]) * h; + y16 += k9 * T::lift(Self::A_EXTRA[50]) * h; + y16 += extra_stages[0] * T::lift(Self::A_EXTRA[51]) * h; // k11 + y16 += extra_stages[1] * T::lift(Self::A_EXTRA[52]) * h; // k12 + y16 += extra_stages[2] * T::lift(Self::A_EXTRA[53]) * h; // k13 + let k16 = (ode.f)(t_start + h * Self::C_EXTRA[5], y16, &ode.params); extra_stages.push(k16); extra_stages @@ -706,7 +723,8 @@ mod tests { while ode.t < t_end { let h_step = h.min(t_end - ode.t); - let (next_y, _, _) = vern7.step(&ode, h_step); + let __step = vern7.step(&ode, h_step); + let (next_y, _, _) = (__step.y, __step.err, __step.dense); ode.y = next_y; ode.t += h_step; } diff --git a/src/lib.rs b/src/lib.rs index 628e72a..6e2e39e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,20 +1,27 @@ #![allow(dead_code)] +// The readme is the crate's front page, and including it here means every +// example in it is compiled and run by `cargo test`. Requires `autodiff` +// because one of those examples differentiates. +#![cfg_attr(feature = "autodiff", doc = include_str!("../readme.md"))] pub mod callback; pub mod controller; pub mod integrator; pub mod ode; pub mod problem; +pub mod scalar; pub mod prelude { - pub use super::callback::{stop, Callback}; + pub use super::callback::{stop, Callback, FnCallback}; pub use super::controller::{PIController, PIDController}; pub use super::integrator::bs3::BS3; pub use super::integrator::dormand_prince::DormandPrince45; pub use super::integrator::rosenbrock::Rosenbrock23; pub use super::integrator::vern7::Vern7; pub use super::ode::ODE; + pub use super::integrator::{Integrator, StepResult}; pub use super::problem::{Problem, Solution}; + pub use super::scalar::Scalar; } #[cfg(test)] @@ -44,13 +51,13 @@ mod tests { let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-6); let controller = PIController::default(); - let value_too_high = Callback { - event: &|t: f64, _y: Vector2, _p: &Params| 5.0 - t, + let value_too_high = FnCallback { + event: &|t: f64, _y: &Vector2, _p: &Params| 5.0 - t, effect: &stop, }; // Solve the problem - let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high); + let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high); let solution = problem.solve(); // Can interpolate solutions to whatever you want @@ -150,3 +157,193 @@ mod tests { ); } } + +/// Tests that the derivatives a dual instantiation produces are the right ones. +#[cfg(all(test, feature = "autodiff"))] +mod autodiff_tests { + use crate::ode::ProblemFunction; + use crate::prelude::*; + use approx::assert_relative_eq; + use nalgebra::{Matrix6, SVector, Vector1, Vector6}; + use num_dual::{Dual64, DualSVec64}; + use std::f64::consts::PI; + + const MU: f64 = 3.98600441500000e14; + const SMA: f64 = 6.7781363e6; + + /// The two-body right-hand side, shared by the orbit tests below. + fn two_body(_t: T, state: Vector6, p: &(T,)) -> Vector6 { + let r = state.fixed_rows::<3>(0); + let acc = -(r * p.0) / r.norm().powi(3); + Vector6::new(state[3], state[4], state[5], acc[0], acc[1], acc[2]) + } + + fn orbit_y0() -> Vector6 { + Vector6::new( + T::lift(4.263868426884883e6), + T::lift(5.146189057155391e6), + T::lift(1.1310208421331816e6), + T::lift(-5923.454461876975), + T::lift(4496.802639690076), + T::lift(1870.3893008991558), + ) + } + + fn solve_with + Copy>( + f: ProblemFunction<'_, T, D, P>, + integrator: S, + t0: T, + t_end: T, + y0: SVector, + params: P, + ) -> SVector { + let ode = ODE::new(f, t0, t_end, y0, params); + let controller = PIController::new(0.37, 0.04, 10.0, 0.2, 1000.0, 0.9, 0.01); + let mut problem = Problem::new(ode, integrator, controller); + let solution = problem.solve(); + *solution.states.last().unwrap() + } + + /// The machine-precision case. + /// + /// A fifth-order Runge-Kutta method integrates `y' = p t^4` exactly, so both + /// the state and its derivative with respect to `p` are exact up to + /// rounding — there is no truncation error to hide behind. The analytic + /// answers are `y(t) = p t^5 / 5` and `dy/dp = t^5 / 5`. + #[test] + fn test_ad_exact_quadrature() { + fn derivative(t: T, _y: Vector1, p: &(T,)) -> Vector1 { + Vector1::new(p.0 * t.powi(4)) + } + + let t_f = 3.0; + // Seed the derivative on the parameter p. + let p = Dual64::from_re(2.0).derivative(); + let y = solve_with( + &derivative, + DormandPrince45::new().a_tol(1e-12).r_tol(1e-12), + Dual64::lift(0.0), + Dual64::lift(t_f), + Vector1::new(Dual64::lift(0.0)), + (p,), + ); + + assert_relative_eq!(y[0].re, 2.0 * t_f.powi(5) / 5.0, max_relative = 1e-14); + assert_relative_eq!(y[0].eps, t_f.powi(5) / 5.0, max_relative = 1e-14); + } + + /// `y' = 5y - 3` has the closed-form solution + /// `y(t) = (y0 - 0.6) e^{5(t - t0)} + 0.6`, so `dy(t_f)/dy0 = e^{5(t_f - t0)}` + /// exactly. Here the comparison is limited by the solver's own truncation + /// error rather than by rounding, which is the point: AD returns the exact + /// derivative of the numerical solution, and the numerical solution is + /// accurate to tolerance. + #[test] + fn test_ad_matches_analytic_sensitivity() { + fn derivative(_t: T, y: Vector1, _p: &()) -> Vector1 { + Vector1::new(y[0] * 5.0 - 3.0) + } + + let y0 = Dual64::from_re(1.0).derivative(); + let y = solve_with( + &derivative, + DormandPrince45::new().a_tol(1e-14).r_tol(1e-14), + Dual64::lift(2.0), + Dual64::lift(3.0), + Vector1::new(y0), + (), + ); + + assert_relative_eq!(y[0].re, 0.4 * 5.0f64.exp() + 0.6, max_relative = 1e-11); + assert_relative_eq!(y[0].eps, 5.0f64.exp(), max_relative = 1e-11); + } + + /// A gradient on the orbit problem, checked against a property of the + /// two-body flow that is known analytically. + /// + /// `DualSVec64<6>` carries all six partials at once, so one propagation + /// yields the full 6x6 state transition matrix over a period. Two-body + /// motion is Hamiltonian, so by Liouville's theorem its flow is + /// volume-preserving and `det(Phi) = 1` exactly — for any interval, for any + /// mu. A finite-difference Jacobian satisfies this only to the accuracy of + /// the difference formula; an exact one satisfies it to the accuracy of the + /// propagation. + #[test] + fn test_ad_orbit_state_transition_matrix() { + let period = 2.0 * PI * (SMA.powi(3) / MU).sqrt(); + let y0: Vector6 = orbit_y0(); + + // One propagation in `DualSVec64<6>` yields the whole state transition matrix. + let (y_end, stm): (Vector6, Matrix6) = num_dual::jacobian( + |y: Vector6>| { + type T = DualSVec64<6>; + solve_with( + &two_body, + DormandPrince45::new().a_tol(1e-12).r_tol(1e-12), + T::lift(0.0), + T::lift(period), + y, + (T::lift(MU),), + ) + }, + &y0, + ); + + // One period returns the state to where it started. + for i in 0..6 { + assert_relative_eq!(y_end[i], y0[i], max_relative = 1e-10); + } + + assert_relative_eq!(stm.determinant(), 1.0, max_relative = 1e-9); + + // And cross-check the whole matrix against central differences, which + // agree only to about the square root of machine epsilon. + let propagate = |y_start: Vector6| { + solve_with( + &two_body, + DormandPrince45::new().a_tol(1e-12).r_tol(1e-12), + 0.0, + period, + y_start, + (MU,), + ) + }; + for j in 0..6 { + let h = 1e-5 * y0[j].abs(); + let mut plus = y0; + let mut minus = y0; + plus[j] += h; + minus[j] -= h; + let fd = (propagate(plus) - propagate(minus)) / (2.0 * h); + for i in 0..6 { + assert_relative_eq!(stm[(i, j)], fd[i], max_relative = 1e-5, epsilon = 1e-7); + } + } + } + + /// The same, through Vern7 rather than Dormand-Prince, which exercises a + /// path Dormand-Prince does not have: Vern7 keeps six of its sixteen stages + /// back and computes them only when something interpolates inside a step. + /// That lazy computation calls the right-hand side after the solve has + /// finished, so it has to work in dual arithmetic too. + #[test] + fn test_ad_vern7_lazy_dense_output() { + fn derivative(_t: T, y: Vector1, _p: &()) -> Vector1 { + -y + } + + let y0 = Dual64::from_re(1.0).derivative(); + let ode = ODE::new(&derivative, Dual64::lift(0.0), Dual64::lift(2.0), Vector1::new(y0), ()); + let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12); + let mut problem = Problem::new(ode, vern7, PIController::default()); + let solution = problem.solve(); + + // Interpolating inside a step is what triggers the extra stages. + let t = Dual64::lift(1.3); + let y = solution.interpolate(t)[0]; + + // y(t) = y0 e^-t, so dy/dy0 = e^-t. + assert_relative_eq!(y.re, (-1.3f64).exp(), max_relative = 1e-10); + assert_relative_eq!(y.eps, (-1.3f64).exp(), max_relative = 1e-10); + } +} diff --git a/src/ode.rs b/src/ode.rs index 6ac6061..0fad207 100644 --- a/src/ode.rs +++ b/src/ode.rs @@ -1,27 +1,37 @@ use nalgebra::SVector; -type ProblemFunction<'a, const D: usize, P> = &'a dyn Fn(f64, SVector, &P) -> SVector; +use crate::scalar::Scalar; -/// The basic ODE object that will be passed around. The type (T) and the size (D) will be +/// The right-hand side of the system: `dy/dt = f(t, y, params)`. +/// +/// Write it generically over the scalar — `fn f(t: T, y: SVector, p: &P)` — +/// and the same function can be instantiated at `f64` to solve the problem and +/// at a dual number to differentiate it. +/// +/// Or just write it for f64 if you want to keep things simple +pub type ProblemFunction<'a, T, const D: usize, P> = + &'a dyn Fn(T, SVector, &P) -> SVector; + +/// The basic ODE object that will be passed around. The scalar (T) and the size (D) will be /// determined upon creation of the object #[derive(Clone, Copy)] -pub struct ODE<'a, const D: usize, P> { - pub f: ProblemFunction<'a, D, P>, - pub y: SVector, - pub t: f64, +pub struct ODE<'a, T: Scalar, const D: usize, P> { + pub f: ProblemFunction<'a, T, D, P>, + pub y: SVector, + pub t: T, pub params: P, - pub t0: f64, - pub t_end: f64, - pub h: f64, + pub t0: T, + pub t_end: T, + pub h: T, pub finished: bool, } -impl<'a, const D: usize, P> ODE<'a, D, P> { +impl<'a, T: Scalar, const D: usize, P> ODE<'a, T, D, P> { pub fn new( - f: ProblemFunction<'a, D, P>, - t0: f64, - t_end: f64, - y0: SVector, + f: ProblemFunction<'a, T, D, P>, + t0: T, + t_end: T, + y0: SVector, params: P, ) -> Self { Self { @@ -31,7 +41,7 @@ impl<'a, const D: usize, P> ODE<'a, D, P> { params, t0, t_end, - h: 0.001, + h: T::lift(0.001), finished: false, } } @@ -45,7 +55,7 @@ mod tests { #[test] fn test_ode_creation() { type Params = (); - fn derivative(_t: f64, y: Vector3, _p: &Params) -> Vector3 { + fn derivative(_t: T, y: Vector3, _p: &Params) -> Vector3 { -y } @@ -64,7 +74,7 @@ mod tests { type Params = (f64, bool); let params = (34.0, true); - fn derivative(t: f64, y: Vector3, p: &Params) -> Vector3 { + fn derivative(t: T, y: Vector3, p: &Params) -> Vector3 { if p.1 { -y } else { diff --git a/src/problem.rs b/src/problem.rs index eccedcf..885a65d 100644 --- a/src/problem.rs +++ b/src/problem.rs @@ -1,28 +1,33 @@ +//! Assembling an ODE, an integrator and a controller + use nalgebra::SVector; use roots::{find_root_brent, SimpleConvergency}; use std::cell::RefCell; use super::callback::Callback; -use super::controller::{Controller, PIController, TryStep}; +use super::controller::{Controller, TryStep}; use super::integrator::Integrator; use super::ode::ODE; +use super::scalar::Scalar; #[derive(Clone)] -pub struct Problem<'a, const D: usize, S, P> +pub struct Problem<'a, T: Scalar, const D: usize, S, P, C> where S: Integrator, + C: Controller, { - ode: ODE<'a, D, P>, + ode: ODE<'a, T, D, P>, integrator: S, - controller: PIController, - callbacks: Vec>, + controller: C, + callbacks: Vec<&'a dyn Callback>, } -impl<'a, const D: usize, S, P> Problem<'a, D, S, P> +impl<'a, T: Scalar, const D: usize, S, P, C> Problem<'a, T, D, S, P, C> where S: Integrator + Copy, + C: Controller, { - pub fn new(ode: ODE<'a, D, P>, integrator: S, controller: PIController) -> Self { + pub fn new(ode: ODE<'a, T, D, P>, integrator: S, controller: C) -> Self { Problem { ode, integrator, @@ -30,69 +35,79 @@ where callbacks: Vec::new(), } } - pub fn solve(&mut self) -> Solution<'_, S, D, P> { + + pub fn solve(&mut self) -> Solution<'_, T, S, D, P> { let mut convergency = SimpleConvergency { - eps: 1e-12, + eps: 1e-12f64, max_iter: 1000, }; - let mut times: Vec = vec![self.ode.t]; - let mut states: Vec> = vec![self.ode.y]; - let mut dense_coefficients: Vec>>> = Vec::new(); - while self.ode.t < self.ode.t_end { - if self.ode.t + self.controller.next_step_guess.extract() > self.ode.t_end { + let mut times: Vec = vec![self.ode.t]; + let mut states: Vec> = vec![self.ode.y]; + let mut dense_coefficients: Vec>>> = Vec::new(); + while self.ode.t.re() < self.ode.t_end.re() { + if (self.ode.t + self.controller.next_step_guess().extract()).re() + > self.ode.t_end.re() + { // If the next step would go past the end, then just set it to the end - self.controller.next_step_guess = TryStep::NotYetAccepted( - self.ode.t_end - self.ode.t, - ); + self.controller + .set_next_step_guess(TryStep::NotYetAccepted(self.ode.t_end - self.ode.t)); } let (mut new_y, mut curr_step, mut dense_option) = if S::ADAPTIVE { // First, we try stepping with the "next step guess" to get the error - let (mut trial_y, mut err_option, mut dense_option) = - self.integrator.step(&self.ode, self.controller.next_step_guess.extract()); - let mut err = err_option.unwrap(); + let mut result = self + .integrator + .step(&self.ode, self.controller.next_step_guess().extract()); + let mut err = result.err.unwrap(); // Then we determine whether we need to reduce the step size or not // If successful, we get the next step guess - let initial_guess = self.controller.next_step_guess.extract(); - let mut next_step_guess = >::determine_step( - &mut self.controller, - initial_guess, - err, - ); + let initial_guess = self.controller.next_step_guess().extract(); + let mut next_step_guess = self.controller.determine_step(initial_guess, err); while !next_step_guess.is_accepted() { // If that step isn't acceptable, then change the step until it is - (trial_y, err_option, dense_option) = - self.integrator.step(&self.ode, next_step_guess.extract()); - next_step_guess = >::determine_step( - &mut self.controller, - next_step_guess.extract(), - err, - ); - err = err_option.unwrap(); + result = self.integrator.step(&self.ode, next_step_guess.extract()); + err = result.err.unwrap(); + next_step_guess = self + .controller + .determine_step(next_step_guess.extract(), err); } // So at this point we can safely assume we have an accepted step - self.controller.next_step_guess = next_step_guess.reset().unwrap(); - (trial_y, next_step_guess.extract(), dense_option) + self.controller + .set_next_step_guess(next_step_guess.reset().unwrap()); + (result.y, next_step_guess.extract(), result.dense) } else { // If fixed time step just step forward one step - let (trial_y, _, dense_option) = self.integrator.step(&self.ode, self.controller.next_step_guess.extract()); - (trial_y, self.controller.next_step_guess.extract(), dense_option) + let result = self + .integrator + .step(&self.ode, self.controller.next_step_guess().extract()); + ( + result.y, + self.controller.next_step_guess().extract(), + result.dense, + ) }; if !self.callbacks.is_empty() { // Check for events occurring for callback in &self.callbacks { - if (callback.event)(self.ode.t, self.ode.y, &self.ode.params) - * (callback.event)(self.ode.t + curr_step, new_y, &self.ode.params) - < 0.0 - { + // The event condition is a real-valued question, so it is + // asked in the reals — see this module's header. + let before = callback.event(self.ode.t, &self.ode.y, &self.ode.params).re(); + let after = callback + .event(self.ode.t + curr_step, &new_y, &self.ode.params) + .re(); + if before * after < 0.0 { // If the event crossed zero, then find the root - let f = |test_t| { - let test_y = self.integrator.step(&self.ode, test_t).0; - (callback.event)(self.ode.t + test_t, test_y, &self.ode.params) + let f = |test_h: f64| { + let h = T::lift(test_h); + let test_y = self.integrator.step(&self.ode, h).y; + callback.event(self.ode.t + h, &test_y, &self.ode.params).re() }; - let root = find_root_brent(0.0, curr_step, &f, &mut convergency).unwrap(); - curr_step = root; - (new_y, _, dense_option) = self.integrator.step(&self.ode, curr_step); - (callback.effect)(&mut self.ode); + let root = + find_root_brent(0.0f64, curr_step.re(), &f, &mut convergency).unwrap(); + curr_step = T::lift(root); + let result = self.integrator.step(&self.ode, curr_step); + new_y = result.y; + dense_option = result.dense; + callback.effect(&mut self.ode); } } } @@ -112,52 +127,48 @@ where } } - pub fn with_callback(mut self, callback: Callback<'a, D, P>) -> Self { + pub fn with_callback(mut self, callback: &'a dyn Callback) -> Self { self.callbacks.push(callback); - Self { - ode: self.ode, - integrator: self.integrator, - controller: self.controller, - callbacks: self.callbacks, - } + self } } -pub struct Solution<'a, S, const D: usize, P> +pub struct Solution<'a, T: Scalar, S, const D: usize, P> where S: Integrator, { - pub ode: &'a ODE<'a, D, P>, + pub ode: &'a ODE<'a, T, D, P>, pub integrator: S, - pub times: Vec, - pub states: Vec>, - pub dense: Vec>>>, + pub times: Vec, + pub states: Vec>, + pub dense: Vec>>>, } -impl<'a, S, const D: usize, P> Solution<'a, S, D, P> +impl Solution<'_, T, S, D, P> where S: Integrator, { - pub fn interpolate(&self, t: f64) -> SVector { + pub fn interpolate(&self, t: T) -> SVector { // First check that the t is within bounds - let last = self.times.last().unwrap(); - let first = self.times.first().unwrap(); + let last = self.times.last().unwrap().re(); + let first = self.times.first().unwrap().re(); // TODO: Improve these errors - let mut times = self.times.clone(); - if *first > *last { + let mut times: Vec = self.times.iter().map(|x| x.re()).collect(); + if first > last { times.reverse(); } - if t < *first || t > *last { + if t.re() < first || t.re() > last { panic!(); } // Then find the two t values closest to the desired t - match times.binary_search_by(|x| x.total_cmp(&t)) { + match times.binary_search_by(|x| x.total_cmp(&t.re())) { Ok(index) => self.states[index], Err(end_index) => { - let t_start = times[end_index - 1]; - let t_end = times[end_index]; + // Then send that to the integrator + let t_start = self.times[end_index - 1]; + let t_end = self.times[end_index]; let y_start = self.states[end_index - 1]; let h = t_end - t_start; @@ -205,7 +216,7 @@ where #[cfg(test)] mod tests { use super::*; - use crate::callback::stop; + use crate::callback::{stop, FnCallback}; use crate::controller::PIController; use crate::integrator::dormand_prince::DormandPrince45; use approx::assert_relative_eq; @@ -247,12 +258,12 @@ mod tests { let dp45 = DormandPrince45::new().a_tol(1e-12).r_tol(1e-5); let controller = PIController::default(); - let value_too_high = Callback { - event: &|_: f64, y: SVector, _: &Params| 10.0 - y[0], + let value_too_high = FnCallback { + event: &|_: f64, y: &SVector, _: &Params| 10.0 - y[0], effect: &stop, }; - let mut problem = Problem::new(ode, dp45, controller).with_callback(value_too_high); + let mut problem = Problem::new(ode, dp45, controller).with_callback(&value_too_high); let solution = problem.solve(); assert_relative_eq!( diff --git a/src/scalar.rs b/src/scalar.rs new file mode 100644 index 0000000..6cade67 --- /dev/null +++ b/src/scalar.rs @@ -0,0 +1,119 @@ +//! The scalar type an ODE is integrated in. +//! +//! All of our problem types: state, time and step size, is generic over +//! [`Scalar`]. Instantiating a problem at `f64` solves it. Instantiating the +//! same problem at a forward-mode dual number solves it *and* propagates exact +//! derivatives alongside the state, at machine precision, with no finite +//! differencing and no extra right-hand-side evaluations beyond the ones the +//! dual arithmetic performs anyway. +//! +//! [`Scalar`] is implemented for `f64` always, and for the dual numbers of the +//! [`num_dual`] crate under the `autodiff` feature (on by default). + +use nalgebra::RealField; +use num_traits::NumOps; + +/// A scalar that an ODE's state, time and step size may take. +pub trait Scalar: RealField + Copy + NumOps { + /// This value with all derivative information discarded. + fn re(&self) -> f64; + + /// A real constant as a value of this type, with zero derivative. + fn lift(x: f64) -> Self; +} + +impl Scalar for f64 { + #[inline(always)] + fn re(&self) -> f64 { + *self + } + + #[inline(always)] + fn lift(x: f64) -> Self { + x + } +} + +#[cfg(feature = "autodiff")] +mod dual { + use super::Scalar; + use nalgebra::{Const, RealField}; + use num_dual::{DualNum, DualStruct}; + + /// Implements `Scalar` for one of `num-dual`'s dual number types. + macro_rules! impl_scalar { + ($ty:ident) => { + impl Scalar for num_dual::$ty + where + Self: RealField + Copy + DualNum + DualStruct, + { + #[inline(always)] + fn re(&self) -> f64 { + DualStruct::re(self) + } + + #[inline(always)] + fn lift(x: f64) -> Self { + Self::from(x) + } + } + }; + ($ty:ident, $($n:ident),+) => { + impl Scalar + for num_dual::$ty),+> + where + Self: RealField + Copy + DualNum + DualStruct, + { + #[inline(always)] + fn re(&self) -> f64 { + DualStruct::re(self) + } + + #[inline(always)] + fn lift(x: f64) -> Self { + Self::from(x) + } + } + }; + } + + impl_scalar!(Dual); + impl_scalar!(Dual2); + impl_scalar!(Dual3); + impl_scalar!(HyperDual); + impl_scalar!(HyperHyperDual); + impl_scalar!(DualVec, N); + impl_scalar!(Dual2Vec, N); + impl_scalar!(HyperDualVec, M, N); +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_f64_is_a_scalar() { + assert_eq!(f64::lift(3.5), 3.5); + assert_eq!(Scalar::re(&3.5), 3.5); + } + + #[cfg(feature = "autodiff")] + #[test] + fn test_dual_projects_to_its_real_part() { + use num_dual::{Dual64, DualNum}; + + // A lifted constant has no derivative; a seeded variable does, and both + // project back to the same real value. + let c = ::lift(3.5); + assert_eq!(c.re, 3.5); + assert_eq!(c.eps, 0.0); + + let x = Dual64::from_re(3.5).derivative(); + assert_eq!(Scalar::re(&x), 3.5); + assert_eq!(x.eps, 1.0); + + // d/dx x^3 = 3x^2 = 36.75 at x = 3.5 + assert_eq!(x.powi(3).eps, 3.0 * 3.5 * 3.5); + } + +}