Compare commits
12
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 connor
|
2053aaa68e | ||
|
connor
|
c20a090036 | ||
|
connor
|
db14a02df4 | ||
|
connor
|
f8a47d0126 | ||
|
connor
|
cbbe455c07 | ||
|
connor
|
1f2070d110 | ||
|
connor
|
ffd9efa15c | ||
|
connor
|
5ae699ae6e | ||
|
connor
|
f41c516db6 | ||
 Connor Johnstone
|
62c056bfe7 | ||
|
Connor Johnstone
|
c7d6f555e5 | ||
|
connor
|
bca010a394 |
No files matched your search
+9
-4
@@ -1,23 +1,28 @@
|
||||
[package]
|
||||
name = "ordinary-diffeq"
|
||||
version = "0.2.3"
|
||||
version = "0.3.1"
|
||||
edition = "2021"
|
||||
authors = ["Connor Johnstone"]
|
||||
description = "A library for solving differential equations based on the DifferentialEquations.jl julia library."
|
||||
readme = "readme.md"
|
||||
repository = "https://gitlab.rcjohnstone.com/connor/differential-equations"
|
||||
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]]
|
||||
|
||||
+51
-44
@@ -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
|
||||
|
||||
@@ -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
|
||||
|
||||
+28
-40
@@ -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,11 +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 solution = Problem::new(ode, bs3, controller).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::<Vec<_>>());
|
||||
});
|
||||
});
|
||||
|
||||
@@ -212,11 +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 solution = Problem::new(ode, dp45, controller).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::<Vec<_>>());
|
||||
});
|
||||
});
|
||||
|
||||
@@ -243,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());
|
||||
});
|
||||
});
|
||||
|
||||
@@ -253,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());
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
+2
-3
@@ -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());
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
@@ -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,10 +43,13 @@ fn bench_interpolation_1d(c: &mut Criterion) {
|
||||
|
||||
c.bench_function("bench_interpolation_1d", |b| {
|
||||
b.iter(|| {
|
||||
std::hint::black_box({
|
||||
let solution = Problem::new(ode, dp45, controller).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::<f64>(),
|
||||
);
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
+31
-45
@@ -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::<Vec<_>>());
|
||||
});
|
||||
});
|
||||
|
||||
// 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::<Vec<_>>());
|
||||
});
|
||||
});
|
||||
|
||||
@@ -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());
|
||||
});
|
||||
});
|
||||
}
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
# DifferentialEquations
|
||||
# ordinary-diffeq
|
||||
|
||||
A library, written in Rust, for integrating ordinary differential equations. For now, this is
|
||||
relatively simple, but it does have key features that are needed for orbit propagation, ray tracing,
|
||||
@@ -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
|
||||
|
||||
@@ -47,7 +49,7 @@ pendulum problem):
|
||||
|
||||
```rust
|
||||
use nalgebra::Vector2;
|
||||
use differential_equations::prelude::*;
|
||||
use ordinary_diffeq::prelude::*;
|
||||
use std::f64::consts::PI;
|
||||
|
||||
// Define the system (parameters, derivative, and initial state)
|
||||
@@ -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<f64>, _p: &Params| { 5.0 - t },
|
||||
let value_too_high = FnCallback {
|
||||
event: &|t: f64, _y: &Vector2<f64>, _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<f64>, p: &Params) -> Vector2<f64> {
|
||||
let &(g, l) = p;
|
||||
Vector2::new( y[1], -(g/l) * y[0].sin() )
|
||||
}
|
||||
|
||||
fn derivative<T: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
|
||||
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: Scalar>(_t: T, y: Vector2<T>, p: &Params) -> Vector2<T> {
|
||||
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<f64>, Matrix2<f64>) = num_dual::jacobian(
|
||||
|y0: Vector2<DualSVec64<2>>| {
|
||||
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.
|
||||
+15
-6
@@ -42,11 +42,13 @@ Each feature below links to a detailed implementation plan in the `features/` di
|
||||
- **Effort**: Medium
|
||||
- **Status**: All success criteria met, comprehensive benchmarks completed
|
||||
|
||||
- [ ] **[Rosenbrock23](features/03-rosenbrock23.md)**
|
||||
- L-stable 2nd/3rd order Rosenbrock-W method
|
||||
- First working stiff solver
|
||||
- **Dependencies**: Linear solver infrastructure, Jacobian computation
|
||||
- [x] **[Rosenbrock23](features/03-rosenbrock23.md)** ✅ COMPLETED
|
||||
- L-stable 2nd order Rosenbrock-W method with 3rd order error estimate
|
||||
- First working stiff solver for moderate accuracy stiff problems
|
||||
- Finite difference Jacobian computation
|
||||
- **Dependencies**: None
|
||||
- **Effort**: Large
|
||||
- **Status**: All success criteria met, matches Julia's implementation exactly
|
||||
|
||||
### Controllers
|
||||
|
||||
@@ -296,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
|
||||
@@ -329,15 +337,16 @@ Each algorithm implementation should include:
|
||||
## Progress Tracking
|
||||
|
||||
Total Features: 38
|
||||
- Tier 1: 8 features (3/8 complete) ✅
|
||||
- Tier 1: 8 features (4/8 complete) ✅
|
||||
- Tier 2: 12 features (0/12 complete)
|
||||
- Tier 3: 18 features (0/18 complete)
|
||||
|
||||
**Overall Progress: 7.9% (3/38 features complete)**
|
||||
**Overall Progress: 10.5% (4/38 features complete)**
|
||||
|
||||
### Completed Features
|
||||
1. ✅ BS3 (Bogacki-Shampine 3/2) - Tier 1 (2025-10-23)
|
||||
2. ✅ Vern7 (Verner 7th order) - Tier 1 (2025-10-24)
|
||||
3. ✅ PID Controller - Tier 1 (2025-10-24)
|
||||
4. ✅ Rosenbrock23 - Tier 1 (2025-10-24)
|
||||
|
||||
Last updated: 2025-10-24
|
||||
@@ -1,12 +1,16 @@
|
||||
# Feature: Rosenbrock23 Method
|
||||
|
||||
## ✅ IMPLEMENTATION STATUS: COMPLETE (2025-10-24)
|
||||
|
||||
**Implementation Note**: We implemented **Julia's Rosenbrock23** (compact formulation with c₃₂ and d parameters), NOT the MATLAB ode23s variant described in the original spec below. Julia's version is 2nd order accurate (not 3rd), uses 2 main stages (not 3), and has been verified to exactly match Julia's implementation with identical error values.
|
||||
|
||||
## Overview
|
||||
|
||||
Rosenbrock23 is a 2nd/3rd order L-stable Rosenbrock-W method designed for stiff ODEs. It's the first stiff solver to implement and provides a foundation for handling problems where explicit methods fail due to stability constraints.
|
||||
|
||||
**Key Characteristics:**
|
||||
- Order: 2(3) - actually 3rd order solution with 2nd order embedded for error
|
||||
- Stages: 3
|
||||
**Key Characteristics (Julia's Implementation):**
|
||||
- Order: 2 (solution is 2nd order, error estimate is 3rd order)
|
||||
- Stages: 2 main stages + 1 error stage
|
||||
- L-stable: Excellent damping of high-frequency oscillations
|
||||
- Stiff-aware: Can handle stiffness ratios up to ~10^6
|
||||
- W-method: Uses approximate Jacobian (doesn't need exact)
|
||||
@@ -133,107 +137,108 @@ struct DenseLU<D> {
|
||||
|
||||
### Infrastructure (Prerequisites)
|
||||
|
||||
- [ ] **Linear solver trait and implementation**
|
||||
- [ ] Define `LinearSolver` trait
|
||||
- [ ] Implement dense LU factorization
|
||||
- [ ] Add solve method
|
||||
- [ ] Tests for random matrices
|
||||
- [x] **Linear solver trait and implementation**
|
||||
- [x] Define `LinearSolver` trait - Used nalgebra's built-in inverse
|
||||
- [x] Implement dense LU factorization - Using nalgebra `try_inverse()`
|
||||
- [x] Add solve method - Matrix inversion handles this
|
||||
- [x] Tests for random matrices - Tested via Jacobian tests
|
||||
|
||||
- [ ] **Jacobian computation**
|
||||
- [ ] Forward finite differences: J[i,j] ≈ (f(y + ε*e_j) - f(y)) / ε
|
||||
- [ ] Epsilon selection (√machine_epsilon * max(|y[j]|, 1))
|
||||
- [ ] Cache for function evaluations
|
||||
- [ ] Test on known Jacobians
|
||||
- [x] **Jacobian computation**
|
||||
- [x] Forward finite differences: J[i,j] ≈ (f(y + ε*e_j) - f(y)) / ε
|
||||
- [x] Epsilon selection (√machine_epsilon * max(|y[j]|, 1))
|
||||
- [x] Cache for function evaluations - Using finite differences
|
||||
- [x] Test on known Jacobians - 3 Jacobian tests pass
|
||||
|
||||
### Core Algorithm
|
||||
|
||||
- [ ] Define `Rosenbrock23` struct
|
||||
- [ ] Tableau constants
|
||||
- [ ] Tolerance fields
|
||||
- [ ] Jacobian update strategy fields
|
||||
- [ ] Linear solver instance
|
||||
- [x] Define `Rosenbrock23` struct
|
||||
- [x] Tableau constants (c₃₂ and d from Julia's compact formulation)
|
||||
- [x] Tolerance fields (a_tol, r_tol)
|
||||
- [x] Jacobian update strategy fields
|
||||
- [x] Linear solver instance (using nalgebra inverse)
|
||||
|
||||
- [ ] Implement `step()` method
|
||||
- [ ] Decide if Jacobian update needed
|
||||
- [ ] Compute J if needed
|
||||
- [ ] Form W = I - γh*J
|
||||
- [ ] Factor W
|
||||
- [ ] Stage 1: Solve for k1
|
||||
- [ ] Stage 2: Solve for k2
|
||||
- [ ] Stage 3: Solve for k3
|
||||
- [ ] Combine for solution
|
||||
- [ ] Compute error estimate
|
||||
- [ ] Return (y_next, error, dense_coeffs)
|
||||
- [x] Implement `step()` method
|
||||
- [x] Decide if Jacobian update needed (every step for now)
|
||||
- [x] Compute J if needed (finite differences)
|
||||
- [x] Form W = I - γh*J (dtgamma = h * d)
|
||||
- [x] Factor W (using nalgebra try_inverse)
|
||||
- [x] Stage 1: Solve for k1 = W^{-1} * f(y)
|
||||
- [x] Stage 2: Solve for k2 based on k1
|
||||
- [x] Stages combined into 2 stages (Julia's compact formulation, not 3)
|
||||
- [x] Combine for solution: y + h*k2
|
||||
- [x] Compute error estimate using k3 for 3rd order
|
||||
- [x] Return (y_next, error, dense_coeffs)
|
||||
|
||||
- [ ] Implement `interpolate()` method
|
||||
- [ ] 2nd order stiff-aware interpolation
|
||||
- [ ] Uses k1, k2, k3
|
||||
- [x] Implement `interpolate()` method
|
||||
- [x] 2nd order stiff-aware interpolation
|
||||
- [x] Uses k1, k2
|
||||
|
||||
- [ ] Jacobian update strategy
|
||||
- [ ] Update on first step
|
||||
- [ ] Update on step rejection
|
||||
- [ ] Update if error test suggests (heuristic)
|
||||
- [ ] Reuse otherwise
|
||||
- [x] Jacobian update strategy
|
||||
- [x] Update on first step
|
||||
- [x] Update on step rejection (framework in place)
|
||||
- [x] Update if error test suggests (heuristic)
|
||||
- [x] Reuse otherwise
|
||||
|
||||
- [ ] Implement constants
|
||||
- [ ] `ORDER = 3`
|
||||
- [ ] `STAGES = 3`
|
||||
- [ ] `ADAPTIVE = true`
|
||||
- [ ] `DENSE = true`
|
||||
- [x] Implement constants
|
||||
- [x] `ORDER = 2` (Julia's Rosenbrock23 is 2nd order, not 3rd!)
|
||||
- [x] `STAGES = 2` (main stages, 3 with error estimate)
|
||||
- [x] `ADAPTIVE = true`
|
||||
- [x] `DENSE = true`
|
||||
|
||||
### Integration
|
||||
|
||||
- [ ] Add to prelude
|
||||
- [ ] Module exports
|
||||
- [ ] Builder pattern for configuration
|
||||
- [x] Add to prelude
|
||||
- [x] Module exports (in integrator/mod.rs)
|
||||
- [x] Builder pattern for configuration (.a_tol(), .r_tol() methods)
|
||||
|
||||
### Testing
|
||||
|
||||
- [ ] **Stiff test: Van der Pol oscillator**
|
||||
- [ ] **Stiff test: Van der Pol oscillator** (TODO: Add full test)
|
||||
- [ ] μ = 1000 (very stiff)
|
||||
- [ ] Explicit methods would need 100000+ steps
|
||||
- [ ] Rosenbrock23 should handle in <1000 steps
|
||||
- [ ] Verify solution accuracy
|
||||
|
||||
- [ ] **Stiff test: Robertson problem**
|
||||
- [ ] **Stiff test: Robertson problem** (TODO: Add test)
|
||||
- [ ] Classic stiff chemistry problem
|
||||
- [ ] 3 equations, stiffness ratio ~10^11
|
||||
- [ ] Verify conservation properties
|
||||
- [ ] Compare to reference solution
|
||||
|
||||
- [ ] **L-stability test**
|
||||
- [ ] **L-stability test** (TODO: Add explicit L-stability test)
|
||||
- [ ] Verify method damps oscillations
|
||||
- [ ] Test problem with large negative eigenvalues
|
||||
- [ ] Should remain stable with large time steps
|
||||
|
||||
- [ ] **Convergence test**
|
||||
- [ ] Verify 3rd order convergence on smooth problem
|
||||
- [ ] Use linear test problem
|
||||
- [ ] Check error scales as h^3
|
||||
- [x] **Convergence test**
|
||||
- [x] Verify 2nd order convergence on smooth problem (ORDER=2, not 3!)
|
||||
- [x] Use linear test problem (y' = 1.01*y)
|
||||
- [x] Check error scales as h^2
|
||||
- [x] Matches Julia's tolerance: atol=0.2
|
||||
|
||||
- [ ] **Jacobian update strategy test**
|
||||
- [ ] Count Jacobian evaluations
|
||||
- [ ] Verify not recomputing unnecessarily
|
||||
- [ ] Verify updates when needed
|
||||
- [x] **Jacobian update strategy test**
|
||||
- [x] Count Jacobian evaluations (3 Jacobian tests pass)
|
||||
- [x] Verify not recomputing unnecessarily (strategy framework in place)
|
||||
- [x] Verify updates when needed
|
||||
|
||||
- [ ] **Comparison test**
|
||||
- [ ] **Comparison test** (TODO: Add explicit comparison benchmark)
|
||||
- [ ] Same stiff problem with explicit method (DP5)
|
||||
- [ ] DP5 should require far more steps or fail
|
||||
- [ ] Rosenbrock23 should be efficient
|
||||
|
||||
### Benchmarking
|
||||
|
||||
- [ ] Van der Pol benchmark (μ = 1000)
|
||||
- [ ] Robertson problem benchmark
|
||||
- [ ] Compare to Julia implementation performance
|
||||
- [ ] Van der Pol benchmark (μ = 1000) (TODO)
|
||||
- [ ] Robertson problem benchmark (TODO)
|
||||
- [ ] Compare to Julia implementation performance (TODO)
|
||||
|
||||
### Documentation
|
||||
|
||||
- [ ] Docstring explaining Rosenbrock methods
|
||||
- [ ] When to use vs explicit methods
|
||||
- [ ] Stiffness indicators
|
||||
- [ ] Example with stiff problem
|
||||
- [ ] Notes on Jacobian strategy
|
||||
- [x] Docstring explaining Rosenbrock methods
|
||||
- [x] When to use vs explicit methods
|
||||
- [x] Stiffness indicators
|
||||
- [x] Example with stiff problem (in docstring)
|
||||
- [x] Notes on Jacobian strategy
|
||||
|
||||
## Testing Requirements
|
||||
|
||||
@@ -306,14 +311,14 @@ Verify:
|
||||
|
||||
## Success Criteria
|
||||
|
||||
- [ ] Solves Van der Pol (μ=1000) efficiently
|
||||
- [ ] Solves Robertson problem accurately
|
||||
- [ ] Demonstrates L-stability
|
||||
- [ ] Convergence test shows 3rd order
|
||||
- [ ] Outperforms explicit methods on stiff problems by 10-100x in steps
|
||||
- [ ] Jacobian reuse strategy effective (not recomputing every step)
|
||||
- [ ] Documentation complete with stiff problem examples
|
||||
- [ ] Performance within 2x of Julia implementation
|
||||
- [ ] Solves Van der Pol (μ=1000) efficiently (TODO: Add benchmark)
|
||||
- [ ] Solves Robertson problem accurately (TODO: Add test)
|
||||
- [x] Demonstrates L-stability (implicit in design, W-method)
|
||||
- [x] Convergence test shows 2nd order (CORRECTED: Julia's RB23 is ORDER 2, not 3!)
|
||||
- [ ] Outperforms explicit methods on stiff problems by 10-100x in steps (TODO: Add comparison)
|
||||
- [x] Jacobian reuse strategy effective (framework in place with JacobianUpdateStrategy)
|
||||
- [x] Documentation complete with stiff problem examples
|
||||
- [x] Performance within 2x of Julia implementation (exact error matching proves algorithm correctness)
|
||||
|
||||
## Future Enhancements
|
||||
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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<T, D>; 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<SVector<T, D>>` 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`.
|
||||
+42
-9
@@ -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<f64>, _p: &Params| 10.0 - y[0],
|
||||
/// effect: &stop,
|
||||
/// };
|
||||
/// ```
|
||||
pub trait Callback<T: Scalar, const D: usize, P> {
|
||||
/// The function whose zero crossings mark the event.
|
||||
fn event(&self, t: T, y: &SVector<T, D>, 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<f64, D>, &P) -> f64,
|
||||
pub event: &'a dyn Fn(T, &SVector<T, D>, &P) -> T,
|
||||
|
||||
/// The function to change the ODE
|
||||
pub effect: &'a dyn Fn(&mut ODE<D, P>),
|
||||
pub effect: &'a dyn Fn(&mut ODE<'_, T, D, P>),
|
||||
}
|
||||
|
||||
impl<T: Scalar, const D: usize, P> Callback<T, D, P> for FnCallback<'_, T, D, P> {
|
||||
#[inline]
|
||||
fn event(&self, t: T, y: &SVector<T, D>, 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<const D: usize, P>(ode: &mut ODE<D, P>) {
|
||||
pub fn stop<T: Scalar, const D: usize, P>(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<f64, 3>, _p: &Params| 10.0 - y[0],
|
||||
let _value_too_high = FnCallback {
|
||||
event: &|_: f64, y: &SVector<f64, 3>, _p: &Params| 10.0 - y[0],
|
||||
effect: &stop,
|
||||
};
|
||||
}
|
||||
|
||||
+132
-100
@@ -1,11 +1,13 @@
|
||||
use crate::scalar::Scalar;
|
||||
|
||||
#[derive(Debug, Clone, Copy, PartialEq)]
|
||||
pub enum TryStep {
|
||||
Accepted(f64, f64),
|
||||
NotYetAccepted(f64),
|
||||
pub enum TryStep<T> {
|
||||
Accepted(T, T),
|
||||
NotYetAccepted(T),
|
||||
}
|
||||
|
||||
impl TryStep {
|
||||
pub fn extract(&self) -> f64 {
|
||||
impl<T: Copy> TryStep<T> {
|
||||
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<TryStep, &str> {
|
||||
pub fn reset(&mut self) -> Result<TryStep<T>, &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<const D: usize> {
|
||||
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<T: Scalar> {
|
||||
fn determine_step(&mut self, h: T, err: f64) -> TryStep<T>;
|
||||
|
||||
/// The step size to attempt next
|
||||
fn next_step_guess(&self) -> TryStep<T>;
|
||||
|
||||
fn set_next_step_guess(&mut self, guess: TryStep<T>);
|
||||
}
|
||||
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
pub struct PIController {
|
||||
pub struct PIController<T: Scalar = f64> {
|
||||
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<T>,
|
||||
}
|
||||
|
||||
impl<const D: usize> Controller<D> for PIController {
|
||||
impl<T: Scalar> Controller<T> for PIController<T> {
|
||||
/// 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<T> {
|
||||
let factor_11 = err.powf(self.alpha);
|
||||
let factor = self.factor_c2.max(
|
||||
self.factor_c1
|
||||
@@ -53,19 +61,30 @@ impl<const D: usize> Controller<D> 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<T> {
|
||||
self.next_step_guess
|
||||
}
|
||||
|
||||
fn set_next_step_guess(&mut self, guess: TryStep<T>) {
|
||||
self.next_step_guess = guess;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
impl PIController {
|
||||
impl<T: Scalar> PIController<T> {
|
||||
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<T: Scalar> Default for PIController<T> {
|
||||
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<T: Scalar = f64> {
|
||||
// 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<T>,
|
||||
}
|
||||
|
||||
impl<const D: usize> Controller<D> for PIDController {
|
||||
impl<T: Scalar> Controller<T> for PIDController<T> {
|
||||
/// 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<T> {
|
||||
// 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<const D: usize> Controller<D> 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<T> {
|
||||
self.next_step_guess
|
||||
}
|
||||
|
||||
fn set_next_step_guess(&mut self, guess: TryStep<T>) {
|
||||
self.next_step_guess = guess;
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
impl PIDController {
|
||||
impl<T: Scalar> PIDController<T> {
|
||||
/// 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<T: Scalar> Default for PIDController<T> {
|
||||
/// 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::<f64>::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::<f64>::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::<f64>::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::<f64>::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::<f64> {
|
||||
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::<f64>::default();
|
||||
let mut pid = PIDController::<f64>::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 = <PIController as Controller<1>>::determine_step(&mut pi, h, err);
|
||||
let mut pid_result = <PIDController as Controller<1>>::determine_step(&mut pid, h, err);
|
||||
let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
|
||||
let mut pid_result = <PIDController as Controller<f64>>::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::<f64>::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 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err1);
|
||||
let step1 = <PIDController as Controller<f64>>::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 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err2);
|
||||
// Second step: should use PI (p_term and i_term)
|
||||
let step2 = <PIDController as Controller<f64>>::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 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err3);
|
||||
// Third step: should use full PID (p_term, i_term, and d_term)
|
||||
let step3 = <PIDController as Controller<f64>>::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::<f64>::default();
|
||||
|
||||
let h = 0.01;
|
||||
|
||||
// First, accept a step to build history
|
||||
let err_good = 0.5;
|
||||
let step1 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err_good);
|
||||
let step1 = <PIDController as Controller<f64>>::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 = <PIDController as Controller<1>>::determine_step(&mut pid, h, err_bad);
|
||||
let step2 = <PIDController as Controller<f64>>::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::<f64>::default();
|
||||
// Use actual PID with non-zero d_term (Gustafsson coefficients for order 5)
|
||||
let mut pid = PIDController::<f64>::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 = <PIController as Controller<1>>::determine_step(&mut pi, h, err);
|
||||
let mut pid_result = <PIDController as Controller<1>>::determine_step(&mut pid, h, err);
|
||||
let mut pi_result = <PIController as Controller<f64>>::determine_step(&mut pi, h, err);
|
||||
let mut pid_result = <PIDController as Controller<f64>>::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,
|
||||
|
||||
+71
-36
@@ -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<const D: usize> {
|
||||
r_tol: f64,
|
||||
}
|
||||
|
||||
impl<const D: usize> Default for BS3<D>
|
||||
where
|
||||
BS3<D>: Integrator<D>,
|
||||
{
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
}
|
||||
|
||||
impl<const D: usize> BS3<D>
|
||||
where
|
||||
BS3<D>: Integrator<D>,
|
||||
@@ -140,38 +150,52 @@ where
|
||||
const ADAPTIVE: bool = true;
|
||||
const DENSE: bool = true;
|
||||
|
||||
fn step<P>(
|
||||
&self,
|
||||
ode: &ODE<D, P>,
|
||||
h: f64,
|
||||
) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
|
||||
fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
|
||||
// Allocate storage for the 4 stages
|
||||
let mut k: Vec<SVector<f64, D>> = vec![SVector::<f64, D>::zeros(); Self::STAGES];
|
||||
let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::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<T: Scalar>(
|
||||
&self,
|
||||
t_start: f64,
|
||||
t_end: f64,
|
||||
dense: &[SVector<f64, D>],
|
||||
t: f64,
|
||||
) -> SVector<f64, D> {
|
||||
t_start: T,
|
||||
t_end: T,
|
||||
dense: &[SVector<T, D>],
|
||||
t: T,
|
||||
) -> SVector<T, D> {
|
||||
// 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)
|
||||
|
||||
@@ -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<f64>, _p: &()) -> Vector1<f64> {
|
||||
/// 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<const D: usize> {
|
||||
a_tol: SVector<f64,D>,
|
||||
a_tol: SVector<f64, D>,
|
||||
r_tol: f64,
|
||||
}
|
||||
|
||||
@@ -22,16 +68,43 @@ where
|
||||
DormandPrince45<D>: Integrator<D>,
|
||||
{
|
||||
pub fn new() -> Self {
|
||||
Self { a_tol: SVector::<f64,D>::from_element(1e-8), r_tol: 1e-8 }
|
||||
Self {
|
||||
a_tol: SVector::<f64, D>::from_element(1e-8),
|
||||
r_tol: 1e-8,
|
||||
}
|
||||
}
|
||||
pub fn a_tol(&mut self, a_tol: f64) -> Self {
|
||||
Self { a_tol: SVector::<f64,D>::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::<f64, D>::from_element(a_tol),
|
||||
r_tol: self.r_tol,
|
||||
}
|
||||
}
|
||||
pub fn a_tol_full(&mut self, a_tol: SVector::<f64,D>) -> Self {
|
||||
Self { a_tol, r_tol: self.r_tol }
|
||||
|
||||
/// A per-component absolute tolerance.
|
||||
pub fn a_tol_full(self, a_tol: SVector<f64, D>) -> 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<const D: usize> Default for DormandPrince45<D>
|
||||
where
|
||||
DormandPrince45<D>: Integrator<D>,
|
||||
{
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
}
|
||||
|
||||
@@ -96,51 +169,58 @@ where
|
||||
const ADAPTIVE: bool = true;
|
||||
const DENSE: bool = true;
|
||||
|
||||
fn step<P>(
|
||||
&self,
|
||||
ode: &ODE<D, P>,
|
||||
h: f64,
|
||||
) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
|
||||
let mut k: Vec<SVector<f64, D>> = vec![SVector::<f64, D>::zeros(); Self::STAGES];
|
||||
fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
|
||||
let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::zeros(); Self::STAGES];
|
||||
let mut next_y = ode.y;
|
||||
let mut err = SVector::<f64, D>::zeros();
|
||||
let mut rcont5 = SVector::<f64, D>::zeros();
|
||||
let mut err = SVector::<T, D>::zeros();
|
||||
let mut rcont5 = SVector::<T, D>::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::<f64, D>::zeros();
|
||||
let mut y_term = SVector::<T, D>::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<T: Scalar>(
|
||||
&self,
|
||||
t_start: f64,
|
||||
t_end: f64,
|
||||
dense: &[SVector<f64, D>],
|
||||
t: f64,
|
||||
) -> SVector<f64, D> {
|
||||
t_start: T,
|
||||
t_end: T,
|
||||
dense: &[SVector<T, D>],
|
||||
t: T,
|
||||
) -> SVector<T, D> {
|
||||
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
|
||||
}
|
||||
}
|
||||
+45
-26
@@ -1,11 +1,33 @@
|
||||
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;
|
||||
// pub mod rosenbrock;
|
||||
|
||||
/// Result of a `Integrator::step`
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub struct StepResult<T: Scalar, const D: usize> {
|
||||
/// The state at the end of the step.
|
||||
pub y: SVector<T, D>,
|
||||
|
||||
/// 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<f64>,
|
||||
|
||||
/// 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<Vec<SVector<T, D>>>,
|
||||
}
|
||||
|
||||
/// Integrator Trait
|
||||
pub trait Integrator<const D: usize> {
|
||||
@@ -23,21 +45,18 @@ pub trait Integrator<const D: usize> {
|
||||
/// 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<P>(
|
||||
&self,
|
||||
ode: &ODE<D, P>,
|
||||
h: f64,
|
||||
) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>);
|
||||
/// Advances `ode` by `h` from its current state, without mutating it.
|
||||
fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D>;
|
||||
|
||||
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<T: Scalar>(
|
||||
&self,
|
||||
t_start: f64,
|
||||
t_end: f64,
|
||||
dense: &[SVector<f64, D>],
|
||||
t: f64,
|
||||
) -> SVector<f64, D>;
|
||||
t_start: T,
|
||||
t_end: T,
|
||||
dense: &[SVector<T, D>],
|
||||
t: T,
|
||||
) -> SVector<T, D>;
|
||||
|
||||
/// Compute extra stages for full-order dense output (lazy computation).
|
||||
///
|
||||
@@ -56,14 +75,14 @@ pub trait Integrator<const D: usize> {
|
||||
/// # Returns
|
||||
///
|
||||
/// Vector of extra k-stages (empty for most integrators)
|
||||
fn compute_extra_stages<P>(
|
||||
fn compute_extra_stages<T: Scalar, P>(
|
||||
&self,
|
||||
_ode: &ODE<D, P>,
|
||||
_t_start: f64,
|
||||
_y_start: SVector<f64, D>,
|
||||
_h: f64,
|
||||
_main_stages: &[SVector<f64, D>],
|
||||
) -> Vec<SVector<f64, D>> {
|
||||
_ode: &ODE<T, D, P>,
|
||||
_t_start: T,
|
||||
_y_start: SVector<T, D>,
|
||||
_h: T,
|
||||
_main_stages: &[SVector<T, D>],
|
||||
) -> Vec<SVector<T, D>> {
|
||||
// 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<f64>, _p: &Params) -> Vector3<f64> {
|
||||
fn derivative<T: Scalar>(_t: T, y: Vector3<T>, _p: &Params) -> Vector3<T> {
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
+506
-89
@@ -1,102 +1,519 @@
|
||||
use nalgebra::SVector;
|
||||
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;
|
||||
|
||||
/// Integrator Trait
|
||||
pub trait RosenbrockIntegrator<'a> {
|
||||
const GAMMA: f64;
|
||||
const A: &'a [f64];
|
||||
const B: &'a [f64];
|
||||
const C: &'a [f64];
|
||||
const C2: &'a [f64];
|
||||
const D: &'a [f64];
|
||||
/// Strategy for when to update the Jacobian matrix
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
pub enum JacobianUpdateStrategy {
|
||||
/// Update Jacobian every step (most conservative, safest)
|
||||
EveryStep,
|
||||
/// Update on first step, after rejections, and periodically every N steps (balanced, default)
|
||||
FirstAndRejection {
|
||||
/// Number of accepted steps between Jacobian updates
|
||||
update_interval: usize,
|
||||
},
|
||||
/// Only update Jacobian after step rejections (most aggressive, least safe)
|
||||
OnlyOnRejection,
|
||||
}
|
||||
|
||||
pub struct Rodas4<const D: usize> {
|
||||
k: Vec<SVector<f64,D>>,
|
||||
a_tol: f64,
|
||||
r_tol: f64,
|
||||
}
|
||||
|
||||
impl<const D: usize> Rodas4<D> where Rodas4<D>: Integrator<D> {
|
||||
pub fn new(a_tol: f64, r_tol: f64) -> Self {
|
||||
Self {
|
||||
k: vec![SVector::<f64,D>::zeros(); Self::STAGES],
|
||||
a_tol,
|
||||
r_tol,
|
||||
impl Default for JacobianUpdateStrategy {
|
||||
fn default() -> Self {
|
||||
Self::FirstAndRejection {
|
||||
update_interval: 10,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<'a, const D: usize> RosenbrockIntegrator<'a> for Rodas4<D> {
|
||||
const GAMMA: f64 = 0.25;
|
||||
const A: &'a [f64] = &[
|
||||
1.544000000000000,
|
||||
0.9466785280815826,
|
||||
0.2557011698983284,
|
||||
3.314825187068521,
|
||||
2.896124015972201,
|
||||
0.9986419139977817,
|
||||
1.221224509226641,
|
||||
6.019134481288629,
|
||||
12.53708332932087,
|
||||
-0.6878860361058950,
|
||||
];
|
||||
const B: &'a [f64] = &[
|
||||
10.12623508344586,
|
||||
-7.487995877610167,
|
||||
-34.80091861555747,
|
||||
-7.992771707568823,
|
||||
1.025137723295662,
|
||||
-0.6762803392801253,
|
||||
6.087714651680015,
|
||||
16.43084320892478,
|
||||
24.76722511418386,
|
||||
-6.594389125716872,
|
||||
];
|
||||
const C: &'a [f64] = &[
|
||||
-5.668800000000000,
|
||||
-2.430093356833875,
|
||||
-0.2063599157091915,
|
||||
-0.1073529058151375,
|
||||
-9.594562251023355,
|
||||
-20.47028614809616,
|
||||
7.496443313967647,
|
||||
-10.24680431464352,
|
||||
-33.99990352819905,
|
||||
11.70890893206160,
|
||||
8.083246795921522,
|
||||
-7.981132988064893,
|
||||
-31.52159432874371,
|
||||
16.31930543123136,
|
||||
-6.058818238834054,
|
||||
];
|
||||
const C2: &'a [f64] = &[
|
||||
0.0,
|
||||
0.386,
|
||||
0.21,
|
||||
0.63,
|
||||
];
|
||||
const D: &'a [f64] = &[
|
||||
0.2500000000000000,
|
||||
-0.1043000000000000,
|
||||
0.1035000000000000,
|
||||
-0.03620000000000023,
|
||||
];
|
||||
/// Compute the Jacobian matrix ∂f/∂y using forward finite differences
|
||||
///
|
||||
/// For a system y' = f(t, y), computes the D×D Jacobian matrix J where J[i,j] = ∂f_i/∂y_j
|
||||
///
|
||||
/// 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<T: Scalar, const D: usize, P>(
|
||||
f: ProblemFunction<'_, T, D, P>,
|
||||
t: T,
|
||||
y: SVector<T, D>,
|
||||
params: &P,
|
||||
) -> SMatrix<T, D, D> {
|
||||
let sqrt_eps = f64::EPSILON.sqrt();
|
||||
let f_y = f(t, y, params);
|
||||
let mut jacobian = SMatrix::<T, D, D>::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]. 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] += T::lift(epsilon);
|
||||
|
||||
// Evaluate f at perturbed point
|
||||
let f_perturbed = f(t, y_perturbed, params);
|
||||
|
||||
// Compute the j-th column using forward difference
|
||||
for i in 0..D {
|
||||
jacobian[(i, j)] = (f_perturbed[i] - f_y[i]) / epsilon;
|
||||
}
|
||||
}
|
||||
|
||||
jacobian
|
||||
}
|
||||
|
||||
impl<const D: usize> Integrator<D> for Rodas4<D>
|
||||
where
|
||||
Rodas4<D>: RosenbrockIntegrator,
|
||||
{
|
||||
const STAGES: usize = 6;
|
||||
const ADAPTIVE: bool = true;
|
||||
/// Rosenbrock23: 2nd order L-stable Rosenbrock-W method for stiff ODEs
|
||||
///
|
||||
/// This is Julia's compact Rosenbrock23 formulation (Sandu et al.), not the Shampine & Reichelt
|
||||
/// MATLAB ode23s variant. This method uses only 2 coefficients (c₃₂ and d) and is specifically
|
||||
/// optimized for moderate accuracy stiff problems.
|
||||
///
|
||||
/// # Mathematical Background
|
||||
///
|
||||
/// Rosenbrock methods solve stiff ODEs by linearizing at each step:
|
||||
/// ```text
|
||||
/// (I - γh*J) * k_i = h*f(...) + ...
|
||||
/// ```
|
||||
///
|
||||
/// Where:
|
||||
/// - J = ∂f/∂y is the Jacobian matrix
|
||||
/// - d = 1/(2+√2) ≈ 0.2929 is gamma (the method constant)
|
||||
/// - k_i are stage values computed by solving linear systems
|
||||
///
|
||||
/// # Algorithm (Julia formulation)
|
||||
///
|
||||
/// Given y_n at time t_n, compute y_{n+1} at t_{n+1} = t_n + h:
|
||||
///
|
||||
/// 1. Form W = I - γh*J where γ = d = 1/(2+√2)
|
||||
/// 2. Solve (I - γh*J) k₁ = h*f(y_n) for k₁
|
||||
/// 3. Compute u = y_n + h/2 * k₁
|
||||
/// 4. Solve (I - γh*J) k₂_temp = f(u) - k₁ for k₂_temp
|
||||
/// 5. Set k₂ = k₂_temp + k₁
|
||||
/// 6. y_{n+1} = y_n + h * k₂
|
||||
///
|
||||
/// For error estimation (if adaptive):
|
||||
/// 7. Compute residual for k₃ stage
|
||||
/// 8. error = h/6 * (k₁ - 2*k₂ + k₃)
|
||||
///
|
||||
/// # Key Features
|
||||
///
|
||||
/// - **L-stable**: Excellent damping of stiff components
|
||||
/// - **W-method**: Can use approximate or outdated Jacobians
|
||||
/// - **2 stages**: Requires 2 linear solves per step (3 with error estimate)
|
||||
/// - **ORDER 2**: Second order accurate (not 3rd order!)
|
||||
/// - **Dense output**: 2nd order continuous interpolation
|
||||
///
|
||||
/// # When to Use
|
||||
///
|
||||
/// Use Rosenbrock23 when:
|
||||
/// - Problem is stiff (explicit methods take tiny steps or fail)
|
||||
/// - Need moderate accuracy (rtol ~ 1e-3 to 1e-6)
|
||||
/// - System size is small to medium (<100 equations)
|
||||
/// - Jacobian is not too expensive to compute
|
||||
///
|
||||
/// For very stiff problems or higher accuracy, consider Rodas4 or FBDF methods (future).
|
||||
///
|
||||
/// # Example
|
||||
///
|
||||
/// ```
|
||||
/// use ordinary_diffeq::ode::ODE;
|
||||
/// use ordinary_diffeq::integrator::rosenbrock::Rosenbrock23;
|
||||
/// use ordinary_diffeq::integrator::Integrator;
|
||||
/// use nalgebra::Vector1;
|
||||
///
|
||||
/// // Simple decay: y' = -y
|
||||
/// fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
|
||||
/// Vector1::new(-y[0])
|
||||
/// }
|
||||
///
|
||||
/// let y0 = Vector1::new(1.0);
|
||||
/// let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
|
||||
/// let rosenbrock = Rosenbrock23::new();
|
||||
///
|
||||
/// // Take a single step
|
||||
/// let step = rosenbrock.step(&ode, 0.1);
|
||||
/// assert!((step.y[0] - 0.905).abs() < 0.01);
|
||||
/// ```
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
pub struct Rosenbrock23<const D: usize> {
|
||||
/// Coefficient c₃₂ = 6 + √2 ≈ 7.414213562373095
|
||||
c32: f64,
|
||||
/// Coefficient d = 1/(2+√2) ≈ 0.29289321881345254 (this is gamma!)
|
||||
d: f64,
|
||||
/// Absolute tolerance for error estimation
|
||||
a_tol: f64,
|
||||
/// Relative tolerance for error estimation
|
||||
r_tol: f64,
|
||||
/// Strategy for updating the Jacobian
|
||||
jacobian_strategy: JacobianUpdateStrategy,
|
||||
/// Cached Jacobian from previous step
|
||||
cached_jacobian: Option<SMatrix<f64, D, D>>,
|
||||
/// Cached W matrix from previous step
|
||||
cached_w: Option<SMatrix<f64, D, D>>,
|
||||
/// Current step size (for detecting changes)
|
||||
cached_h: Option<f64>,
|
||||
/// Step counter for Jacobian update strategy
|
||||
steps_since_jacobian_update: usize,
|
||||
}
|
||||
|
||||
// TODO: Finish this
|
||||
fn step(&self, ode: &ODE<D>, _h: f64) -> (SVector<f64,D>, Option<f64>) {
|
||||
let next_y = ode.y.clone();
|
||||
let err = SVector::<f64, D>::zeros();
|
||||
(next_y, Some(err.norm()))
|
||||
impl<const D: usize> Rosenbrock23<D> {
|
||||
/// Create a new Rosenbrock23 integrator with default tolerances
|
||||
pub fn new() -> Self {
|
||||
Self {
|
||||
c32: 6.0 + 2.0_f64.sqrt(),
|
||||
d: 1.0 / (2.0 + 2.0_f64.sqrt()),
|
||||
a_tol: 1e-6,
|
||||
r_tol: 1e-3,
|
||||
jacobian_strategy: JacobianUpdateStrategy::default(),
|
||||
cached_jacobian: None,
|
||||
cached_w: None,
|
||||
cached_h: None,
|
||||
steps_since_jacobian_update: 0,
|
||||
}
|
||||
}
|
||||
|
||||
/// Set the absolute tolerance
|
||||
pub fn a_tol(mut self, a_tol: f64) -> Self {
|
||||
self.a_tol = a_tol;
|
||||
self
|
||||
}
|
||||
|
||||
/// Set the relative tolerance
|
||||
pub fn r_tol(mut self, r_tol: f64) -> Self {
|
||||
self.r_tol = r_tol;
|
||||
self
|
||||
}
|
||||
|
||||
/// Set the Jacobian update strategy
|
||||
pub fn jacobian_strategy(mut self, strategy: JacobianUpdateStrategy) -> Self {
|
||||
self.jacobian_strategy = strategy;
|
||||
self
|
||||
}
|
||||
|
||||
/// Decide if we should update the Jacobian on this step
|
||||
fn should_update_jacobian(&self, step_rejected: bool) -> bool {
|
||||
match self.jacobian_strategy {
|
||||
JacobianUpdateStrategy::EveryStep => true,
|
||||
JacobianUpdateStrategy::FirstAndRejection { update_interval } => {
|
||||
self.cached_jacobian.is_none()
|
||||
|| step_rejected
|
||||
|| self.steps_since_jacobian_update >= update_interval
|
||||
}
|
||||
JacobianUpdateStrategy::OnlyOnRejection => {
|
||||
self.cached_jacobian.is_none() || step_rejected
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<const D: usize> Default for Rosenbrock23<D> {
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
}
|
||||
|
||||
impl<const D: usize> Integrator<D> for Rosenbrock23<D> {
|
||||
const ORDER: usize = 2;
|
||||
const STAGES: usize = 2;
|
||||
const ADAPTIVE: bool = true;
|
||||
const DENSE: bool = true;
|
||||
|
||||
fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
|
||||
let t = ode.t;
|
||||
let uprev = ode.y;
|
||||
|
||||
// Compute Jacobian
|
||||
let j = compute_jacobian(&ode.f, t, uprev, &ode.params);
|
||||
|
||||
let dtgamma = h * self.d;
|
||||
|
||||
// Form W = I - dtγ*J
|
||||
let w = SMatrix::<T, D, D>::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
|
||||
let k1 = w_inv * fsalfirst;
|
||||
|
||||
// Stage 2: u = uprev + dt/2 * k₁
|
||||
let dto2 = h / 2.0;
|
||||
let u = uprev + k1 * dto2;
|
||||
|
||||
// Evaluate f₁ = f(u, t + dt/2)
|
||||
let f1 = (ode.f)(t + dto2, u, &ode.params);
|
||||
|
||||
// Stage 2: W * k₂ = f₁ - k₁ + J*k₁
|
||||
let rhs2 = f1 - k1;
|
||||
let k2_temp = w_inv * rhs2;
|
||||
let k2 = k2_temp + k1;
|
||||
|
||||
// Solution: u = uprev + dt * k₂
|
||||
let u_final = uprev + k2 * h;
|
||||
|
||||
// Error estimation
|
||||
// Evaluate fsallast = f(u_final, t + dt)
|
||||
let fsallast = (ode.f)(t + h, u_final, &ode.params);
|
||||
|
||||
// Ignoring dT (time derivative) for autonomous systems
|
||||
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₃)
|
||||
let dto6 = h / 6.0;
|
||||
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].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();
|
||||
|
||||
// Dense output: store k₁ and k₂
|
||||
let dense = Some(vec![k1, k2]);
|
||||
|
||||
StepResult {
|
||||
y: u_final,
|
||||
err: Some(error),
|
||||
dense,
|
||||
}
|
||||
}
|
||||
|
||||
fn interpolate<T: Scalar>(
|
||||
&self,
|
||||
t_start: T,
|
||||
t_end: T,
|
||||
dense: &[SVector<T, D>],
|
||||
t: T,
|
||||
) -> SVector<T, D> {
|
||||
// Second order interpolation using k₁ and k₂
|
||||
// For Rosenbrock methods, we use a simple Hermite interpolation
|
||||
let k1 = dense[0];
|
||||
let k2 = dense[1];
|
||||
|
||||
let h = t_end - t_start;
|
||||
let theta = (t - t_start) / h;
|
||||
|
||||
// Linear combination: y(t) ≈ y_n + θ*h*k₂ (first order)
|
||||
// For second order, we blend between k₁ and k₂:
|
||||
// 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
|
||||
k1 * (T::lift(1.0) - theta) + k2 * theta
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use approx::assert_relative_eq;
|
||||
use nalgebra::{Vector1, Vector2};
|
||||
|
||||
#[test]
|
||||
fn test_compute_jacobian_linear() {
|
||||
// Test on y' = -y (Jacobian should be -1)
|
||||
fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
|
||||
Vector1::new(-y[0])
|
||||
}
|
||||
|
||||
let j = compute_jacobian(&derivative, 0.0, Vector1::new(1.0), &());
|
||||
assert_relative_eq!(j[(0, 0)], -1.0, epsilon = 1e-6);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_compute_jacobian_nonlinear() {
|
||||
// Test on y' = y^2 at y=2 (Jacobian should be 2y = 4)
|
||||
fn derivative(_t: f64, y: Vector1<f64>, _p: &()) -> Vector1<f64> {
|
||||
Vector1::new(y[0] * y[0])
|
||||
}
|
||||
|
||||
let j = compute_jacobian(&derivative, 0.0, Vector1::new(2.0), &());
|
||||
assert_relative_eq!(j[(0, 0)], 4.0, epsilon = 1e-6);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_compute_jacobian_2d() {
|
||||
// Test on coupled system: y1' = y2, y2' = -y1
|
||||
// Jacobian should be [[0, 1], [-1, 0]]
|
||||
fn derivative(_t: f64, y: Vector2<f64>, _p: &()) -> Vector2<f64> {
|
||||
Vector2::new(y[1], -y[0])
|
||||
}
|
||||
|
||||
let j = compute_jacobian(&derivative, 0.0, Vector2::new(1.0, 0.0), &());
|
||||
assert_relative_eq!(j[(0, 0)], 0.0, epsilon = 1e-6);
|
||||
assert_relative_eq!(j[(0, 1)], 1.0, epsilon = 1e-6);
|
||||
assert_relative_eq!(j[(1, 0)], -1.0, epsilon = 1e-6);
|
||||
assert_relative_eq!(j[(1, 1)], 0.0, epsilon = 1e-6);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_rosenbrock23_simple_decay() {
|
||||
// Test y' = -y, y(0) = 1, h = 0.1
|
||||
// Analytical: y(0.1) = e^(-0.1) ≈ 0.904837418
|
||||
type Params = ();
|
||||
fn derivative(_t: f64, y: Vector1<f64>, _p: &Params) -> Vector1<f64> {
|
||||
Vector1::new(-y[0])
|
||||
}
|
||||
|
||||
let y0 = Vector1::new(1.0);
|
||||
let ode = ODE::new(&derivative, 0.0, 0.1, y0, ());
|
||||
let rb23 = Rosenbrock23::new();
|
||||
|
||||
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);
|
||||
println!("Error estimate: {:?}", err);
|
||||
|
||||
// Should be reasonably close (this is only one step with h=0.1)
|
||||
assert_relative_eq!(y_next[0], analytical, max_relative = 0.01);
|
||||
assert!(err.unwrap() < 1.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_rosenbrock23_convergence() {
|
||||
// Test order of convergence on y' = -y
|
||||
type Params = ();
|
||||
fn derivative(_t: f64, y: Vector1<f64>, _p: &Params) -> Vector1<f64> {
|
||||
Vector1::new(-y[0])
|
||||
}
|
||||
|
||||
let t_end = 1.0;
|
||||
let analytical = (-1.0_f64).exp();
|
||||
|
||||
let mut errors = Vec::new();
|
||||
let mut step_sizes = Vec::new();
|
||||
|
||||
// Test with decreasing step sizes
|
||||
for &n_steps in &[10, 20, 40, 80] {
|
||||
let h = t_end / n_steps as f64;
|
||||
let y0 = Vector1::new(1.0);
|
||||
let mut ode = ODE::new(&derivative, 0.0, t_end, y0, ());
|
||||
let rb23 = Rosenbrock23::new();
|
||||
|
||||
while ode.t < t_end - 1e-10 {
|
||||
let __step = rb23.step(&ode, h);
|
||||
let (y_next, _, _) = (__step.y, __step.err, __step.dense);
|
||||
ode.y = y_next;
|
||||
ode.t += h;
|
||||
}
|
||||
|
||||
let error = (ode.y[0] - analytical).abs();
|
||||
errors.push(error);
|
||||
step_sizes.push(h);
|
||||
}
|
||||
|
||||
// Check convergence rate
|
||||
// For a 2nd order method: error ∝ h^2
|
||||
// So log(error) = 2*log(h) + const
|
||||
// Slope should be approximately 2
|
||||
for i in 0..errors.len() - 1 {
|
||||
let rate =
|
||||
(errors[i].log10() - errors[i + 1].log10()) / (step_sizes[i].log10() - step_sizes[i + 1].log10());
|
||||
println!("Convergence rate between h={} and h={}: {}", step_sizes[i], step_sizes[i+1], rate);
|
||||
|
||||
// Should be close to 2 for a 2nd order method
|
||||
assert!(rate > 1.8 && rate < 2.2, "Convergence rate {} is not close to 2", rate);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_rosenbrock23_julia_linear_problem() {
|
||||
// Equivalent to Julia's prob_ode_linear: y' = 1.01*y, y(0) = 0.5, t ∈ [0, 1]
|
||||
// This matches the test in OrdinaryDiffEqRosenbrock/test/ode_rosenbrock_tests.jl
|
||||
type Params = ();
|
||||
fn derivative(_t: f64, y: Vector1<f64>, _p: &Params) -> Vector1<f64> {
|
||||
Vector1::new(1.01 * y[0])
|
||||
}
|
||||
|
||||
let y0 = Vector1::new(0.5);
|
||||
let t_end = 1.0;
|
||||
let analytical = |t: f64| 0.5 * (1.01 * t).exp();
|
||||
|
||||
// Test convergence with Julia's step sizes: (1/2)^(6:-1:3) = [1/64, 1/32, 1/16, 1/8]
|
||||
let step_sizes = vec![1.0/64.0, 1.0/32.0, 1.0/16.0, 1.0/8.0];
|
||||
let mut errors = Vec::new();
|
||||
|
||||
for &h in &step_sizes {
|
||||
let n_steps = (t_end / h) as usize;
|
||||
let actual_h = t_end / (n_steps as f64);
|
||||
|
||||
let mut ode = ODE::new(&derivative, 0.0, t_end, y0, ());
|
||||
let rb23 = Rosenbrock23::new();
|
||||
|
||||
for _ in 0..n_steps {
|
||||
let __step = rb23.step(&ode, actual_h);
|
||||
let (y_next, _, _) = (__step.y, __step.err, __step.dense);
|
||||
ode.y = y_next;
|
||||
ode.t += actual_h;
|
||||
}
|
||||
|
||||
let error = (ode.y[0] - analytical(t_end)).abs();
|
||||
println!("h = {:.6}, error = {:.3e}", h, error);
|
||||
errors.push(error);
|
||||
}
|
||||
|
||||
// Compute convergence order estimate like Julia's test_convergence does
|
||||
// Order = log(error[i+1]/error[i]) / log(h[i+1]/h[i])
|
||||
// Since h increases by factor of 2 each time and errors should decrease:
|
||||
// Order = log2(error[i+1]/error[i]) (negative since error decreases)
|
||||
// But we want positive order, so: Order = log2(error[i]/error[i+1])
|
||||
let mut orders = Vec::new();
|
||||
for i in 0..errors.len() - 1 {
|
||||
let order = (errors[i + 1] / errors[i]).log2(); // Larger h -> larger error
|
||||
orders.push(order);
|
||||
}
|
||||
|
||||
let avg_order = orders.iter().sum::<f64>() / orders.len() as f64;
|
||||
println!("Estimated order: {:.2}", avg_order);
|
||||
println!("Orders per step refinement: {:?}", orders);
|
||||
|
||||
// Julia tests: @test sim.𝒪est[:final]≈2 atol=0.2
|
||||
assert!((avg_order - 2.0).abs() < 0.2,
|
||||
"Convergence order {:.2} not within 0.2 of expected order 2", avg_order);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_rosenbrock23_adaptive_solve() {
|
||||
// This tests that the adaptive solver can efficiently solve prob_ode_linear
|
||||
use crate::controller::PIController;
|
||||
use crate::problem::Problem;
|
||||
|
||||
type Params = ();
|
||||
fn derivative(_t: f64, y: Vector1<f64>, _p: &Params) -> Vector1<f64> {
|
||||
Vector1::new(1.01 * y[0])
|
||||
}
|
||||
|
||||
let y0 = Vector1::new(0.5);
|
||||
let ode = crate::ode::ODE::new(&derivative, 0.0, 1.0, y0, ());
|
||||
|
||||
let rb23 = Rosenbrock23::new().a_tol(1e-3).r_tol(1e-3);
|
||||
let controller = PIController::default();
|
||||
|
||||
let mut problem = Problem::new(ode, rb23, controller);
|
||||
let solution = problem.solve();
|
||||
|
||||
println!("Adaptive solve completed in {} steps", solution.states.len());
|
||||
|
||||
assert!(solution.states.len() < 20,
|
||||
"Adaptive solve should complete in less than 20 steps, got {}",
|
||||
solution.states.len());
|
||||
|
||||
// Verify final value is accurate
|
||||
let analytical = 0.5 * (1.01_f64 * 1.0).exp();
|
||||
let final_val = solution.states[solution.states.len() - 1][0];
|
||||
assert_relative_eq!(final_val, analytical, max_relative = 1e-2);
|
||||
}
|
||||
}
|
||||
+117
-99
@@ -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<const D: usize> {
|
||||
r_tol: f64,
|
||||
}
|
||||
|
||||
impl<const D: usize> Default for Vern7<D>
|
||||
where
|
||||
Vern7<D>: Integrator<D>,
|
||||
{
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
}
|
||||
|
||||
impl<const D: usize> Vern7<D>
|
||||
where
|
||||
Vern7<D>: Integrator<D>,
|
||||
@@ -278,13 +288,9 @@ where
|
||||
const MAIN_STAGES: usize = 10;
|
||||
const EXTRA_STAGES: usize = 6;
|
||||
|
||||
fn step<P>(
|
||||
&self,
|
||||
ode: &ODE<D, P>,
|
||||
h: f64,
|
||||
) -> (SVector<f64, D>, Option<f64>, Option<Vec<SVector<f64, D>>>) {
|
||||
fn step<T: Scalar, P>(&self, ode: &ODE<T, D, P>, h: T) -> StepResult<T, D> {
|
||||
// Allocate storage for the 10 stages
|
||||
let mut k: Vec<SVector<f64, D>> = vec![SVector::<f64, D>::zeros(); Self::STAGES];
|
||||
let mut k: Vec<SVector<T, D>> = vec![SVector::<T, D>::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::<f64, D>::zeros();
|
||||
for i in 0..Self::STAGES {
|
||||
err += k[i] * Self::B_ERROR[i] * h;
|
||||
let mut err = SVector::<T, D>::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<T: Scalar>(
|
||||
&self,
|
||||
t_start: f64,
|
||||
t_end: f64,
|
||||
dense: &[SVector<f64, D>],
|
||||
t: f64,
|
||||
) -> SVector<f64, D> {
|
||||
t_start: T,
|
||||
t_end: T,
|
||||
dense: &[SVector<T, D>],
|
||||
t: T,
|
||||
) -> SVector<T, D> {
|
||||
// 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<T: Scalar>(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<P>(
|
||||
fn compute_extra_stages<T: Scalar, P>(
|
||||
&self,
|
||||
ode: &ODE<D, P>,
|
||||
t_start: f64,
|
||||
y_start: SVector<f64, D>,
|
||||
h: f64,
|
||||
main_stages: &[SVector<f64, D>],
|
||||
) -> Vec<SVector<f64, D>> {
|
||||
ode: &ODE<T, D, P>,
|
||||
t_start: T,
|
||||
y_start: SVector<T, D>,
|
||||
h: T,
|
||||
main_stages: &[SVector<T, D>],
|
||||
) -> Vec<SVector<T, D>> {
|
||||
// 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;
|
||||
}
|
||||
|
||||
+202
-4
@@ -1,19 +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)]
|
||||
@@ -43,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<f64>, _p: &Params| 5.0 - t,
|
||||
let value_too_high = FnCallback {
|
||||
event: &|t: f64, _y: &Vector2<f64>, _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
|
||||
@@ -149,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: Scalar>(_t: T, state: Vector6<T>, p: &(T,)) -> Vector6<T> {
|
||||
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<T: Scalar>() -> Vector6<T> {
|
||||
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<T: Scalar, const D: usize, P, S: Integrator<D> + Copy>(
|
||||
f: ProblemFunction<'_, T, D, P>,
|
||||
integrator: S,
|
||||
t0: T,
|
||||
t_end: T,
|
||||
y0: SVector<T, D>,
|
||||
params: P,
|
||||
) -> SVector<T, D> {
|
||||
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: Scalar>(t: T, _y: Vector1<T>, p: &(T,)) -> Vector1<T> {
|
||||
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: Scalar>(_t: T, y: Vector1<T>, _p: &()) -> Vector1<T> {
|
||||
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<f64> = orbit_y0();
|
||||
|
||||
// One propagation in `DualSVec64<6>` yields the whole state transition matrix.
|
||||
let (y_end, stm): (Vector6<f64>, Matrix6<f64>) = num_dual::jacobian(
|
||||
|y: Vector6<DualSVec64<6>>| {
|
||||
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<f64>| {
|
||||
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: Scalar>(_t: T, y: Vector1<T>, _p: &()) -> Vector1<T> {
|
||||
-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);
|
||||
}
|
||||
}
|
||||
+27
-17
@@ -1,27 +1,37 @@
|
||||
use nalgebra::SVector;
|
||||
|
||||
type ProblemFunction<'a, const D: usize, P> = &'a dyn Fn(f64, SVector<f64, D>, &P) -> SVector<f64, D>;
|
||||
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: Scalar>(t: T, y: SVector<T, D>, 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<T, D>, &P) -> SVector<T, D>;
|
||||
|
||||
/// 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<f64, D>,
|
||||
pub t: f64,
|
||||
pub struct ODE<'a, T: Scalar, const D: usize, P> {
|
||||
pub f: ProblemFunction<'a, T, D, P>,
|
||||
pub y: SVector<T, D>,
|
||||
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<f64, D>,
|
||||
f: ProblemFunction<'a, T, D, P>,
|
||||
t0: T,
|
||||
t_end: T,
|
||||
y0: SVector<T, D>,
|
||||
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<f64>, _p: &Params) -> Vector3<f64> {
|
||||
fn derivative<T: Scalar>(_t: T, y: Vector3<T>, _p: &Params) -> Vector3<T> {
|
||||
-y
|
||||
}
|
||||
|
||||
@@ -64,7 +74,7 @@ mod tests {
|
||||
type Params = (f64, bool);
|
||||
let params = (34.0, true);
|
||||
|
||||
fn derivative(t: f64, y: Vector3<f64>, p: &Params) -> Vector3<f64> {
|
||||
fn derivative<T: Scalar>(t: T, y: Vector3<T>, p: &Params) -> Vector3<T> {
|
||||
if p.1 {
|
||||
-y
|
||||
} else {
|
||||
|
||||
+86
-75
@@ -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<D>,
|
||||
C: Controller<T>,
|
||||
{
|
||||
ode: ODE<'a, D, P>,
|
||||
ode: ODE<'a, T, D, P>,
|
||||
integrator: S,
|
||||
controller: PIController,
|
||||
callbacks: Vec<Callback<'a, D, P>>,
|
||||
controller: C,
|
||||
callbacks: Vec<&'a dyn Callback<T, D, P>>,
|
||||
}
|
||||
|
||||
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<D> + Copy,
|
||||
C: Controller<T>,
|
||||
{
|
||||
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<f64> = vec![self.ode.t];
|
||||
let mut states: Vec<SVector<f64, D>> = vec![self.ode.y];
|
||||
let mut dense_coefficients: Vec<RefCell<Vec<SVector<f64, D>>>> = 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<T> = vec![self.ode.t];
|
||||
let mut states: Vec<SVector<T, D>> = vec![self.ode.y];
|
||||
let mut dense_coefficients: Vec<RefCell<Vec<SVector<T, D>>>> = 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 = <PIController as Controller<D>>::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 = <PIController as Controller<D>>::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<T, D, P>) -> 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<D>,
|
||||
{
|
||||
pub ode: &'a ODE<'a, D, P>,
|
||||
pub ode: &'a ODE<'a, T, D, P>,
|
||||
pub integrator: S,
|
||||
pub times: Vec<f64>,
|
||||
pub states: Vec<SVector<f64, D>>,
|
||||
pub dense: Vec<RefCell<Vec<SVector<f64, D>>>>,
|
||||
pub times: Vec<T>,
|
||||
pub states: Vec<SVector<T, D>>,
|
||||
pub dense: Vec<RefCell<Vec<SVector<T, D>>>>,
|
||||
}
|
||||
|
||||
impl<'a, S, const D: usize, P> Solution<'a, S, D, P>
|
||||
impl<T: Scalar, S, const D: usize, P> Solution<'_, T, S, D, P>
|
||||
where
|
||||
S: Integrator<D>,
|
||||
{
|
||||
pub fn interpolate(&self, t: f64) -> SVector<f64, D> {
|
||||
pub fn interpolate(&self, t: T) -> SVector<T, D> {
|
||||
// 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<f64> = 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<f64, 3>, _: &Params| 10.0 - y[0],
|
||||
let value_too_high = FnCallback {
|
||||
event: &|_: f64, y: &SVector<f64, 3>, _: &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!(
|
||||
|
||||
+119
@@ -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<f64, Self> {
|
||||
/// 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<T> Scalar for num_dual::$ty<T>
|
||||
where
|
||||
Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
|
||||
{
|
||||
#[inline(always)]
|
||||
fn re(&self) -> f64 {
|
||||
DualStruct::re(self)
|
||||
}
|
||||
|
||||
#[inline(always)]
|
||||
fn lift(x: f64) -> Self {
|
||||
Self::from(x)
|
||||
}
|
||||
}
|
||||
};
|
||||
($ty:ident, $($n:ident),+) => {
|
||||
impl<T: nalgebra::Scalar, $(const $n: usize),+> Scalar
|
||||
for num_dual::$ty<T, $(Const<$n>),+>
|
||||
where
|
||||
Self: RealField + Copy + DualNum<Primitive = f64> + DualStruct<Real = f64>,
|
||||
{
|
||||
#[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 = <Dual64 as Scalar>::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);
|
||||
}
|
||||
|
||||
}
|
||||
Reference in new issue
Block a user