Git hooks, format changes, etc.
This commit is contained in:
1 parent
328abfde1a
commit
cd3c40940e
15 files changed
+213
-110
No files matched your search
+10
-2
@@ -194,7 +194,11 @@ fn bench_interpolation(c: &mut Criterion) {
|
||||
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<_>>());
|
||||
black_box(
|
||||
(0..100)
|
||||
.map(|i| solution.interpolate(i as f64 * 0.05))
|
||||
.collect::<Vec<_>>(),
|
||||
);
|
||||
});
|
||||
});
|
||||
|
||||
@@ -206,7 +210,11 @@ fn bench_interpolation(c: &mut Criterion) {
|
||||
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<_>>());
|
||||
black_box(
|
||||
(0..100)
|
||||
.map(|i| solution.interpolate(i as f64 * 0.05))
|
||||
.collect::<Vec<_>>(),
|
||||
);
|
||||
});
|
||||
});
|
||||
|
||||
|
||||
+1
-1
@@ -1,7 +1,7 @@
|
||||
use criterion::{criterion_group, criterion_main, Criterion};
|
||||
|
||||
use ordinary_diffeq::prelude::*;
|
||||
use nalgebra::Vector6;
|
||||
use ordinary_diffeq::prelude::*;
|
||||
|
||||
fn bench_orbit(c: &mut Criterion) {
|
||||
let mu = 3.98600441500000e14;
|
||||
|
||||
@@ -1,7 +1,7 @@
|
||||
use criterion::{criterion_group, criterion_main, Criterion};
|
||||
|
||||
use ordinary_diffeq::prelude::*;
|
||||
use nalgebra::Vector1;
|
||||
use ordinary_diffeq::prelude::*;
|
||||
|
||||
fn bench_simple_1d(c: &mut Criterion) {
|
||||
type Params = (f64,);
|
||||
|
||||
@@ -172,7 +172,11 @@ fn bench_vern7_interpolation(c: &mut Criterion) {
|
||||
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<_>>());
|
||||
black_box(
|
||||
(0..100)
|
||||
.map(|i| solution.interpolate(i as f64 * 0.05))
|
||||
.collect::<Vec<_>>(),
|
||||
);
|
||||
});
|
||||
});
|
||||
|
||||
@@ -183,7 +187,11 @@ fn bench_vern7_interpolation(c: &mut Criterion) {
|
||||
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<_>>());
|
||||
black_box(
|
||||
(0..100)
|
||||
.map(|i| solution.interpolate(i as f64 * 0.05))
|
||||
.collect::<Vec<_>>(),
|
||||
);
|
||||
});
|
||||
});
|
||||
|
||||
|
||||
Executable
+11
@@ -0,0 +1,11 @@
|
||||
#!/bin/sh
|
||||
# Fast checks only. On a warm tree these are about a second between them; a hook
|
||||
# that takes a minute is a hook that gets --no-verify. The doctests are in
|
||||
# pre-push instead, since each one links its own binary and they cost ~15s.
|
||||
#
|
||||
# Install: git config core.hooksPath hooks
|
||||
set -e
|
||||
|
||||
cargo fmt --check
|
||||
cargo clippy --all-targets -- -D warnings
|
||||
cargo test --quiet --lib
|
||||
Executable
+9
@@ -0,0 +1,9 @@
|
||||
#!/bin/sh
|
||||
# Everything the pre-commit hook is too slow to run. The doctests matter here
|
||||
# because the readme is the crate's front page, so its examples are compiled.
|
||||
#
|
||||
# Install: git config core.hooksPath hooks
|
||||
set -e
|
||||
|
||||
cargo test --quiet
|
||||
cargo doc --quiet --no-deps
|
||||
@@ -127,6 +127,27 @@ This is a solid foundation! The roadmap builds on this.
|
||||
|
||||
## Development Tips
|
||||
|
||||
### Checks, and the Hooks That Run Them
|
||||
|
||||
The tree is rustfmt-clean, clippy-clean under `-D warnings`, and builds its docs
|
||||
without warnings. All three stay that way, and a pre-existing warning is not an
|
||||
exemption — fix it or mark it deliberately, as the Butcher tableaux are marked
|
||||
`#[rustfmt::skip]` so rustfmt leaves their hand alignment alone.
|
||||
|
||||
The hooks in `hooks/` run them. Git does not track where it looks for hooks, so
|
||||
each clone needs one command:
|
||||
|
||||
```sh
|
||||
git config core.hooksPath hooks
|
||||
```
|
||||
|
||||
`pre-commit` is `cargo fmt --check`, `cargo clippy --all-targets -- -D warnings`
|
||||
and `cargo test --lib`, about a second between them on a warm tree. `pre-push`
|
||||
adds the full `cargo test` and `cargo doc --no-deps`, which cost around 25
|
||||
seconds — the doctests are slow because each one links its own binary, and they
|
||||
matter because the readme is the crate's front page and its examples are
|
||||
compiled.
|
||||
|
||||
### Testing Strategy
|
||||
|
||||
Every feature should have:
|
||||
|
||||
+44
-48
@@ -76,9 +76,7 @@ impl<T: Scalar> Controller<T> for PIController<T> {
|
||||
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))
|
||||
}
|
||||
}
|
||||
|
||||
@@ -124,7 +122,6 @@ impl<T: Scalar> Default for PIController<T> {
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// PID (Proportional-Integral-Derivative) step size controller
|
||||
///
|
||||
/// The PID controller is an advanced adaptive time-stepping controller that provides
|
||||
@@ -171,20 +168,20 @@ impl<T: Scalar> Default for PIController<T> {
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
pub struct PIDController<T: Scalar = f64> {
|
||||
// PID Coefficients
|
||||
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
|
||||
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)
|
||||
pub factor_c2: f64, // 1/max_factor (maximum step increase)
|
||||
pub h_max: f64, // Maximum allowed step size
|
||||
pub safety_factor: f64, // Safety factor (typically 0.9)
|
||||
pub factor_c1: f64, // 1/min_factor (maximum step decrease)
|
||||
pub factor_c2: f64, // 1/max_factor (maximum step increase)
|
||||
pub h_max: f64, // Maximum allowed step size
|
||||
pub safety_factor: f64, // Safety factor (typically 0.9)
|
||||
|
||||
// Error history for PID control
|
||||
pub err_old: f64, // ε_{n-1}: previous step error
|
||||
pub err_older: f64, // ε_{n-2}: error two steps ago
|
||||
pub h_old: f64, // h_{n-1}: previous step size
|
||||
pub err_old: f64, // ε_{n-1}: previous step error
|
||||
pub err_older: f64, // ε_{n-2}: error two steps ago
|
||||
pub h_old: f64, // h_{n-1}: previous step size
|
||||
|
||||
// Next step guess
|
||||
pub next_step_guess: TryStep<T>,
|
||||
@@ -199,15 +196,15 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
|
||||
let factor = if self.err_old <= 0.0 {
|
||||
// First step: use only proportional control
|
||||
let factor_11 = err.powf(self.p_term);
|
||||
self.factor_c2.max(
|
||||
self.factor_c1.min(factor_11 / self.safety_factor)
|
||||
)
|
||||
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.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)
|
||||
self.factor_c1
|
||||
.min(factor_11 * factor_12 / self.safety_factor),
|
||||
)
|
||||
} else {
|
||||
// Full PID control (proportional + integral + derivative)
|
||||
@@ -220,7 +217,8 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
|
||||
1.0
|
||||
};
|
||||
self.factor_c2.max(
|
||||
self.factor_c1.min(factor_11 * factor_12 * factor_13 / self.safety_factor)
|
||||
self.factor_c1
|
||||
.min(factor_11 * factor_12 * factor_13 / self.safety_factor),
|
||||
)
|
||||
};
|
||||
|
||||
@@ -230,7 +228,7 @@ impl<T: Scalar> Controller<T> for PIDController<T> {
|
||||
|
||||
// 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.err_old = err.max(1.0e-4); // Prevent very small values
|
||||
self.h_old = h.re();
|
||||
|
||||
// Apply maximum step size limit
|
||||
@@ -297,9 +295,9 @@ impl<T: Scalar> PIDController<T> {
|
||||
factor_c2: 1.0 / max_factor,
|
||||
h_max: h_max.abs(),
|
||||
safety_factor,
|
||||
err_old: 0.0, // No history initially
|
||||
err_older: 0.0, // No history initially
|
||||
h_old: 0.0, // No history initially
|
||||
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(T::lift(initial_h)),
|
||||
}
|
||||
}
|
||||
@@ -317,14 +315,14 @@ impl<T: Scalar> PIDController<T> {
|
||||
pub fn for_order(order: usize) -> Self {
|
||||
let k_plus_1 = (order + 1) as f64;
|
||||
Self::new(
|
||||
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
|
||||
0.9, // safety_factor
|
||||
1e-4, // initial_h
|
||||
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
|
||||
0.9, // safety_factor
|
||||
1e-4, // initial_h
|
||||
)
|
||||
}
|
||||
}
|
||||
@@ -342,19 +340,18 @@ impl<T: Scalar> Default for PIDController<T> {
|
||||
/// For method-specific tuning, use `PIDController::for_order(order)` instead.
|
||||
fn default() -> Self {
|
||||
Self::new(
|
||||
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
|
||||
0.9, // safety_factor
|
||||
1e-4, // initial_h
|
||||
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
|
||||
0.9, // safety_factor
|
||||
1e-4, // initial_h
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
@@ -406,7 +403,7 @@ mod tests {
|
||||
// Default uses conservative PID coefficients
|
||||
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
|
||||
assert_eq!(controller.d_term, 0.01); // True PID with derivative term
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -443,7 +440,8 @@ mod tests {
|
||||
|
||||
for &err in &errors {
|
||||
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);
|
||||
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());
|
||||
@@ -555,7 +553,8 @@ mod tests {
|
||||
|
||||
for &err in &errors {
|
||||
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);
|
||||
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();
|
||||
@@ -574,11 +573,8 @@ mod tests {
|
||||
|
||||
// Compute standard deviation of log step size ratios
|
||||
let pi_mean: f64 = pi_ratios.iter().sum::<f64>() / pi_ratios.len() as f64;
|
||||
let pi_variance: f64 = pi_ratios
|
||||
.iter()
|
||||
.map(|r| (r - pi_mean).powi(2))
|
||||
.sum::<f64>()
|
||||
/ pi_ratios.len() as f64;
|
||||
let pi_variance: f64 =
|
||||
pi_ratios.iter().map(|r| (r - pi_mean).powi(2)).sum::<f64>() / pi_ratios.len() as f64;
|
||||
let pi_std = pi_variance.sqrt();
|
||||
|
||||
let pid_mean: f64 = pid_ratios.iter().sum::<f64>() / pid_ratios.len() as f64;
|
||||
|
||||
+20
-8
@@ -64,7 +64,7 @@ pub struct BS3<const D: usize> {
|
||||
impl<const D: usize> Default for BS3<D>
|
||||
where
|
||||
BS3<D>: Integrator<D>,
|
||||
{
|
||||
{
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
@@ -120,6 +120,10 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
|
||||
// Row 2: [1/2]
|
||||
// Row 3: [0, 3/4]
|
||||
// Row 4: [2/9, 1/3, 4/9]
|
||||
// The tableau is aligned by hand, so rustfmt leaves it alone.
|
||||
#[rustfmt::skip]
|
||||
// The tableau is aligned by hand, so rustfmt leaves it alone.
|
||||
#[rustfmt::skip]
|
||||
const A: &'a [f64] = &[
|
||||
1.0 / 2.0, // a[1,0]
|
||||
0.0, // a[2,0]
|
||||
@@ -130,6 +134,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
|
||||
];
|
||||
|
||||
// Solution weights (3rd order)
|
||||
#[rustfmt::skip]
|
||||
#[rustfmt::skip]
|
||||
const B: &'a [f64] = &[
|
||||
2.0 / 9.0, // b[0]
|
||||
1.0 / 3.0, // b[1]
|
||||
@@ -138,6 +144,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
|
||||
];
|
||||
|
||||
// Error estimate weights (difference between 3rd and 2nd order)
|
||||
#[rustfmt::skip]
|
||||
#[rustfmt::skip]
|
||||
const B_ERROR: &'a [f64] = &[
|
||||
2.0 / 9.0 - 7.0 / 24.0, // b[0] - b*[0]
|
||||
1.0 / 3.0 - 1.0 / 4.0, // b[1] - b*[1]
|
||||
@@ -146,6 +154,8 @@ impl<'a, const D: usize> BS3Integrator<'a> for BS3<D> {
|
||||
];
|
||||
|
||||
// Stage times
|
||||
#[rustfmt::skip]
|
||||
#[rustfmt::skip]
|
||||
const C: &'a [f64] = &[
|
||||
0.0, // c[0]
|
||||
1.0 / 2.0, // c[1]
|
||||
@@ -191,7 +201,9 @@ where
|
||||
|
||||
// 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
|
||||
+ (k[0] * T::lift(Self::A[3]) + k[1] * T::lift(Self::A[4]) + k[2] * T::lift(Self::A[5]))
|
||||
+ (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);
|
||||
|
||||
@@ -225,10 +237,10 @@ where
|
||||
// BS3 uses standard cubic Hermite interpolation with derivatives at endpoints
|
||||
// Store: y0, y1, f0=k[0], f1=k[3] (FSAL)
|
||||
let dense_coeffs = vec![
|
||||
ode.y, // y0 at start of step
|
||||
next_y, // y1 at end of step
|
||||
k[0], // f(t0, y0) - derivative at start
|
||||
k[3], // f(t1, y1) - derivative at end (FSAL)
|
||||
ode.y, // y0 at start of step
|
||||
next_y, // y1 at end of step
|
||||
k[0], // f(t0, y0) - derivative at start
|
||||
k[3], // f(t1, y1) - derivative at end (FSAL)
|
||||
];
|
||||
|
||||
StepResult {
|
||||
@@ -304,7 +316,7 @@ mod tests {
|
||||
let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
|
||||
|
||||
let bs3 = BS3::new();
|
||||
let h = 0.001; // Smaller step size for tighter tolerances
|
||||
let h = 0.001; // Smaller step size for tighter tolerances
|
||||
|
||||
let __step = bs3.step(&ode, h);
|
||||
let (y_next, err, dense) = (__step.y, __step.err, __step.dense);
|
||||
@@ -334,7 +346,7 @@ mod tests {
|
||||
let ode = ODE::new(&derivative, 0.0, 1.0, y0, ());
|
||||
|
||||
let bs3 = BS3::new();
|
||||
let h = 0.001; // Smaller step size
|
||||
let h = 0.001; // Smaller step size
|
||||
|
||||
let __step = bs3.step(&ode, h);
|
||||
let (_y_next, _err, dense) = (__step.y, __step.err, __step.dense);
|
||||
|
||||
@@ -212,7 +212,11 @@ where
|
||||
for (j, item) in k.iter().enumerate().take(i) {
|
||||
y_term += item * T::lift(Self::A[(i * (i - 1)) / 2 + j]);
|
||||
}
|
||||
k[i] = (ode.f)(ode.t + T::lift(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 * T::lift(Self::B[i]);
|
||||
|
||||
@@ -28,10 +28,11 @@ impl Default for JacobianUpdateStrategy {
|
||||
|
||||
/// 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
|
||||
/// 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)
|
||||
/// 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,
|
||||
@@ -443,12 +444,21 @@ mod tests {
|
||||
// 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);
|
||||
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);
|
||||
assert!(
|
||||
rate > 1.8 && rate < 2.2,
|
||||
"Convergence rate {} is not close to 2",
|
||||
rate
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -466,7 +476,7 @@ mod tests {
|
||||
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 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 {
|
||||
@@ -504,8 +514,11 @@ mod tests {
|
||||
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);
|
||||
assert!(
|
||||
(avg_order - 2.0).abs() < 0.2,
|
||||
"Convergence order {:.2} not within 0.2 of expected order 2",
|
||||
avg_order
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -528,11 +541,16 @@ mod tests {
|
||||
let mut problem = Problem::new(ode, rb23, controller);
|
||||
let solution = problem.solve();
|
||||
|
||||
println!("Adaptive solve completed in {} steps", solution.states.len());
|
||||
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());
|
||||
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();
|
||||
|
||||
+16
-14
@@ -80,7 +80,7 @@ pub struct Vern7<const D: usize> {
|
||||
impl<const D: usize> Default for Vern7<D>
|
||||
where
|
||||
Vern7<D>: Integrator<D>,
|
||||
{
|
||||
{
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
@@ -141,6 +141,8 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
|
||||
// Stage 8: [a81, 0, a83, a84, a85, a86, a87]
|
||||
// Stage 9: [a91, 0, a93, a94, a95, a96, a97, a98]
|
||||
// Stage 10: [a101, 0, a103, a104, a105, a106, a107, 0, 0]
|
||||
// The tableau is aligned by hand, so rustfmt leaves it alone.
|
||||
#[rustfmt::skip]
|
||||
const A: &'a [f64] = &[
|
||||
// Stage 2
|
||||
0.005,
|
||||
@@ -177,6 +179,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
|
||||
];
|
||||
|
||||
// Error estimate weights (difference between 7th and 6th order: b - b*)
|
||||
#[rustfmt::skip]
|
||||
const B_ERROR: &'a [f64] = &[
|
||||
0.002547011879931045, // b1 - b*1
|
||||
0.0, // b2 - b*2
|
||||
@@ -191,6 +194,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
|
||||
];
|
||||
|
||||
// Time nodes (c coefficients)
|
||||
#[rustfmt::skip]
|
||||
const C: &'a [f64] = &[
|
||||
0.0, // c1
|
||||
0.005, // c2
|
||||
@@ -210,6 +214,7 @@ impl<'a, const D: usize> Vern7Integrator<'a> for Vern7<D> {
|
||||
|
||||
impl<'a, const D: usize> Vern7ExtraStages<'a> for Vern7<D> {
|
||||
// Time nodes for extra stages
|
||||
#[rustfmt::skip]
|
||||
const C_EXTRA: &'a [f64] = &[
|
||||
1.0, // c11
|
||||
0.29, // c12
|
||||
@@ -375,6 +380,8 @@ where
|
||||
}
|
||||
}
|
||||
|
||||
// One interpolation coefficient per line, likewise.
|
||||
#[rustfmt::skip]
|
||||
fn interpolate<T: Scalar>(
|
||||
&self,
|
||||
t_start: T,
|
||||
@@ -744,7 +751,7 @@ mod tests {
|
||||
}
|
||||
|
||||
let y0 = Vector1::new(1.0);
|
||||
let t_end: f64 = 1.0; // Longer interval to get larger errors
|
||||
let t_end: f64 = 1.0; // Longer interval to get larger errors
|
||||
let exact = t_end.exp();
|
||||
|
||||
let step_sizes: [f64; 3] = [0.2, 0.1, 0.05];
|
||||
@@ -771,16 +778,8 @@ mod tests {
|
||||
let ratio2 = errors[1] / errors[2];
|
||||
|
||||
// Allow some tolerance (expect ratio between 64 and 256)
|
||||
assert!(
|
||||
ratio1 > 64.0 && ratio1 < 256.0,
|
||||
"First ratio: {}",
|
||||
ratio1
|
||||
);
|
||||
assert!(
|
||||
ratio2 > 64.0 && ratio2 < 256.0,
|
||||
"Second ratio: {}",
|
||||
ratio2
|
||||
);
|
||||
assert!(ratio1 > 64.0 && ratio1 < 256.0, "First ratio: {}", ratio1);
|
||||
assert!(ratio2 > 64.0 && ratio2 < 256.0, "Second ratio: {}", ratio2);
|
||||
}
|
||||
|
||||
#[test]
|
||||
@@ -828,7 +827,7 @@ mod tests {
|
||||
Vector2::new(y[1], -y[0])
|
||||
}
|
||||
|
||||
let y0 = Vector2::new(1.0, 0.0); // Start at maximum displacement, zero velocity
|
||||
let y0 = Vector2::new(1.0, 0.0); // Start at maximum displacement, zero velocity
|
||||
|
||||
// Period of harmonic oscillator is 2π
|
||||
let period = 2.0 * std::f64::consts::PI;
|
||||
@@ -854,7 +853,10 @@ mod tests {
|
||||
|
||||
println!("Initial energy: {}", energy_initial);
|
||||
println!("Final energy: {}", energy_final);
|
||||
println!("Energy drift after {} periods: {:.2e}", num_periods, energy_drift);
|
||||
println!(
|
||||
"Energy drift after {} periods: {:.2e}",
|
||||
num_periods, energy_drift
|
||||
);
|
||||
println!("Number of steps: {}", solution.times.len());
|
||||
|
||||
// Energy should be conserved to high precision (< 1e-7 relative error over 1000 periods)
|
||||
|
||||
+8
-2
@@ -19,8 +19,8 @@ pub mod prelude {
|
||||
pub use super::integrator::method::{Method, MethodName};
|
||||
pub use super::integrator::rosenbrock::Rosenbrock23;
|
||||
pub use super::integrator::vern7::Vern7;
|
||||
pub use super::ode::ODE;
|
||||
pub use super::integrator::{Integrator, StepResult, Tolerances};
|
||||
pub use super::ode::ODE;
|
||||
pub use super::problem::{Problem, Solution};
|
||||
pub use super::scalar::Scalar;
|
||||
}
|
||||
@@ -376,7 +376,13 @@ mod autodiff_tests {
|
||||
}
|
||||
|
||||
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 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();
|
||||
|
||||
+25
-16
@@ -91,8 +91,7 @@ where
|
||||
// crossing starts from rather than at — see `EVENT_PROBE_FRACTION`.
|
||||
let mut last_fired: Vec<Option<f64>> = vec![None; self.callbacks.len()];
|
||||
while !self.ode.finished && 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 (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
|
||||
@@ -149,7 +148,9 @@ where
|
||||
// The event condition is a real-valued question, so it is
|
||||
// asked in the reals — see this module's header.
|
||||
let before = if search_start.re() == 0.0 {
|
||||
callback.event(self.ode.t, &self.ode.y, &self.ode.params).re()
|
||||
callback
|
||||
.event(self.ode.t, &self.ode.y, &self.ode.params)
|
||||
.re()
|
||||
} else {
|
||||
let probe = self.integrator.step(&self.ode, search_start).y;
|
||||
callback
|
||||
@@ -171,7 +172,9 @@ where
|
||||
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()
|
||||
callback
|
||||
.event(self.ode.t + h, &test_y, &self.ode.params)
|
||||
.re()
|
||||
};
|
||||
Some(
|
||||
find_root_brent(search_start.re(), curr_step.re(), &f, &mut convergency)
|
||||
@@ -366,6 +369,7 @@ mod tests {
|
||||
use crate::callback::{stop, FnCallback};
|
||||
use crate::controller::PIController;
|
||||
use crate::integrator::dormand_prince::DormandPrince45;
|
||||
use crate::integrator::vern7::Vern7;
|
||||
use crate::integrator::Tolerances;
|
||||
use approx::assert_relative_eq;
|
||||
use nalgebra::{Matrix3, Vector1, Vector2, Vector3};
|
||||
@@ -440,8 +444,7 @@ mod tests {
|
||||
effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
|
||||
};
|
||||
|
||||
let mut problem =
|
||||
Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
|
||||
let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
|
||||
let solution = problem.solve();
|
||||
|
||||
// Five units of integration before the jump, five after.
|
||||
@@ -464,8 +467,7 @@ mod tests {
|
||||
effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = Vector1::new(100.0),
|
||||
};
|
||||
|
||||
let mut problem =
|
||||
Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
|
||||
let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut jump);
|
||||
let solution = problem.solve();
|
||||
|
||||
let at_crossing: Vec<usize> = (0..solution.times.len())
|
||||
@@ -650,8 +652,7 @@ mod tests {
|
||||
effect: &|ode: &mut ODE<'_, f64, 1, ()>| ode.y = -ode.y,
|
||||
};
|
||||
|
||||
let mut problem =
|
||||
Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip);
|
||||
let mut problem = Problem::new(ode, dp45, PIController::default()).with_callback(&mut flip);
|
||||
let solution = problem.solve();
|
||||
|
||||
// y rises to 1 at t = 1 and flips to -1, which the effect leaves on the
|
||||
@@ -665,7 +666,11 @@ mod tests {
|
||||
assert_relative_eq!(jumps[0], 1.0, max_relative = 1e-9);
|
||||
assert_relative_eq!(jumps[1], 3.0, max_relative = 1e-9);
|
||||
assert_relative_eq!(solution.times.last().unwrap(), &3.0, max_relative = 1e-12);
|
||||
assert_relative_eq!(solution.states.last().unwrap()[0], -1.0, max_relative = 1e-8);
|
||||
assert_relative_eq!(
|
||||
solution.states.last().unwrap()[0],
|
||||
-1.0,
|
||||
max_relative = 1e-8
|
||||
);
|
||||
}
|
||||
|
||||
/// The controller's error history is discarded when an effect changes the
|
||||
@@ -708,8 +713,8 @@ mod tests {
|
||||
effect,
|
||||
};
|
||||
let ode = ODE::new(&unit_rate, 0.0, 10.0, Vector1::new(0.0), ());
|
||||
let mut problem = Problem::new(ode, DormandPrince45::new(), spy)
|
||||
.with_callback(&mut callback);
|
||||
let mut problem =
|
||||
Problem::new(ode, DormandPrince45::new(), spy).with_callback(&mut callback);
|
||||
problem.solve();
|
||||
resets.get()
|
||||
}
|
||||
@@ -739,7 +744,7 @@ mod tests {
|
||||
}
|
||||
|
||||
let ode = ODE::new(&free_fall, 0.0, 2.0, Vector2::new(1.0, 0.0), G);
|
||||
let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12);
|
||||
let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
|
||||
let mut bounce = FnCallback {
|
||||
event: &|_t: f64, y: &Vector2<f64>, _p: &f64| y[0],
|
||||
effect: &|ode: &mut ODE<'_, f64, 2, f64>| {
|
||||
@@ -803,7 +808,7 @@ mod tests {
|
||||
let omega = Vector3::new(0.0, 0.0, 1.0);
|
||||
let t_end = 5.0; // past half a turn, short of a full one
|
||||
let ode = ODE::new(&mrp_rate, 0.0, t_end, Vector3::zeros(), omega);
|
||||
let vern7 = crate::integrator::vern7::Vern7::new().a_tol(1e-12).r_tol(1e-12);
|
||||
let vern7 = Vern7::new().a_tol(1e-12).r_tol(1e-12);
|
||||
|
||||
let mut shadow = FnCallback {
|
||||
event: &|_t: f64, s: &Vector3<f64>, _p: &Vector3<f64>| 1.0 - s.norm(),
|
||||
@@ -830,7 +835,11 @@ mod tests {
|
||||
// And the attitude at the end is the one the body is actually in.
|
||||
let expected = -Vector3::new(0.0, 0.0, 1.0 / (t_end / 4.0).tan());
|
||||
for i in 0..3 {
|
||||
assert_relative_eq!(solution.states.last().unwrap()[i], expected[i], epsilon = 1e-9);
|
||||
assert_relative_eq!(
|
||||
solution.states.last().unwrap()[i],
|
||||
expected[i],
|
||||
epsilon = 1e-9
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@@ -115,5 +115,4 @@ mod tests {
|
||||
// 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