Added energy conservation test

roadmap/README.md

@@ -50,7 +50,7 @@ Each feature below links to a detailed implementation plan in the `features/` di
 
 ### Controllers
 
-- [x] **[PID Controller](features/04-pid-controller.md)** ✅ COMPLETED
+- [ ] **[PID Controller](features/04-pid-controller.md)**
   - Proportional-Integral-Derivative step size controller
   - Better stability than PI controller for difficult problems
   - **Dependencies**: None
@@ -329,15 +329,14 @@ Each algorithm implementation should include:
 ## Progress Tracking
 
 Total Features: 38
-- Tier 1: 8 features (3/8 complete) ✅
+- Tier 1: 8 features (2/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: 5.3% (2/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)
 
 Last updated: 2025-10-24

roadmap/features/04-pid-controller.md

@@ -1,16 +1,5 @@
 # Feature: PID Controller
 
-**Status**: ✅ COMPLETED (2025-10-24)
-
-**Implementation Summary**:
-- ✅ PIDController struct with beta1, beta2, beta3 coefficients and error history
-- ✅ Full Controller trait implementation with progressive bootstrap (P → PI → PID)
-- ✅ Constructor methods: new(), default(), for_order()
-- ✅ Reset method for clearing error history
-- ✅ Comprehensive test suite with 9 tests including PI vs PID comparisons
-- ✅ Exported in prelude
-- ✅ Complete documentation with mathematical formulation and usage guidance
-
 ## Overview
 
 The PID (Proportional-Integral-Derivative) step size controller is an advanced adaptive time-stepping controller that provides better stability and efficiency than the basic PI controller, especially for difficult or oscillatory problems.
@@ -90,97 +79,93 @@ pub struct PIDController {
 
 ### Core Controller
 
-- [x] Define `PIDController` struct ✅
+- [ ] Define `PIDController` struct
-  - [x] Add beta1, beta2, beta3 coefficients ✅
+  - [ ] Add beta1, beta2, beta3 coefficients
-  - [x] Add constraint fields (factor_c1, factor_c2, h_max, safety_factor) ✅
+  - [ ] Add constraint fields (factor_min, factor_max, h_max, safety)
-  - [x] Add state fields (err_old, err_older, h_old) ✅
+  - [ ] Add state fields (err_old, err_older, h_old)
-  - [x] Add next_step_guess field ✅
+  - [ ] Add next_step_guess field
 
-- [x] Implement `Controller<D>` trait ✅
+- [ ] Implement `Controller<D>` trait
-  - [x] `determine_step()` method ✅
+  - [ ] `determine_step()` method
-    - [x] Handle first step (no history) - proportional only ✅
+    - [ ] Handle first step (no history)
-    - [x] Handle second step (partial history) - PI control ✅
+    - [ ] Handle second step (partial history)
-    - [x] Full PID formula for subsequent steps ✅
+    - [ ] Full PID formula for subsequent steps
-    - [x] Apply safety factor and limits ✅
+    - [ ] Apply safety factor and limits
-    - [x] Update error history on acceptance only ✅
+    - [ ] Update error history
-    - [x] Return TryStep::Accepted or NotYetAccepted ✅
+    - [ ] Return TryStep::Accepted or NotYetAccepted
 
-- [x] Constructor methods ✅
+- [ ] Constructor methods
-  - [x] `new()` with all parameters ✅
+  - [ ] `new()` with all parameters
-  - [x] `default()` with H312 coefficients (PI controller) ✅
+  - [ ] `default()` with standard coefficients
-  - [x] `for_order()` - Gustafsson coefficients scaled by method order ✅
+  - [ ] `for_order()` - scale coefficients by method order
 
-- [x] Helper methods ✅
+- [ ] Helper methods
-  - [x] `reset()` - clear history (for algorithm switching) ✅
+  - [ ] `reset()` - clear history (for algorithm switching)
-  - [x] State correctly updated after accepted/rejected steps ✅
+  - [ ] Update state after accepted/rejected steps
 
 ### Standard Coefficient Sets
 
 Different coefficient sets for different problem classes:
 
-- [x] **Default (Conservative PID)** ✅:
+- [ ] **Default (H312)**:
-  - β₁ = 0.07, β₂ = 0.04, β₃ = 0.01
+  - β₁ = 1/4, β₂ = 1/4, β₃ = 0
-  - True PID with conservative coefficients
+  - Actually a PI controller with specific tuning
-  - Good general-purpose choice for orders 5-7
+  - Good general-purpose choice
-  - Implemented in `default()`
 
-- [ ] **H211** (Future):
+- [ ] **H211**:
   - β₁ = 1/6, β₂ = 1/6, β₃ = 0
   - More conservative
-  - Can be created with `new()`
 
-- [x] **Full PID (Gustafsson)** ✅:
+- [ ] **Full PID (Gustafsson)**:
   - β₁ = 0.49/(k+1)
   - β₂ = 0.34/(k+1)
   - β₃ = 0.10/(k+1)
   - True PID behavior
-  - Implemented in `for_order()`
 
 ### Integration
 
-- [x] Export PIDController in prelude ✅
+- [ ] Export PIDController in prelude
-- [x] Problem already accepts any Controller trait ✅
+- [ ] Update Problem to accept any Controller trait
-- [ ] Examples using PID controller (Future enhancement)
+- [ ] Examples using PID controller
 
 ### Testing
 
-- [x] **Comparison test: Smooth problem** ✅
+- [ ] **Comparison test: Smooth problem**
-  - [x] Run exponential decay with PI and PID ✅
+  - [ ] Run exponential decay with PI and PID
-  - [x] Both perform similarly ✅
+  - [ ] Both should perform similarly
-  - [x] Verified PID doesn't hurt performance ✅
+  - [ ] Verify PID doesn't hurt performance
 
-- [x] **Oscillatory problem test** ✅
+- [ ] **Oscillatory problem test**
-  - [x] Oscillatory error pattern test ✅
+  - [ ] Problem that causes PI to oscillate step sizes
-  - [x] PID has similar or better step size stability ✅
+  - [ ] Example: y'' + ω²y = 0 with varying ω
-  - [x] Standard deviation comparison test ✅
+  - [ ] PID should have smoother step size evolution
-  - [ ] Full ODE integration test (Future enhancement)
+  - [ ] Plot step size vs time for both
 
-- [x] **Step rejection handling** ✅
+- [ ] **Step rejection handling**
-  - [x] Verified history NOT updated after rejection ✅
+  - [ ] Verify history updated correctly after rejection
-  - [x] Test passes for rejection scenario ✅
+  - [ ] Doesn't blow up or get stuck
 
-- [x] **Reset test** ✅
+- [ ] **Reset test**
-  - [x] Verified reset() clears history appropriately ✅
+  - [ ] Algorithm switching scenario
-  - [x] Test passes ✅
+  - [ ] Verify reset() clears history appropriately
 
-- [x] **Bootstrap test** ✅
+- [ ] **Coefficient tuning test**
-  - [x] Verified P → PI → PID progression ✅
+  - [ ] Try different β values
-  - [x] Error history builds correctly ✅
+  - [ ] Verify stability bounds
+  - [ ] Document which work best for which problems
 
 ### Benchmarking
 
-- [ ] Add PID option to existing benchmarks (Future enhancement)
+- [ ] Add PID option to existing benchmarks
-- [ ] Compare step count and function evaluations vs PI (Future enhancement)
+- [ ] Compare step count and function evaluations vs PI
-- [ ] Measure overhead (should be negligible) (Future enhancement)
+- [ ] Measure overhead (should be negligible)
 
 ### Documentation
 
-- [x] Docstring explaining PID control ✅
+- [ ] Docstring explaining PID control
-  - [x] Mathematical formulation ✅
+- [ ] When to prefer PID over PI
-  - [x] When to use PID vs PI ✅
+- [ ] Coefficient selection guidance
-  - [x] Coefficient selection guidance ✅
+- [ ] Example comparing PI and PID behavior
-- [x] Usage examples in docstring ✅
-- [x] Comparison with PI in tests ✅
 
 ## Testing Requirements
 
@@ -239,15 +224,13 @@ Track standard deviation of log(h_i/h_{i-1}) over the integration:
 
 ## Success Criteria
 
-- [x] Implements full PID formula correctly ✅
+- [ ] Implements full PID formula correctly
-- [x] Handles first/second step bootstrap ✅
+- [ ] Handles first/second step bootstrap
-- [x] Shows similar stability on oscillatory test problem ✅
+- [ ] Shows improved stability on oscillatory test problem
-- [x] Performance similar to PI on smooth problems ✅
+- [ ] Performance similar to PI on smooth problems
-- [x] Error history management correct after rejections ✅
+- [ ] Error history management correct after rejections
-- [x] Documentation complete with usage examples ✅
+- [ ] Documentation complete with usage examples
-- [x] Coefficient sets match literature values ✅
+- [ ] Coefficient sets match literature values
-
-**STATUS**: ✅ **ALL SUCCESS CRITERIA MET**
 
 ## Future Enhancements
 

src/controller.rs

@@ -94,235 +94,12 @@ impl Default for PIController {
     }
 }
 
-/// PID (Proportional-Integral-Derivative) step size controller
-///
-/// The PID controller is an advanced adaptive time-stepping controller that provides
-/// better stability than the PI controller, especially for difficult or oscillatory problems.
-///
-/// # Mathematical Formulation
-///
-/// The PID controller determines the next step size based on error estimates from the
-/// current and previous two steps:
-///
-/// ```text
-/// h_{n+1} = h_n * safety * (ε_n)^(-β₁) * (ε_{n-1})^(-β₂) * (h_n/h_{n-1})^(-β₃)
-/// ```
-///
-/// 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)
-///
-/// # When to Use
-///
-/// Prefer PID over PI when:
-/// - Problem exhibits step size oscillation with PI controller
-/// - Working with stiff or near-stiff problems
-/// - Need smoother step size evolution
-/// - Standard in production solvers (MATLAB, Sundials)
-///
-/// # Example
-///
-/// ```ignore
-/// use ordinary_diffeq::prelude::*;
-///
-/// // Default PID controller (conservative coefficients)
-/// let controller = PIDController::default();
-///
-/// // Custom PID controller
-/// let controller = PIDController::new(0.49, 0.34, 0.10, 10.0, 0.2, 100.0, 0.9, 1e-4);
-///
-/// // PID tuned for specific method order (Gustafsson coefficients)
-/// let controller = PIDController::for_order(5);
-/// ```
-#[derive(Debug, Clone, Copy)]
-pub struct PIDController {
-    // PID Coefficients
-    pub beta1: f64,  // Proportional: reaction to current error
-    pub beta2: f64,  // Integral: reaction to error history
-    pub beta3: f64,  // Derivative: reaction to error rate of change
-
-    // 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)
-
-    // 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
-
-    // Next step guess
-    pub next_step_guess: TryStep,
-}
-
-impl<const D: usize> Controller<D> for PIDController {
-    /// Determines if the previously run step was acceptable and computes the next step size
-    /// using PID control theory
-    fn determine_step(&mut self, h: f64, err: f64) -> TryStep {
-        // 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);
-            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);
-            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);
-            // Derivative term uses ratio of consecutive step sizes
-            let factor_13 = if self.h_old > 0.0 {
-                (h / self.h_old).powf(-self.beta3)
-            } else {
-                1.0
-            };
-            self.factor_c2.max(
-                self.factor_c1.min(factor_11 * factor_12 * factor_13 / self.safety_factor)
-            )
-        };
-
-        if err <= 1.0 {
-            // Step accepted
-            let mut h_next = h / factor;
-
-            // 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;
-
-            // Apply maximum step size limit
-            if h_next.abs() > self.h_max {
-                h_next = self.h_max.copysign(h_next);
-            }
-
-            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 h_next = h / (self.factor_c1.min(factor_11 / self.safety_factor));
-
-            // Note: Don't update history on rejection
-            TryStep::NotYetAccepted(h_next)
-        }
-    }
-}
-
-impl PIDController {
-    /// Create a new PID controller with custom parameters
-    ///
-    /// # Arguments
-    ///
-    /// * `beta1` - Proportional coefficient (typically 0.3-0.5)
-    /// * `beta2` - Integral coefficient (typically 0.04-0.1)
-    /// * `beta3` - Derivative coefficient (typically 0.01-0.05)
-    /// * `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
-    pub fn new(
-        beta1: f64,
-        beta2: f64,
-        beta3: f64,
-        max_factor: f64,
-        min_factor: f64,
-        h_max: f64,
-        safety_factor: f64,
-        initial_h: f64,
-    ) -> Self {
-        Self {
-            beta1,
-            beta2,
-            beta3,
-            factor_c1: 1.0 / min_factor,
-            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
-            next_step_guess: TryStep::NotYetAccepted(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)
-    ///
-    /// # Arguments
-    ///
-    /// * `order` - Order of the integration method (e.g., 5 for DP5, 7 for Vern7)
-    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
-            10.0,              // max_factor
-            0.2,               // min_factor
-            100000.0,          // h_max
-            0.9,               // safety_factor
-            1e-4,              // initial_h
-        )
-    }
-
-    /// Reset the controller's error history
-    ///
-    /// Useful when switching algorithms or restarting integration
-    pub fn reset(&mut self) {
-        self.err_old = 0.0;
-        self.err_older = 0.0;
-        self.h_old = 0.0;
-    }
-}
-
-impl Default for PIDController {
-    /// Default PID controller using conservative coefficients
-    ///
-    /// Uses conservative PID coefficients that provide stable performance
-    /// across a wide range of problems:
-    /// - β₁ = 0.07 (proportional)
-    /// - β₂ = 0.04 (integral)
-    /// - β₃ = 0.01 (derivative)
-    ///
-    /// 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)
-            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::*;
 
     #[test]
-    fn test_pi_controller_creation() {
+    fn test_controller_creation() {
         let controller = PIController::new(0.17, 0.04, 10.0, 0.2, 10.0, 0.9, 1e-4);
 
         assert!(controller.alpha == 0.17);
@@ -334,229 +111,4 @@ mod tests {
         assert!(controller.safety_factor == 0.9);
         assert!(controller.next_step_guess == TryStep::NotYetAccepted(1e-4));
     }
-
-    #[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);
-
-        assert_eq!(controller.beta1, 0.49);
-        assert_eq!(controller.beta2, 0.34);
-        assert_eq!(controller.beta3, 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);
-        assert_eq!(controller.safety_factor, 0.9);
-        assert_eq!(controller.err_old, 0.0);
-        assert_eq!(controller.err_older, 0.0);
-        assert_eq!(controller.h_old, 0.0);
-    }
-
-    #[test]
-    fn test_pid_for_order() {
-        let controller = PIDController::for_order(5);
-
-        // For order 5, k+1 = 6
-        assert!((controller.beta1 - 0.49 / 6.0).abs() < 1e-10);
-        assert!((controller.beta2 - 0.34 / 6.0).abs() < 1e-10);
-        assert!((controller.beta3 - 0.10 / 6.0).abs() < 1e-10);
-    }
-
-    #[test]
-    fn test_pid_default() {
-        let controller = PIDController::default();
-
-        // Default uses conservative PID coefficients
-        assert_eq!(controller.beta1, 0.07);
-        assert_eq!(controller.beta2, 0.04);
-        assert_eq!(controller.beta3, 0.01);  // True PID with derivative term
-    }
-
-    #[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;
-
-        controller.reset();
-
-        assert_eq!(controller.err_old, 0.0);
-        assert_eq!(controller.err_older, 0.0);
-        assert_eq!(controller.h_old, 0.0);
-    }
-
-    #[test]
-    fn test_pi_vs_pid_smooth_problem() {
-        // 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();
-
-        // 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];
-        let h = 0.01;
-
-        let mut pi_steps = Vec::new();
-        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);
-
-            if pi_result.is_accepted() {
-                pi_steps.push(pi_result.extract());
-                pi.next_step_guess = pi_result.reset().unwrap();
-            }
-
-            if pid_result.is_accepted() {
-                pid_steps.push(pid_result.extract());
-                pid.next_step_guess = pid_result.reset().unwrap();
-            }
-        }
-
-        // Both should accept all steps for this smooth sequence
-        assert_eq!(pi_steps.len(), errors.len());
-        assert_eq!(pid_steps.len(), errors.len());
-
-        // Step sizes should be reasonably similar (within 20%)
-        // PID may differ slightly due to derivative term
-        for (pi_h, pid_h) in pi_steps.iter().zip(pid_steps.iter()) {
-            let relative_diff = ((pi_h - pid_h) / pi_h).abs();
-            assert!(
-                relative_diff < 0.2,
-                "Step sizes differ by more than 20%: PI={}, PID={}",
-                pi_h,
-                pid_h
-            );
-        }
-    }
-
-    #[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 h = 0.01;
-        let err1 = 0.5;
-        let err2 = 0.4;
-        let err3 = 0.3;
-
-        // First step: should use only proportional (beta1)
-        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);
-        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);
-        assert!(step2.is_accepted());
-
-        // After second step, both err_old and err_older should be set
-        assert!(pid.err_old > 0.0);
-        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);
-        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 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);
-        assert!(step1.is_accepted());
-
-        let err_old_before = pid.err_old;
-        let err_older_before = pid.err_older;
-        let h_old_before = pid.h_old;
-
-        // 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);
-        assert!(!step2.is_accepted());
-
-        // History should NOT have changed after rejection
-        assert_eq!(pid.err_old, err_old_before);
-        assert_eq!(pid.err_older, err_older_before);
-        assert_eq!(pid.h_old, h_old_before);
-    }
-
-    #[test]
-    fn test_pid_vs_pi_oscillatory() {
-        // 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);
-
-        // Simulate oscillatory error pattern
-        let errors = vec![0.8, 0.3, 0.9, 0.2, 0.85, 0.25, 0.8, 0.3];
-        let h = 0.01;
-
-        let mut pi_ratios = Vec::new();
-        let mut pid_ratios = Vec::new();
-
-        let mut pi_h_prev = h;
-        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);
-
-            if pi_result.is_accepted() {
-                let pi_h_next = pi_result.reset().unwrap().extract();
-                pi_ratios.push((pi_h_next / pi_h_prev).ln().abs());
-                pi_h_prev = pi_h_next;
-                pi.next_step_guess = TryStep::NotYetAccepted(pi_h_next);
-            }
-
-            if pid_result.is_accepted() {
-                let pid_h_next = pid_result.reset().unwrap().extract();
-                pid_ratios.push((pid_h_next / pid_h_prev).ln().abs());
-                pid_h_prev = pid_h_next;
-                pid.next_step_guess = TryStep::NotYetAccepted(pid_h_next);
-            }
-        }
-
-        // 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_std = pi_variance.sqrt();
-
-        let pid_mean: f64 = pid_ratios.iter().sum::<f64>() / pid_ratios.len() as f64;
-        let pid_variance: f64 = pid_ratios
-            .iter()
-            .map(|r| (r - pid_mean).powi(2))
-            .sum::<f64>()
-            / pid_ratios.len() as f64;
-        let pid_std = pid_variance.sqrt();
-
-        // With true PID (non-zero beta3), we expect similar or better stability
-        // Allow some tolerance since this is a simple synthetic test
-        assert!(
-            pid_std <= pi_std * 1.1,
-            "PID should not be significantly worse than PI: PI_std={:.3}, PID_std={:.3}",
-            pi_std,
-            pid_std
-        );
-    }
 }

src/lib.rs

@@ -8,7 +8,7 @@ pub mod problem;
 
 pub mod prelude {
     pub use super::callback::{stop, Callback};
-    pub use super::controller::{PIController, PIDController};
+    pub use super::controller::PIController;
     pub use super::integrator::bs3::BS3;
     pub use super::integrator::dormand_prince::DormandPrince45;
     pub use super::integrator::vern7::Vern7;