diff --git a/roadmap/FEATURE_TEMPLATES.md b/roadmap/FEATURE_TEMPLATES.md new file mode 100644 index 0000000..aacfa5c --- /dev/null +++ b/roadmap/FEATURE_TEMPLATES.md @@ -0,0 +1,237 @@ +# Feature File Templates + +This document contains brief summaries for features 6-38. Detailed feature files should be created when you're ready to implement each one, using the detailed examples in features 01-05 and 12 as templates. + +## How to Use This Document + +When ready to implement a feature: +1. Copy the template structure from features/01-bs3-method.md or similar +2. Fill in the details from the summary below +3. Add implementation-specific details +4. Create comprehensive testing requirements + +--- + +## Feature 06: CallbackSet +**Description**: Compose multiple callbacks with ordering +**Dependencies**: Discrete callbacks +**Effort**: Small +**Key Points**: Builder pattern, execution priority, enable/disable + +## Feature 07: Saveat Functionality +**Description**: Save solution at specific timepoints +**Dependencies**: None +**Effort**: Medium +**Key Points**: Interpolation to exact times, dense vs sparse saving, memory efficiency + +## Feature 08: Solution Derivatives +**Description**: Access derivatives at any time via interpolation +**Dependencies**: None +**Effort**: Small +**Key Points**: `solution.derivative(t)` interface, use dense output or finite differences + +## Feature 09: DP8 Method +**Description**: Dormand-Prince 8th order method +**Dependencies**: None +**Effort**: Medium +**Key Points**: 13 stages, very high accuracy, tableau from literature + +## Feature 10: FBDF Method +**Description**: Fixed-leading-coefficient BDF multistep method +**Dependencies**: Linear solver, Nordsieck representation +**Effort**: Large +**Key Points**: Variable order (1-5), excellent for very stiff problems, complex state management + +## Feature 11: Rodas4/Rodas5P +**Description**: Higher-order Rosenbrock methods +**Dependencies**: Rosenbrock23 +**Effort**: Medium +**Key Points**: 4th/5th order accuracy, more stages, better for higher accuracy stiff problems + +## Feature 13: Default Algorithm Selection +**Description**: Smart defaults based on problem characteristics +**Dependencies**: Auto-switching, multiple algorithms +**Effort**: Medium +**Key Points**: Analyze tolerance, problem size, choose algorithms automatically + +## Feature 14: Automatic Initial Stepsize +**Description**: Algorithm to compute good initial dt +**Dependencies**: None +**Effort**: Medium +**Key Points**: Based on Hairer & Wanner algorithm, uses local Lipschitz estimate + +## Feature 15: PresetTimeCallback +**Description**: Callbacks at predetermined times +**Dependencies**: Discrete callbacks +**Effort**: Small +**Key Points**: Efficient time-based events, integration with tstops + +## Feature 16: TerminateSteadyState +**Description**: Auto-detect when solution reaches steady state +**Dependencies**: Discrete callbacks +**Effort**: Small +**Key Points**: Monitor du/dt, terminate when small enough + +## Feature 17: SavingCallback +**Description**: Custom saving logic beyond default +**Dependencies**: CallbackSet +**Effort**: Small +**Key Points**: User-defined save conditions, memory-efficient for large problems + +## Feature 18: Linear Solver Infrastructure +**Description**: Generic linear solver interface and dense LU +**Dependencies**: None +**Effort**: Large +**Key Points**: +- Trait-based design for flexibility +- Dense LU factorization with partial pivoting +- Solve Ax = b efficiently +- Foundation for all implicit methods +- Consider using nalgebra's built-in LU or implement custom + +## Feature 19: Jacobian Computation +**Description**: Finite difference and auto-diff Jacobians +**Dependencies**: None +**Effort**: Large +**Key Points**: +- Forward finite differences: (f(y+εe_j) - f(y))/ε +- Epsilon selection: √eps * max(|y_j|, 1) +- Sparse Jacobian support (future) +- Integration with AD crates (future) + +## Feature 20: Low-Storage Runge-Kutta +**Description**: 2N/3N/4N storage variants for large systems +**Dependencies**: None +**Effort**: Medium +**Key Points**: Specialized RK methods that reuse storage, critical for PDEs via method-of-lines + +## Feature 21: SSP Methods +**Description**: Strong Stability Preserving RK methods +**Dependencies**: None +**Effort**: Medium +**Key Points**: SSPRK22, SSPRK33, SSPRK43, SSPRK53, preserve TVD/monotonicity, for hyperbolic PDEs + +## Feature 22: Symplectic Integrators +**Description**: Verlet, Leapfrog, KahanLi for Hamiltonian systems +**Dependencies**: None (second-order ODE support already exists) +**Effort**: Medium +**Key Points**: Preserve energy/symplectic structure, special for p,q formulation + +## Feature 23: Verner Methods Suite +**Description**: Complete Verner family (Vern6, Vern8, Vern9) +**Dependencies**: Vern7 +**Effort**: Medium +**Key Points**: Different orders for different accuracy needs, all highly efficient + +## Feature 24: SDIRK Methods +**Description**: Singly Diagonally Implicit RK (KenCarp3/4/5) +**Dependencies**: Linear solver, nonlinear solver +**Effort**: Large +**Key Points**: IMEX methods, good for semi-stiff problems, L-stable + +## Feature 25: Exponential Integrators +**Description**: Exp4, EPIRK4, EXPRB53 for semi-linear stiff +**Dependencies**: Matrix exponential computation +**Effort**: Large +**Key Points**: For du/dt = Lu + N(u), where L is linear stiff part + +## Feature 26: Extrapolation Methods +**Description**: Richardson extrapolation with adaptive order +**Dependencies**: Linear solver +**Effort**: Large +**Key Points**: High accuracy from low-order methods, variable order selection + +## Feature 27: Stabilized Methods +**Description**: ROCK2, ROCK4, RKC for mildly stiff +**Dependencies**: None +**Effort**: Medium +**Key Points**: Extended stability regions, good for PDEs with explicit time-stepping + +## Feature 28: I Controller +**Description**: Basic integral controller +**Dependencies**: None +**Effort**: Small +**Key Points**: Simplest adaptive controller, mainly for comparison/testing + +## Feature 29: Predictive Controller +**Description**: Advanced predictive step size control +**Dependencies**: None +**Effort**: Medium +**Key Points**: Predicts future error, more sophisticated than PID + +## Feature 30: VectorContinuousCallback +**Description**: Multiple simultaneous event detection +**Dependencies**: CallbackSet +**Effort**: Medium +**Key Points**: More efficient than separate callbacks, shared root-finding + +## Feature 31: PositiveDomain +**Description**: Enforce positivity constraints +**Dependencies**: CallbackSet +**Effort**: Small +**Key Points**: Ensures solution stays positive, important for physical systems + +## Feature 32: ManifoldProjection +**Description**: Project solution onto constraint manifolds +**Dependencies**: CallbackSet +**Effort**: Medium +**Key Points**: For constrained mechanical systems, projection step after integration + +## Feature 33: Nonlinear Solver Infrastructure +**Description**: Newton and quasi-Newton methods +**Dependencies**: Linear solver +**Effort**: Large +**Key Points**: +- Newton's method for implicit stages +- Line search +- Convergence criteria +- Foundation for SDIRK, FIRK methods + +## Feature 34: Krylov Linear Solvers +**Description**: GMRES, BiCGStab for large sparse systems +**Dependencies**: Linear solver infrastructure +**Effort**: Large +**Key Points**: Iterative solvers for when LU factorization too expensive + +## Feature 35: Preconditioners +**Description**: ILU, Jacobi, custom preconditioners +**Dependencies**: Krylov solvers +**Effort**: Large +**Key Points**: Accelerate Krylov methods, essential for large sparse systems + +## Feature 36: FSAL Optimization +**Description**: First-Same-As-Last function evaluation reuse +**Dependencies**: None +**Effort**: Small +**Key Points**: Reduce function evaluations by ~14% for FSAL methods (DP5, Tsit5, etc.) + +## Feature 37: Custom Norms +**Description**: User-definable error norms +**Dependencies**: None +**Effort**: Small +**Key Points**: L2, Linf, weighted norms, custom user functions + +## Feature 38: Step/Stage Limiting +**Description**: Limit state values during integration +**Dependencies**: None +**Effort**: Small +**Key Points**: Enforce bounds on solution, prevent non-physical values + +--- + +## Creating Detailed Feature Files + +When you're ready to work on a feature, create a detailed file following this structure: + +1. **Overview**: What is it, key characteristics +2. **Why This Feature Matters**: Motivation, use cases +3. **Dependencies**: What must be built first +4. **Implementation Approach**: Algorithm details, design decisions +5. **Implementation Tasks**: Detailed checklist with subtasks +6. **Testing Requirements**: Specific tests with expected results +7. **References**: Papers, Julia code, textbooks +8. **Complexity Estimate**: Effort and risk assessment +9. **Success Criteria**: How to know it's done right +10. **Future Enhancements**: What could be added later + +See `features/01-bs3-method.md`, `features/02-vern7-method.md`, `features/03-rosenbrock23.md`, etc. for complete examples. diff --git a/roadmap/GETTING_STARTED.md b/roadmap/GETTING_STARTED.md new file mode 100644 index 0000000..1f5bf50 --- /dev/null +++ b/roadmap/GETTING_STARTED.md @@ -0,0 +1,249 @@ +# Getting Started with the Roadmap + +This guide helps you navigate the development roadmap for the Rust ODE library. + +## Roadmap Structure + +``` +roadmap/ +├── README.md # Master overview with all features +├── GETTING_STARTED.md # This file +├── FEATURE_TEMPLATES.md # Brief summaries of features 6-38 +└── features/ + ├── 01-bs3-method.md # Detailed implementation plan (example) + ├── 02-vern7-method.md # Detailed implementation plan + ├── 03-rosenbrock23.md # Detailed implementation plan + ├── 04-pid-controller.md # Detailed implementation plan + ├── 05-discrete-callbacks.md # Detailed implementation plan + ├── 06-callback-set.md # Brief outline + └── 12-auto-switching.md # Detailed implementation plan + └── ... (create detailed files as needed) +``` + +## How to Use This Roadmap + +### 1. Review the Master Plan + +Start with `README.md` to see: +- All 38 planned features organized by tier +- Dependencies between features +- Current completion status +- Overall progress tracking + +### 2. Choose Your Next Feature + +**Recommended Order for Beginners:** +1. Start with Tier 1 features (essential) +2. Follow dependency chains +3. Mix difficulty levels (alternate hard and easy) + +**Suggested First 5 Features:** +1. **BS3 Method** (feature #1) - Easy, builds confidence +2. **PID Controller** (feature #4) - Easy, immediate value +3. **Discrete Callbacks** (feature #5) - Easy, useful capability +4. **Vern7** (feature #2) - Medium, important algorithm +5. **Linear Solver Infrastructure** (feature #18) - Hard but foundational + +### 3. Read the Detailed Feature File + +Each detailed feature file contains: +- **Overview**: Quick introduction +- **Why It Matters**: Motivation +- **Dependencies**: What you need first +- **Implementation Approach**: Algorithm details +- **Implementation Tasks**: Detailed checklist +- **Testing Requirements**: How to verify it works +- **References**: Where to learn more +- **Complexity Estimate**: Time and difficulty +- **Success Criteria**: Definition of done + +### 4. Implement the Feature + +Follow the detailed task checklist: +- [ ] Read references and understand algorithm +- [ ] Implement core algorithm +- [ ] Write tests +- [ ] Document +- [ ] Benchmark +- [ ] Check off tasks as you complete them + +### 5. Update the Roadmap + +When you complete a feature: +1. Check the box in `README.md` +2. Update completion statistics +3. Note any lessons learned or deviations from plan + +## Current State (Baseline) + +Your library already has: +- ✅ Dormand-Prince 4(5) with dense output +- ✅ Tsit5 with dense output +- ✅ PI Controller +- ✅ Continuous callbacks with zero-crossing detection +- ✅ Solution interpolation interface +- ✅ Generic over compile-time array dimensions +- ✅ Support for second-order ODE problems + +This is a solid foundation! The roadmap builds on this. + +## Recommended Development Path + +### Phase 1: Core Algorithm Diversity (Tier 1) +*Goal: Give users algorithm choices* + +1. BS3 - Easy, quick win +2. Vern7 - High accuracy option +3. Build linear solver infrastructure +4. Rosenbrock23 - First stiff solver +5. PID Controller - Better adaptive stepping +6. Discrete Callbacks - More event types + +**Milestone**: Can handle both non-stiff and stiff problems efficiently. + +### Phase 2: Robustness & Automation (Tier 2) +*Goal: Make the library production-ready* + +7. Auto-switching/stiffness detection +8. Automatic initial step size +9. More Rosenbrock methods (Rodas4) +10. BDF method +11. CallbackSet and advanced callbacks +12. Saveat functionality + +**Milestone**: Library can solve most problems automatically with minimal user input. + +### Phase 3: Specialization & Performance (Tier 3) +*Goal: Optimize for specific problem classes* + +13. Low-storage RK for large systems +14. Symplectic integrators for Hamiltonian systems +15. SSP methods for hyperbolic PDEs +16. Verner suite completion +17. Advanced linear/nonlinear solvers +18. Performance optimizations (FSAL, custom norms) + +**Milestone**: Best-in-class performance for specialized problem types. + +## Development Tips + +### Testing Strategy + +Every feature should have: +1. **Convergence test**: Verify order of accuracy +2. **Correctness test**: Compare to known solutions +3. **Edge case tests**: Boundary conditions, error handling +4. **Benchmark**: Performance measurement + +### Reference Material + +When implementing a feature: +1. Read the Julia implementation for guidance +2. Check original papers for algorithm details +3. Verify tableau/coefficients from authoritative sources +4. Test against reference solutions from DiffEqDevDocs + +### Common Pitfalls + +- **Don't skip testing**: Numerical bugs are subtle +- **Verify tableau coefficients**: Transcription errors are common +- **Check interpolation**: Easy to get wrong +- **Test stiff problems**: If implementing stiff solvers +- **Benchmark early**: Performance problems easier to fix early + +## Getting Help + +### Resources + +1. **Julia's OrdinaryDiffEq.jl**: Reference implementation + - Location: `/tmp/diffeq_copy/OrdinaryDiffEq.jl/` + - Well-tested, can compare behavior + +2. **Hairer & Wanner textbooks**: + - "Solving ODEs I: Nonstiff Problems" + - "Solving ODEs II: Stiff and DAE Problems" + +3. **DiffEqDevDocs**: Developer documentation + - https://docs.sciml.ai/DiffEqDevDocs/stable/ + +4. **Test problems**: Standard ODE test suite + - Van der Pol, Robertson, Pleiades, etc. + - Reference solutions available + +### Creating New Detailed Feature Files + +When ready to work on a feature that only has a brief summary: + +1. Copy structure from `features/01-bs3-method.md` +2. Fill in details from `FEATURE_TEMPLATES.md` +3. Add algorithm-specific information +4. Create comprehensive task checklist +5. Define specific test requirements +6. Estimate complexity honestly + +## Tracking Progress + +### In README.md + +Update the checkboxes as features are completed: +- [ ] Incomplete +- [x] Complete + +Update completion statistics at bottom: +``` +## Progress Tracking + +Total Features: 38 +- Tier 1: 8 features (3/8 complete) # Update these +- Tier 2: 12 features (0/12 complete) +- Tier 3: 18 features (0/18 complete) +``` + +### Optional: Keep a CHANGELOG.md + +Document major milestones: +```markdown +# Changelog + +## 2025-01-XX +- Completed BS3 method +- Completed PID controller +- Started Vern7 implementation + +## 2025-01-YY +- Completed Vern7 +- Linear solver infrastructure in progress +``` + +## Questions to Ask Before Starting + +Before implementing a feature: + +1. **Do I understand the algorithm?** + - Read the papers + - Understand the math + - Know the use cases + +2. **Are dependencies satisfied?** + - Check the dependency list + - Make sure infrastructure exists + +3. **Do I have test cases ready?** + - Know how to verify correctness + - Have reference solutions + +4. **What's the success criteria?** + - Clear definition of "done" + - Performance targets + +## Next Steps + +1. Read `README.md` to see the full roadmap +2. Pick a feature to start with (suggest: BS3 or PID Controller) +3. Read its detailed feature file +4. Implement following the task checklist +5. Test thoroughly +6. Update the roadmap +7. Move to next feature! + +Good luck! You're building something great. 🚀 diff --git a/roadmap/README.md b/roadmap/README.md new file mode 100644 index 0000000..089e541 --- /dev/null +++ b/roadmap/README.md @@ -0,0 +1,334 @@ +# Ordinary Differential Equations Library - Development Roadmap + +This roadmap outlines the planned features for developing a comprehensive Rust-based ODE solver library, inspired by Julia's OrdinaryDiffEq.jl but adapted for Rust's strengths and idioms. + +## Current Foundation + +The library currently has: +- ✅ Dormand-Prince 4(5) adaptive explicit RK method with dense output +- ✅ Tsit5 explicit RK method with dense output +- ✅ PI step size controller +- ✅ Basic continuous callbacks with zero-crossing detection +- ✅ Solution interpolation interface +- ✅ Generic over array dimensions at compile-time +- ✅ Support for ordinary and second-order ODE problems + +## Roadmap Organization + +Features are organized into three tiers based on priority and dependencies: +- **Tier 1**: Essential features for general-purpose ODE solving +- **Tier 2**: Important features for robustness and broader applicability +- **Tier 3**: Advanced/specialized features for specific problem classes + +Each feature below links to a detailed implementation plan in the `features/` directory. + +--- + +## Tier 1: Essential Features + +### Algorithms + +- [ ] **[BS3 (Bogacki-Shampine 3/2)](features/01-bs3-method.md)** + - 3rd order explicit RK method with 2nd order error estimate + - Good for moderate accuracy, lower cost than DP5 + - **Dependencies**: None + - **Effort**: Small + +- [ ] **[Vern7 (Verner 7th order)](features/02-vern7-method.md)** + - 7th order explicit RK method for high-accuracy non-stiff problems + - Efficient for tight tolerances + - **Dependencies**: None + - **Effort**: Medium + +- [ ] **[Rosenbrock23](features/03-rosenbrock23.md)** + - L-stable 2nd/3rd order Rosenbrock-W method + - First working stiff solver + - **Dependencies**: Linear solver infrastructure, Jacobian computation + - **Effort**: Large + +### Controllers + +- [ ] **[PID Controller](features/04-pid-controller.md)** + - Proportional-Integral-Derivative step size controller + - Better stability than PI controller for difficult problems + - **Dependencies**: None + - **Effort**: Small + +### Callbacks + +- [ ] **[Discrete Callbacks](features/05-discrete-callbacks.md)** + - Event detection based on conditions (not zero-crossings) + - Useful for time-based events, iteration counts, etc. + - **Dependencies**: None + - **Effort**: Small + +- [ ] **[CallbackSet](features/06-callback-set.md)** + - Compose multiple callbacks with ordering + - Essential for complex simulations + - **Dependencies**: Discrete callbacks + - **Effort**: Small + +### Solution Interface + +- [ ] **[Saveat Functionality](features/07-saveat.md)** + - Save solution at specific timepoints + - Dense vs sparse saving strategies + - **Dependencies**: None + - **Effort**: Medium + +- [ ] **[Solution Derivatives](features/08-solution-derivatives.md)** + - Access derivatives at any time point via interpolation + - `solution.derivative(t)` interface + - **Dependencies**: None + - **Effort**: Small + +--- + +## Tier 2: Important for Robustness + +### Algorithms + +- [ ] **[DP8 (Dormand-Prince 8th order)](features/09-dp8-method.md)** + - 8th order explicit RK for very high accuracy + - Complements Vern7 for algorithm selection + - **Dependencies**: None + - **Effort**: Medium + +- [ ] **[FBDF (Fixed-leading-coefficient BDF)](features/10-fbdf-method.md)** + - Multistep method for very stiff problems + - More robust than Rosenbrock for some problem classes + - **Dependencies**: Linear solver infrastructure, Nordsieck representation + - **Effort**: Large + +- [ ] **[Rodas4/Rodas5P](features/11-rodas-methods.md)** + - Higher-order Rosenbrock methods (4th/5th order) + - Better accuracy for stiff problems + - **Dependencies**: Rosenbrock23 + - **Effort**: Medium + +### Algorithm Selection + +- [ ] **[Auto-Switching / Stiffness Detection](features/12-auto-switching.md)** + - Automatic detection of stiffness + - Switch between non-stiff and stiff solvers + - Composite algorithm infrastructure + - **Dependencies**: At least one stiff solver (Rosenbrock23 or FBDF) + - **Effort**: Large + +- [ ] **[Default Algorithm Selection](features/13-default-algorithm.md)** + - Smart defaults based on problem characteristics + - `solve(problem)` without specifying algorithm + - **Dependencies**: Auto-switching, multiple algorithms + - **Effort**: Medium + +### Initialization + +- [ ] **[Automatic Initial Step Size](features/14-initial-stepsize.md)** + - Algorithm to determine good initial dt + - Based on local Lipschitz estimate + - **Dependencies**: None + - **Effort**: Medium + +### Callbacks + +- [ ] **[PresetTimeCallback](features/15-preset-time-callback.md)** + - Trigger callbacks at specific predetermined times + - Important for time-varying forcing functions + - **Dependencies**: Discrete callbacks + - **Effort**: Small + +- [ ] **[TerminateSteadyState](features/16-terminate-steady-state.md)** + - Auto-detect when solution reaches steady state + - Stop integration early + - **Dependencies**: Discrete callbacks + - **Effort**: Small + +- [ ] **[SavingCallback](features/17-saving-callback.md)** + - Custom saving logic beyond default + - For memory-efficient large-scale simulations + - **Dependencies**: CallbackSet + - **Effort**: Small + +### Infrastructure + +- [ ] **[Linear Solver Infrastructure](features/18-linear-solver-infrastructure.md)** + - Generic linear solver interface + - Dense LU factorization + - Foundation for implicit methods + - **Dependencies**: None + - **Effort**: Large + +- [ ] **[Jacobian Computation](features/19-jacobian-computation.md)** + - Finite difference Jacobians + - Forward-mode automatic differentiation + - Sparse Jacobian support + - **Dependencies**: None + - **Effort**: Large + +--- + +## Tier 3: Advanced & Specialized + +### Algorithms + +- [ ] **[Low-Storage Runge-Kutta](features/20-low-storage-rk.md)** + - 2N, 3N, 4N storage variants + - Critical for large-scale problems + - **Dependencies**: None + - **Effort**: Medium + +- [ ] **[SSP (Strong Stability Preserving) Methods](features/21-ssp-methods.md)** + - SSPRK22, SSPRK33, SSPRK43, SSPRK53, etc. + - For hyperbolic PDEs and method-of-lines + - **Dependencies**: None + - **Effort**: Medium + +- [ ] **[Symplectic Integrators](features/22-symplectic-integrators.md)** + - Velocity Verlet, Leapfrog, KahanLi6, KahanLi8 + - For Hamiltonian systems, preserves energy + - **Dependencies**: Second-order ODE infrastructure (already exists) + - **Effort**: Medium + +- [ ] **[Verner Methods Suite](features/23-verner-methods.md)** + - Vern6, Vern8, Vern9 + - Complete Verner family for different accuracy needs + - **Dependencies**: Vern7 + - **Effort**: Medium + +- [ ] **[SDIRK Methods](features/24-sdirk-methods.md)** + - KenCarp3, KenCarp4, KenCarp5 + - Singly Diagonally Implicit RK for stiff problems + - **Dependencies**: Linear solver infrastructure, nonlinear solver + - **Effort**: Large + +- [ ] **[Exponential Integrators](features/25-exponential-integrators.md)** + - Exp4, EPIRK4, EXPRB53 + - For semi-linear stiff problems + - **Dependencies**: Matrix exponential computation + - **Effort**: Large + +- [ ] **[Extrapolation Methods](features/26-extrapolation-methods.md)** + - Implicit Euler/Midpoint with Richardson extrapolation + - Adaptive order selection + - **Dependencies**: Linear solver infrastructure + - **Effort**: Large + +- [ ] **[Stabilized Methods](features/27-stabilized-methods.md)** + - ROCK2, ROCK4, RKC + - For mildly stiff problems with large systems + - **Dependencies**: None + - **Effort**: Medium + +### Controllers + +- [ ] **[I Controller](features/28-i-controller.md)** + - Basic integral controller + - For completeness and testing + - **Dependencies**: None + - **Effort**: Small + +- [ ] **[Predictive Controller](features/29-predictive-controller.md)** + - Advanced controller with prediction + - For challenging adaptive stepping scenarios + - **Dependencies**: None + - **Effort**: Medium + +### Advanced Callbacks + +- [ ] **[VectorContinuousCallback](features/30-vector-continuous-callback.md)** + - Multiple simultaneous event detection + - More efficient than multiple callbacks + - **Dependencies**: CallbackSet + - **Effort**: Medium + +- [ ] **[PositiveDomain](features/31-positive-domain.md)** + - Enforce positivity constraints + - Important for physical systems + - **Dependencies**: CallbackSet + - **Effort**: Small + +- [ ] **[ManifoldProjection](features/32-manifold-projection.md)** + - Project solution onto constraint manifolds + - For constrained mechanical systems + - **Dependencies**: CallbackSet + - **Effort**: Medium + +### Infrastructure + +- [ ] **[Nonlinear Solver Infrastructure](features/33-nonlinear-solver.md)** + - Newton's method + - Quasi-Newton methods + - Generic nonlinear solver interface + - **Dependencies**: Linear solver infrastructure + - **Effort**: Large + +- [ ] **[Krylov Linear Solvers](features/34-krylov-solvers.md)** + - GMRES, BiCGStab + - For large sparse systems + - **Dependencies**: Linear solver infrastructure + - **Effort**: Large + +- [ ] **[Preconditioners](features/35-preconditioners.md)** + - ILU, Jacobi, custom preconditioners + - Accelerate Krylov methods + - **Dependencies**: Krylov solvers + - **Effort**: Large + +### Performance & Optimization + +- [ ] **[FSAL Optimization](features/36-fsal-optimization.md)** + - First-Same-As-Last reuse + - Reduce function evaluations + - **Dependencies**: None + - **Effort**: Small + +- [ ] **[Custom Norms](features/37-custom-norms.md)** + - User-definable error norms + - L2, Linf, weighted norms + - **Dependencies**: None + - **Effort**: Small + +- [ ] **[Step/Stage Limiting](features/38-step-stage-limiting.md)** + - Limit state values during integration + - For bounded problems + - **Dependencies**: None + - **Effort**: Small + +--- + +## Implementation Notes + +### General Principles + +1. **Rust-first design**: Leverage Rust's type system, zero-cost abstractions, and safety guarantees +2. **Compile-time optimization**: Use const generics for array sizes where beneficial +3. **Trait-based abstraction**: Generic over array types, number types, and algorithm components +4. **Comprehensive testing**: Each feature needs convergence tests and comparison to known solutions +5. **Benchmarking**: Track performance as features are added + +### Testing Strategy + +Each algorithm implementation should include: +- **Convergence tests**: Verify order of accuracy +- **Correctness tests**: Compare to analytical solutions +- **Stiffness tests**: For stiff solvers, test on Van der Pol, Robertson, etc. +- **Callback tests**: Verify event detection accuracy +- **Regression tests**: Prevent performance degradation + +### Documentation Requirements + +- Public API documentation +- Algorithm descriptions and references +- Usage examples +- Performance characteristics + +--- + +## Progress Tracking + +Total Features: 38 +- Tier 1: 8 features (0/8 complete) +- Tier 2: 12 features (0/12 complete) +- Tier 3: 18 features (0/18 complete) + +Last updated: 2025-01-XX diff --git a/roadmap/features/01-bs3-method.md b/roadmap/features/01-bs3-method.md new file mode 100644 index 0000000..1a0d5b5 --- /dev/null +++ b/roadmap/features/01-bs3-method.md @@ -0,0 +1,191 @@ +# Feature: BS3 (Bogacki-Shampine 3/2) Method + +## Overview + +The Bogacki-Shampine 3/2 method is a 3rd order explicit Runge-Kutta method with an embedded 2nd order method for error estimation. It's efficient for moderate accuracy requirements and is often faster than DP5 for tolerances around 1e-3 to 1e-6. + +**Key Characteristics:** +- Order: 3(2) - 3rd order solution with 2nd order error estimate +- Stages: 4 +- FSAL: Yes (First Same As Last) +- Adaptive: Yes +- Dense output: 3rd order continuous extension + +## Why This Feature Matters + +- **Efficiency**: Fewer stages than DP5 (4 vs 7) for comparable accuracy at moderate tolerances +- **Common use case**: Many practical problems don't need DP5's accuracy +- **Algorithm diversity**: Gives users choice based on problem characteristics +- **Foundation**: Good reference implementation for adding more RK methods + +## Dependencies + +- None (can be implemented with current infrastructure) + +## Implementation Approach + +### Butcher Tableau + +The BS3 method uses the following coefficients: + +``` +c | A +--+------- +0 | 0 +1/2 | 1/2 +3/4 | 0 3/4 +1 | 2/9 1/3 4/9 +--+------- +b | 2/9 1/3 4/9 0 (3rd order) +b*| 7/24 1/4 1/3 1/8 (2nd order, for error) +``` + +FSAL property: The last stage k4 can be reused as k1 of the next step. + +### Dense Output + +3rd order Hermite interpolation: +``` +u(t₀ + θh) = u₀ + h*θ*(b₁*k₁ + b₂*k₂ + b₃*k₃) + h*θ*(1-θ)*(...additional terms) +``` + +Coefficients from Bogacki & Shampine 1989 paper. + +### Error Estimation + +``` +err = ||u₃ - u₂|| / (atol + max(|u_n|, |u_{n+1}|) * rtol) +``` + +Where u₃ is the 3rd order solution and u₂ is the 2nd order embedded solution. + +## Implementation Tasks + +### Core Algorithm + +- [ ] Define `BS3` struct implementing `Integrator` trait + - [ ] Add tableau constants (A, b, b_error, c) + - [ ] Add tolerance fields (a_tol, r_tol) + - [ ] Add builder methods for setting tolerances + +- [ ] Implement `step()` method + - [ ] Compute k1 = f(t, y) + - [ ] Compute k2 = f(t + c[1]*h, y + h*a[0,0]*k1) + - [ ] Compute k3 = f(t + c[2]*h, y + h*(a[1,0]*k1 + a[1,1]*k2)) + - [ ] Compute k4 = f(t + c[3]*h, y + h*(a[2,0]*k1 + a[2,1]*k2 + a[2,2]*k3)) + - [ ] Compute 3rd order solution: y_next = y + h*(b[0]*k1 + b[1]*k2 + b[2]*k3 + b[3]*k4) + - [ ] Compute error estimate: err = h*(b[0]-b*[0])*k1 + ... (for all ki) + - [ ] Store dense output coefficients [k1, k2, k3, k4] + - [ ] Return (y_next, Some(error_norm), Some(dense_coeffs)) + +- [ ] Implement `interpolate()` method + - [ ] Calculate θ = (t - t_start) / (t_end - t_start) + - [ ] Evaluate 3rd order interpolation polynomial + - [ ] Return interpolated state + +- [ ] Implement constants + - [ ] `ORDER = 3` + - [ ] `STAGES = 4` + - [ ] `ADAPTIVE = true` + - [ ] `DENSE = true` + +### Integration with Problem + +- [ ] Export BS3 in prelude +- [ ] Add to `integrator/mod.rs` module exports + +### Testing + +- [ ] **Convergence test**: Linear problem (y' = λy) + - [ ] Run with decreasing tolerances + - [ ] Verify 3rd order convergence rate + - [ ] Compare to analytical solution + +- [ ] **Accuracy test**: Exponential decay + - [ ] y' = -y, y(0) = 1 + - [ ] Verify error < tolerance at t=5 + - [ ] Check intermediate points via interpolation + +- [ ] **FSAL test**: Verify function evaluation count + - [ ] Count evaluations for multi-step integration + - [ ] Should be ~4n for n accepted steps (plus rejections) + +- [ ] **Dense output test**: + - [ ] Interpolate at multiple points + - [ ] Verify 3rd order accuracy of interpolation + - [ ] Compare to fine-step reference solution + +- [ ] **Comparison test**: Run same problem with DP5 and BS3 + - [ ] Verify both achieve required tolerance + - [ ] BS3 should use fewer steps at moderate tolerances + +### Benchmarking + +- [ ] Add benchmark in `benches/` + - [ ] Simple 1D problem + - [ ] 6D orbital mechanics problem + - [ ] Compare to DP5 performance + +### Documentation + +- [ ] Add docstring to BS3 struct + - [ ] Explain when to use BS3 vs DP5 + - [ ] Note FSAL property + - [ ] Reference original paper + +- [ ] Add usage example + - [ ] Show tolerance selection + - [ ] Demonstrate interpolation + +## Testing Requirements + +### Convergence Test Details + +Standard test problem: y' = -5y, y(0) = 1, exact solution: y(t) = e^(-5t) + +Run from t=0 to t=1 with tolerances: [1e-3, 1e-4, 1e-5, 1e-6, 1e-7] + +Expected: Error ∝ tolerance^3 (3rd order convergence) + +### Stiffness Note + +BS3 is an explicit method and will struggle with stiff problems. Include a test that demonstrates this limitation (e.g., Van der Pol oscillator with large μ should require many steps). + +## References + +1. **Original Paper**: + - Bogacki, P. and Shampine, L.F. (1989), "A 3(2) pair of Runge-Kutta formulas", + Applied Mathematics Letters, Vol. 2, No. 4, pp. 321-325 + - DOI: 10.1016/0893-9659(89)90079-7 + +2. **Dense Output**: + - Same paper, Section 3 + +3. **Julia Implementation**: + - `OrdinaryDiffEq.jl/lib/OrdinaryDiffEqLowOrderRK/src/low_order_rk_perform_step.jl` + - Look for `perform_step!` for `BS3` cache + +4. **Textbook Reference**: + - Hairer, Nørsett, Wanner (2008), "Solving Ordinary Differential Equations I: Nonstiff Problems" + - Chapter II.4 on embedded methods + +## Complexity Estimate + +**Effort**: Small (2-4 hours) +- Straightforward explicit RK implementation +- Similar structure to existing DP5 +- Main work is getting tableau coefficients correct and testing + +**Risk**: Low +- Well-understood algorithm +- No new infrastructure needed +- Easy to validate against reference solutions + +## Success Criteria + +- [ ] Passes convergence test with 3rd order rate +- [ ] Passes all accuracy tests within specified tolerances +- [ ] FSAL optimization verified via function evaluation count +- [ ] Dense output achieves 3rd order interpolation accuracy +- [ ] Performance comparable to Julia implementation for similar problems +- [ ] Documentation complete with examples diff --git a/roadmap/features/02-vern7-method.md b/roadmap/features/02-vern7-method.md new file mode 100644 index 0000000..43ec278 --- /dev/null +++ b/roadmap/features/02-vern7-method.md @@ -0,0 +1,243 @@ +# Feature: Vern7 (Verner 7th Order) Method + +## Overview + +Verner's 7th order method is a high-efficiency explicit Runge-Kutta method designed by Jim Verner. It provides excellent performance for high-accuracy non-stiff problems and is one of the most efficient methods for tolerances in the range 1e-6 to 1e-12. + +**Key Characteristics:** +- Order: 7(6) - 7th order solution with 6th order error estimate +- Stages: 9 +- FSAL: Yes +- Adaptive: Yes +- Dense output: 7th order continuous extension +- Optimized for minimal error coefficients + +## Why This Feature Matters + +- **High accuracy**: Essential for tight tolerance requirements (1e-8 to 1e-12) +- **Efficiency**: More efficient than repeatedly refining lower-order methods +- **Astronomical/orbital mechanics**: Common accuracy requirement +- **Auto-switching foundation**: Needed for intelligent algorithm selection (pairs with Tsit5 for tolerance-based switching) + +## Dependencies + +- None (can be implemented with current infrastructure) + +## Implementation Approach + +### Butcher Tableau + +Vern7 has a 9-stage explicit RK tableau. The full coefficients are extensive (45 A-matrix entries). + +Key properties: +- c values: [0, 0.05, 0.1, 0.25, 0.5, 0.75, 1, 1, 1] +- FSAL: k9 = k1 for next step +- Optimized for small error coefficients + +### Dense Output + +7th order Hermite interpolation using all 9 stage values. + +Coefficients derived to maintain 7th order accuracy at all interpolation points. + +### Error Estimation + +``` +err = ||u₇ - u₆|| / (atol + max(|u_n|, |u_{n+1}|) * rtol) +``` + +Where the embedded 6th order method shares most stages with the 7th order method. + +## Implementation Tasks + +### Core Algorithm + +- [ ] Define `Vern7` struct implementing `Integrator` trait + - [ ] Add tableau constants as static arrays + - [ ] A matrix (lower triangular, 9x9, only 45 non-zero entries) + - [ ] b vector (9 elements) for 7th order solution + - [ ] b* vector (9 elements) for 6th order embedded solution + - [ ] c vector (9 elements) for stage times + - [ ] Add tolerance fields (a_tol, r_tol) + - [ ] Add builder methods + - [ ] Add optional `lazy` flag for lazy interpolation (future enhancement) + +- [ ] Implement `step()` method + - [ ] Pre-allocate k array: `Vec>` with capacity 9 + - [ ] Compute k1 = f(t, y) + - [ ] Loop through stages 2-9: + - [ ] Compute stage value using appropriate A-matrix entries + - [ ] Evaluate ki = f(t + c[i]*h, y + h*sum(A[i,j]*kj)) + - [ ] Compute 7th order solution using b weights + - [ ] Compute error using (b - b*) weights + - [ ] Store all k values for dense output + - [ ] Return (y_next, Some(error_norm), Some(k_stages)) + +- [ ] Implement `interpolate()` method + - [ ] Calculate θ = (t - t_start) / (t_end - t_start) + - [ ] Use 7th order interpolation polynomial with all 9 k values + - [ ] Return interpolated state + +- [ ] Implement constants + - [ ] `ORDER = 7` + - [ ] `STAGES = 9` + - [ ] `ADAPTIVE = true` + - [ ] `DENSE = true` + +### Tableau Coefficients + +The full Vern7 tableau is complex. Options: + +1. **Extract from Julia source**: + - File: `OrdinaryDiffEq.jl/lib/OrdinaryDiffEqVerner/src/verner_tableaus.jl` + - Look for `Vern7ConstantCache` or similar + +2. **Use Verner's original coefficients**: + - Available in Verner's published papers + - Verify rational arithmetic for exact representation + +- [ ] Transcribe A matrix coefficients + - [ ] Use Rust rational literals or high-precision floats + - [ ] Add comments indicating matrix structure + +- [ ] Transcribe b and b* vectors + +- [ ] Transcribe c vector + +- [ ] Transcribe dense output coefficients (binterp) + +- [ ] Add test to verify tableau satisfies order conditions + +### Integration with Problem + +- [ ] Export Vern7 in prelude +- [ ] Add to `integrator/mod.rs` module exports + +### Testing + +- [ ] **Convergence test**: Verify 7th order convergence + - [ ] Use y' = -y with known solution + - [ ] Run with tolerances [1e-8, 1e-9, 1e-10, 1e-11, 1e-12] + - [ ] Plot log(error) vs log(tolerance) + - [ ] Verify slope ≈ 7 + +- [ ] **High accuracy test**: Orbital mechanics + - [ ] Two-body problem with known period + - [ ] Integrate for 100 orbits + - [ ] Verify position error < 1e-10 with rtol=1e-12 + +- [ ] **FSAL verification**: + - [ ] Count function evaluations + - [ ] Should be ~9n for n accepted steps (plus rejections) + - [ ] With FSAL optimization active + +- [ ] **Dense output accuracy**: + - [ ] Verify 7th order interpolation between steps + - [ ] Interpolate at 100 points between saved states + - [ ] Error should scale with h^7 + +- [ ] **Comparison with DP5**: + - [ ] Same problem, tight tolerance (1e-10) + - [ ] Vern7 should take significantly fewer steps + - [ ] Both should achieve accuracy, Vern7 should be faster + +- [ ] **Comparison with Tsit5**: + - [ ] Vern7 should be better at tight tolerances + - [ ] Tsit5 may be competitive at moderate tolerances + +### Benchmarking + +- [ ] Add to benchmark suite + - [ ] 3D Kepler problem (orbital mechanics) + - [ ] Pleiades problem (N-body) + - [ ] Compare wall-clock time vs DP5, Tsit5 at various tolerances + +- [ ] Memory usage profiling + - [ ] Verify efficient storage of 9 k-stages + - [ ] Check for unnecessary allocations + +### Documentation + +- [ ] Comprehensive docstring + - [ ] When to use Vern7 (high accuracy, tight tolerances) + - [ ] Performance characteristics + - [ ] Comparison to other methods + - [ ] Note: not suitable for stiff problems + +- [ ] Usage example + - [ ] High-precision orbital mechanics + - [ ] Show tolerance selection guidance + +- [ ] Add to README comparison table + +## Testing Requirements + +### Standard Test: Pleiades Problem + +The Pleiades problem (7-body gravitational system) is a standard benchmark: + +```rust +// 14 equations (7 bodies × 2D positions and velocities) +// Known to require high accuracy +// Non-stiff but requires many function evaluations with low-order methods +``` + +Run from t=0 to t=3 with rtol=1e-10, atol=1e-12 + +Expected: Vern7 should complete in <2000 steps while DP5 might need >10000 steps + +### Energy Conservation Test + +For Hamiltonian systems, verify energy drift is minimal: +- Simple pendulum or harmonic oscillator +- Integrate for long time (1000 periods) +- Measure energy drift at rtol=1e-10 +- Should be < 1e-9 relative error + +## References + +1. **Original Paper**: + - Verner, J.H. (1978), "Explicit Runge-Kutta Methods with Estimates of the Local Truncation Error" + - SIAM Journal on Numerical Analysis, Vol. 15, No. 4, pp. 772-790 + +2. **Coefficients**: + - Verner's website: https://www.sfu.ca/~jverner/ + - Or extract from Julia implementation + +3. **Julia Implementation**: + - `OrdinaryDiffEq.jl/lib/OrdinaryDiffEqVerner/src/` + - Files: `verner_tableaus.jl`, `verner_perform_step.jl`, `verner_caches.jl` + +4. **Comparison Studies**: + - Hairer, Nørsett, Wanner (2008), "Solving ODEs I", Section II.5 + - Performance comparisons with other high-order methods + +## Complexity Estimate + +**Effort**: Medium (6-10 hours) +- Tableau transcription is tedious but straightforward +- More stages than previous methods means more careful indexing +- Dense output coefficients are complex +- Extensive testing needed for verification + +**Risk**: Medium +- Getting tableau coefficients exactly right is crucial +- Numerical precision matters more at 7th order +- Need to verify against trusted reference + +## Success Criteria + +- [ ] Passes 7th order convergence test +- [ ] Pleiades problem completes with expected step count +- [ ] Energy conservation test shows minimal drift +- [ ] FSAL optimization verified +- [ ] Dense output achieves 7th order accuracy +- [ ] Outperforms DP5 at tight tolerances in benchmarks +- [ ] Documentation explains when to use Vern7 +- [ ] All tests pass with rtol down to 1e-14 + +## Future Enhancements + +- [ ] Lazy interpolation option (compute dense output only when needed) +- [ ] Vern6, Vern8, Vern9 for complete family +- [ ] Optimized implementation for small systems (compile-time specialization) diff --git a/roadmap/features/03-rosenbrock23.md b/roadmap/features/03-rosenbrock23.md new file mode 100644 index 0000000..558fede --- /dev/null +++ b/roadmap/features/03-rosenbrock23.md @@ -0,0 +1,324 @@ +# Feature: Rosenbrock23 Method + +## 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 +- 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) +- Stiff-aware interpolation: 2nd order dense output + +## Why This Feature Matters + +- **Stiff problems**: Many real-world ODEs are stiff (chemistry, circuit simulation, etc.) +- **Completes basic toolkit**: With DP5/Tsit5 for non-stiff + Rosenbrock23 for stiff, can handle most problems +- **Foundation for auto-switching**: Enables automatic stiffness detection and algorithm selection +- **Widely used**: MATLAB's ode23s is based on this method + +## Dependencies + +### Required Infrastructure + +- **Linear solver** (see feature #18) + - LU factorization for dense systems + - Generic `LinearSolver` trait + +- **Jacobian computation** (see feature #19) + - Finite difference approximation + - User-provided analytical Jacobian (optional) + - Auto-diff integration (future) + +### Recommended to Implement First + +1. Linear solver infrastructure +2. Jacobian computation +3. Then Rosenbrock23 + +## Implementation Approach + +### Mathematical Background + +Rosenbrock methods solve: +``` +(I - γh*J) * k_i = h*f(y_n + Σa_ij*k_j) + h*J*Σc_ij*k_j +``` + +Where: +- J is the Jacobian ∂f/∂y +- γ is a method-specific constant +- Stages k_i are computed by solving linear systems + +For Rosenbrock23: +- γ = 1/(2 + √2) ≈ 0.2928932188 +- 3 stages requiring 3 linear solves per step +- W-method: J can be approximate or outdated + +### Algorithm Structure + +```rust +struct Rosenbrock23 { + a_tol: SVector, + r_tol: f64, + // Tableau coefficients (as constants) + // Linear solver (to be injected or created) +} +``` + +### Step Procedure + +1. Compute or reuse Jacobian J = ∂f/∂y(t_n, y_n) +2. Form W = I - γh*J +3. Factor W (LU decomposition) +4. For each stage i=1,2,3: + - Compute RHS based on previous stages + - Solve W*k_i = RHS +5. Compute solution: y_{n+1} = y_n + b1*k1 + b2*k2 + b3*k3 +6. Compute error: err = e1*k1 + e2*k2 + e3*k3 +7. Store dense output coefficients + +### Tableau Coefficients + +From Shampine & Reichelt (1997) - MATLAB's ode23s: + +``` +γ = 1/(2 + √2) + +Stage matrix for ki calculations: +a21 = 1.0 +a31 = 1.0 +a32 = 0.0 + +c21 = -1.0707963267948966 +c31 = -0.31381995116890154 +c32 = 1.3846153846153846 + +Solution weights: +b1 = 0.5 +b2 = 0.5 +b3 = 0.0 + +Error estimate: +d1 = -0.17094382871185335 +d2 = 0.17094382871185335 +d3 = 0.0 +``` + +### Jacobian Management + +Key decisions: +- **When to update J**: Error signal, step rejection, every N steps +- **Finite difference formula**: Forward or central differences +- **Sparsity**: Dense for now, sparse in future +- **Storage**: Cache J and LU factorization + +### Linear Solver Integration + +```rust +trait LinearSolver { + fn factor(&mut self, matrix: &SMatrix) -> Result<(), Error>; + fn solve(&self, rhs: &SVector) -> Result, Error>; +} + +struct DenseLU { + lu: SMatrix, + pivots: [usize; D], +} +``` + +## Implementation Tasks + +### Infrastructure (Prerequisites) + +- [ ] **Linear solver trait and implementation** + - [ ] Define `LinearSolver` trait + - [ ] Implement dense LU factorization + - [ ] Add solve method + - [ ] Tests for random matrices + +- [ ] **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 + +### Core Algorithm + +- [ ] Define `Rosenbrock23` struct + - [ ] Tableau constants + - [ ] Tolerance fields + - [ ] Jacobian update strategy fields + - [ ] Linear solver instance + +- [ ] 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) + +- [ ] Implement `interpolate()` method + - [ ] 2nd order stiff-aware interpolation + - [ ] Uses k1, k2, k3 + +- [ ] Jacobian update strategy + - [ ] Update on first step + - [ ] Update on step rejection + - [ ] Update if error test suggests (heuristic) + - [ ] Reuse otherwise + +- [ ] Implement constants + - [ ] `ORDER = 3` + - [ ] `STAGES = 3` + - [ ] `ADAPTIVE = true` + - [ ] `DENSE = true` + +### Integration + +- [ ] Add to prelude +- [ ] Module exports +- [ ] Builder pattern for configuration + +### Testing + +- [ ] **Stiff test: Van der Pol oscillator** + - [ ] μ = 1000 (very stiff) + - [ ] Explicit methods would need 100000+ steps + - [ ] Rosenbrock23 should handle in <1000 steps + - [ ] Verify solution accuracy + +- [ ] **Stiff test: Robertson problem** + - [ ] Classic stiff chemistry problem + - [ ] 3 equations, stiffness ratio ~10^11 + - [ ] Verify conservation properties + - [ ] Compare to reference solution + +- [ ] **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 + +- [ ] **Jacobian update strategy test** + - [ ] Count Jacobian evaluations + - [ ] Verify not recomputing unnecessarily + - [ ] Verify updates when needed + +- [ ] **Comparison test** + - [ ] 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 + +### Documentation + +- [ ] Docstring explaining Rosenbrock methods +- [ ] When to use vs explicit methods +- [ ] Stiffness indicators +- [ ] Example with stiff problem +- [ ] Notes on Jacobian strategy + +## Testing Requirements + +### Van der Pol Oscillator + +```rust +// y1' = y2 +// y2' = μ(1 - y1²)y2 - y1 +// Initial: y1(0) = 2, y2(0) = 0 +// μ = 1000 (very stiff) +``` + +Integrate from t=0 to t=2000. + +Expected behavior: +- Explicit method: >100,000 steps or fails +- Rosenbrock23: ~500-2000 steps depending on tolerance +- Should track limit cycle accurately + +### Robertson Problem + +```rust +// y1' = -0.04*y1 + 1e4*y2*y3 +// y2' = 0.04*y1 - 1e4*y2*y3 - 3e7*y2² +// y3' = 3e7*y2² +// Conservation: y1 + y2 + y3 = 1 +``` + +Integrate from t=0 to t=1e11 (log scale output) + +Verify: +- Conservation law maintained +- Correct steady-state behavior +- Handles extreme stiffness ratio + +## References + +1. **Original Method**: + - Shampine, L.F. and Reichelt, M.W. (1997) + - "The MATLAB ODE Suite", SIAM J. Sci. Computing, 18(1), 1-22 + - DOI: 10.1137/S1064827594276424 + +2. **Rosenbrock Methods Theory**: + - Hairer, E. and Wanner, G. (1996) + - "Solving Ordinary Differential Equations II: Stiff and DAE Problems" + - Chapter IV.7 + +3. **Julia Implementation**: + - `OrdinaryDiffEq.jl/lib/OrdinaryDiffEqRosenbrock/` + - Files: `rosenbrock_perform_step.jl`, `rosenbrock_caches.jl` + +4. **W-methods**: + - Steihaug, T. and Wolfbrandt, A. (1979) + - "An attempt to avoid exact Jacobian and nonlinear equations in the numerical solution of stiff differential equations" + - Math. Comp., 33, 521-534 + +## Complexity Estimate + +**Effort**: Large (20-30 hours) +- Linear solver: 8-10 hours +- Jacobian computation: 4-6 hours +- Rosenbrock23 core: 6-8 hours +- Testing and debugging: 6-8 hours + +**Risk**: High +- First implicit method - new complexity +- Linear algebra integration +- Numerical stability issues possible +- Jacobian update strategy tuning needed + +## 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 + +## Future Enhancements + +- [ ] User-provided analytical Jacobians +- [ ] Sparse Jacobian support +- [ ] More sophisticated update strategies +- [ ] Rodas4, Rodas5P (higher-order Rosenbrock methods) +- [ ] Krylov linear solvers for large systems diff --git a/roadmap/features/04-pid-controller.md b/roadmap/features/04-pid-controller.md new file mode 100644 index 0000000..c12348e --- /dev/null +++ b/roadmap/features/04-pid-controller.md @@ -0,0 +1,240 @@ +# Feature: PID Controller + +## 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. + +**Key Characteristics:** +- Three-term control: proportional, integral, and derivative components +- More stable than PI for challenging problems +- Standard in production ODE solvers +- Can prevent oscillatory step size behavior + +## Why This Feature Matters + +- **Robustness**: Handles difficult problems that cause PI controller to oscillate +- **Industry standard**: Used in MATLAB, Sundials, and other production solvers +- **Minimal overhead**: Small computational cost for significant stability improvement +- **Smooth stepping**: Reduces erratic step size changes + +## Dependencies + +- None (extends current controller infrastructure) + +## Implementation Approach + +### Mathematical Formulation + +The PID controller determines the next step size based on error estimates from the current and previous steps: + +``` +h_{n+1} = h_n * (ε_n)^(-β₁) * (ε_{n-1})^(-β₂) * (ε_{n-2})^(-β₃) +``` + +Where: +- ε_i = error estimate at step i (normalized by tolerance) +- β₁ = proportional coefficient (typically ~0.3 to 0.5) +- β₂ = integral coefficient (typically ~0.04 to 0.1) +- β₃ = derivative coefficient (typically ~0.01 to 0.05) + +Standard formula (Hairer & Wanner): +``` +h_{n+1} = h_n * safety * (ε_n)^(-β₁/(k+1)) * (ε_{n-1})^(-β₂/(k+1)) * (h_n/h_{n-1})^(-β₃/(k+1)) +``` + +Where k is the order of the method. + +### Advantages Over PI + +- **PI controller**: Uses only current and previous error (2 terms) +- **PID controller**: Also uses rate of change of error (3 terms) +- **Result**: Anticipates trends, prevents overshoot + +### Implementation Design + +```rust +pub struct PIDController { + // Coefficients + pub beta1: f64, // Proportional + pub beta2: f64, // Integral + pub beta3: f64, // Derivative + + // Constraints + pub factor_min: f64, // qmax inverse + pub factor_max: f64, // qmin inverse + pub h_max: f64, + pub safety_factor: f64, + + // State (error history) + pub err_old: f64, // ε_{n-1} + pub err_older: f64, // ε_{n-2} + pub h_old: f64, // h_{n-1} + + // Next step guess + pub next_step_guess: TryStep, +} +``` + +## Implementation Tasks + +### Core Controller + +- [ ] Define `PIDController` struct + - [ ] Add beta1, beta2, beta3 coefficients + - [ ] Add constraint fields (factor_min, factor_max, h_max, safety) + - [ ] Add state fields (err_old, err_older, h_old) + - [ ] Add next_step_guess field + +- [ ] Implement `Controller` trait + - [ ] `determine_step()` method + - [ ] Handle first step (no history) + - [ ] Handle second step (partial history) + - [ ] Full PID formula for subsequent steps + - [ ] Apply safety factor and limits + - [ ] Update error history + - [ ] Return TryStep::Accepted or NotYetAccepted + +- [ ] Constructor methods + - [ ] `new()` with all parameters + - [ ] `default()` with standard coefficients + - [ ] `for_order()` - scale coefficients by method order + +- [ ] Helper methods + - [ ] `reset()` - clear history (for algorithm switching) + - [ ] Update state after accepted/rejected steps + +### Standard Coefficient Sets + +Different coefficient sets for different problem classes: + +- [ ] **Default (H312)**: + - β₁ = 1/4, β₂ = 1/4, β₃ = 0 + - Actually a PI controller with specific tuning + - Good general-purpose choice + +- [ ] **H211**: + - β₁ = 1/6, β₂ = 1/6, β₃ = 0 + - More conservative + +- [ ] **Full PID (Gustafsson)**: + - β₁ = 0.49/(k+1) + - β₂ = 0.34/(k+1) + - β₃ = 0.10/(k+1) + - True PID behavior + +### Integration + +- [ ] Export PIDController in prelude +- [ ] Update Problem to accept any Controller trait +- [ ] Examples using PID controller + +### Testing + +- [ ] **Comparison test: Smooth problem** + - [ ] Run exponential decay with PI and PID + - [ ] Both should perform similarly + - [ ] Verify PID doesn't hurt performance + +- [ ] **Oscillatory problem test** + - [ ] Problem that causes PI to oscillate step sizes + - [ ] Example: y'' + ω²y = 0 with varying ω + - [ ] PID should have smoother step size evolution + - [ ] Plot step size vs time for both + +- [ ] **Step rejection handling** + - [ ] Verify history updated correctly after rejection + - [ ] Doesn't blow up or get stuck + +- [ ] **Reset test** + - [ ] Algorithm switching scenario + - [ ] Verify reset() clears history appropriately + +- [ ] **Coefficient tuning test** + - [ ] Try different β values + - [ ] Verify stability bounds + - [ ] Document which work best for which problems + +### Benchmarking + +- [ ] Add PID option to existing benchmarks +- [ ] Compare step count and function evaluations vs PI +- [ ] Measure overhead (should be negligible) + +### Documentation + +- [ ] Docstring explaining PID control +- [ ] When to prefer PID over PI +- [ ] Coefficient selection guidance +- [ ] Example comparing PI and PID behavior + +## Testing Requirements + +### Oscillatory Test Problem + +Problem designed to expose step size oscillation: + +```rust +// Prothero-Robinson equation +// y' = λ(y - φ(t)) + φ'(t) +// where φ(t) = sin(ωt), λ << 0 (stiff), ω moderate +// +// This problem can cause step size oscillation with PI +``` + +Expected: PID should maintain more stable step sizes. + +### Step Size Stability Metric + +Track standard deviation of log(h_i/h_{i-1}) over the integration: +- PI controller: may have σ > 0.5 +- PID controller: should have σ < 0.3 + +## References + +1. **PID Controllers for ODE**: + - Gustafsson, K., Lundh, M., and Söderlind, G. (1988) + - "A PI stepsize control for the numerical solution of ordinary differential equations" + - BIT Numerical Mathematics, 28, 270-287 + +2. **Implementation Details**: + - Hairer, E., Nørsett, S.P., and Wanner, G. (1993) + - "Solving Ordinary Differential Equations I", Section II.4 + - PID controller discussion + +3. **Coefficient Selection**: + - Söderlind, G. (2002) + - "Automatic Control and Adaptive Time-Stepping" + - Numerical Algorithms, 31, 281-310 + +4. **Julia Implementation**: + - `OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl` + - Look for `PIDController` + +## Complexity Estimate + +**Effort**: Small (3-5 hours) +- Straightforward extension of PI controller +- Main work is getting coefficients right +- Testing requires careful problem selection + +**Risk**: Low +- Well-understood algorithm +- Minimal code changes +- Easy to validate + +## Success Criteria + +- [ ] Implements full PID formula correctly +- [ ] Handles first/second step bootstrap +- [ ] Shows improved stability on oscillatory test problem +- [ ] Performance similar to PI on smooth problems +- [ ] Error history management correct after rejections +- [ ] Documentation complete with usage examples +- [ ] Coefficient sets match literature values + +## Future Enhancements + +- [ ] Automatic coefficient selection based on problem characteristics +- [ ] More sophisticated controllers (H0211b, predictive) +- [ ] Limiter functions to prevent extreme changes +- [ ] Per-algorithm default coefficients diff --git a/roadmap/features/05-discrete-callbacks.md b/roadmap/features/05-discrete-callbacks.md new file mode 100644 index 0000000..768619c --- /dev/null +++ b/roadmap/features/05-discrete-callbacks.md @@ -0,0 +1,244 @@ +# Feature: Discrete Callbacks + +## Overview + +Discrete callbacks trigger at discrete events based on conditions that don't require zero-crossing detection. Unlike continuous callbacks which detect sign changes, discrete callbacks check conditions at specific points (e.g., after each step, at specific times, when certain criteria are met). + +**Key Characteristics:** +- Condition-based (not zero-crossing) +- Evaluated at discrete points (typically end of each step) +- No interpolation or root-finding needed +- Can trigger multiple times or once +- Complementary to continuous callbacks + +## Why This Feature Matters + +- **Common use cases**: Time-based events, iteration limits, convergence criteria +- **Simpler than continuous**: No root-finding overhead +- **Essential for many simulations**: Parameter updates, logging, termination conditions +- **Foundation for advanced callbacks**: Basis for SavingCallback, TerminateSteadyState, etc. + +## Dependencies + +- Existing callback infrastructure (continuous callbacks already implemented) + +## Implementation Approach + +### Callback Structure + +```rust +pub struct DiscreteCallback<'a, const D: usize, P> { + /// Condition function: returns true when callback should fire + pub condition: &'a dyn Fn(f64, SVector, &P) -> bool, + + /// Effect function: modifies ODE state + pub effect: &'a dyn Fn(&mut ODE), + + /// Fire only once, or every time condition is true + pub single_trigger: bool, + + /// Has this callback already fired? (for single_trigger) + pub has_fired: bool, +} +``` + +### Evaluation Points + +Discrete callbacks are checked: +1. After each successful step +2. Before continuous callback interpolation +3. Can also check before step (for preset times) + +### Interaction with Continuous Callbacks + +Priority order: +1. Discrete callbacks (checked first) +2. Continuous callbacks (if any triggered, may interpolate backward) + +### Key Differences from Continuous + +| Aspect | Continuous | Discrete | +|--------|-----------|----------| +| Detection | Zero-crossing with root-finding | Boolean condition | +| Timing | Exact (via interpolation) | At step boundaries | +| Cost | Higher (root-finding) | Lower (simple check) | +| Use case | Physical events | Logic-based events | + +## Implementation Tasks + +### Core Structure + +- [ ] Define `DiscreteCallback` struct + - [ ] Condition function field + - [ ] Effect function field + - [ ] `single_trigger` flag + - [ ] `has_fired` state (if single_trigger) + - [ ] Constructor + +- [ ] Convenience constructors + - [ ] `new()` - full specification + - [ ] `repeating()` - always repeat + - [ ] `single()` - fire once only + +### Integration with Problem + +- [ ] Update `Problem` to handle both callback types + - [ ] Separate storage: `Vec` and `Vec` + - [ ] Or unified `Callback` enum: + ```rust + pub enum Callback<'a, const D: usize, P> { + Continuous(ContinuousCallback<'a, D, P>), + Discrete(DiscreteCallback<'a, D, P>), + } + ``` + +- [ ] Update solver loop in `Problem::solve()` + - [ ] After each successful step: + - [ ] Check all discrete callbacks + - [ ] If condition true and (!single_trigger || !has_fired): + - [ ] Apply effect + - [ ] Mark as fired if single_trigger + - [ ] Then check continuous callbacks + +### Standard Discrete Callbacks + +Pre-built common callbacks: + +- [ ] **`stop_at_time(t_stop)`** + - [ ] Condition: `t >= t_stop` + - [ ] Effect: `stop` + - [ ] Single trigger: true + +- [ ] **`max_iterations(n)`** + - [ ] Requires iteration counter in Problem + - [ ] Condition: `iteration >= n` + - [ ] Effect: `stop` + +- [ ] **`periodic(interval, effect)`** + - [ ] Fires every `interval` time units + - [ ] Requires state to track last fire time + +### Testing + +- [ ] **Basic discrete callback test** + - [ ] Simple ODE + - [ ] Callback that stops at t=5.0 + - [ ] Verify integration stops exactly at step containing t=5.0 + +- [ ] **Single trigger test** + - [ ] Callback with single_trigger=true + - [ ] Condition that becomes true, false, true again + - [ ] Verify fires only once + +- [ ] **Multiple triggers test** + - [ ] Callback with single_trigger=false + - [ ] Condition that oscillates + - [ ] Verify fires each time condition is true + +- [ ] **Combined callbacks test** + - [ ] Both discrete and continuous callbacks + - [ ] Verify both types work together + - [ ] Discrete should fire first + +- [ ] **State modification test** + - [ ] Callback that modifies ODE parameters + - [ ] Verify effect persists + - [ ] Integration continues correctly + +### Benchmarking + +- [ ] Compare overhead vs no callbacks + - [ ] Should be minimal (just boolean check) +- [ ] Compare vs continuous callback for same logical event + - [ ] Discrete should be faster + +### Documentation + +- [ ] Docstring explaining discrete vs continuous +- [ ] When to use each type +- [ ] Examples: + - [ ] Stop at specific time + - [ ] Parameter update every N time units + - [ ] Terminate when condition met +- [ ] Integration with CallbackSet (future) + +## Testing Requirements + +### Stop at Time Test + +```rust +fn test_stop_at_time() { + let params = (); + fn derivative(_t: f64, y: Vector1, _p: &()) -> Vector1 { + Vector1::new(y[0]) + } + + let ode = ODE::new(&derivative, 0.0, 10.0, Vector1::new(1.0), ()); + let dp45 = DormandPrince45::new(); + let controller = PIController::default(); + + let stop_callback = DiscreteCallback::single( + &|t: f64, _y, _p| t >= 5.0, + &stop, + ); + + let mut problem = Problem::new(ode, dp45, controller) + .with_discrete_callback(stop_callback); + let solution = problem.solve(); + + // Should stop at first step after t=5.0 + assert!(solution.times.last().unwrap() >= &5.0); + assert!(solution.times.last().unwrap() < &5.5); // Reasonable step size +} +``` + +### Parameter Modification Test + +```rust +// Callback that changes parameter at t=5.0 +// Verify slope of solution changes at that point +``` + +## References + +1. **Julia Implementation**: + - `DiffEqCallbacks.jl/src/discrete_callbacks.jl` + - `OrdinaryDiffEq.jl` - check order of callback evaluation + +2. **Design Patterns**: + - "Event Handling in DifferentialEquations.jl" + - DifferentialEquations.jl documentation on callback types + +3. **Use Cases**: + - Sundials documentation on user-supplied functions + - MATLAB ODE event handling + +## Complexity Estimate + +**Effort**: Small (4-6 hours) +- Relatively simple addition +- Similar structure to existing continuous callbacks +- Main work is integration and testing + +**Risk**: Low +- Straightforward concept +- Minimal changes to solver core +- Easy to test + +## Success Criteria + +- [ ] DiscreteCallback struct defined and documented +- [ ] Integrated into Problem solve loop +- [ ] Single-trigger functionality works correctly +- [ ] Can combine with continuous callbacks +- [ ] All tests pass +- [ ] Performance overhead < 5% +- [ ] Documentation with examples + +## Future Enhancements + +- [ ] CallbackSet for managing multiple callbacks +- [ ] Priority/ordering for callback execution +- [ ] PresetTimeCallback (fires at specific predetermined times) +- [ ] Integration with save points (saveat) +- [ ] Callback composition and chaining diff --git a/roadmap/features/06-callback-set.md b/roadmap/features/06-callback-set.md new file mode 100644 index 0000000..a815d01 --- /dev/null +++ b/roadmap/features/06-callback-set.md @@ -0,0 +1,41 @@ +# Feature: CallbackSet + +## Overview + +CallbackSet allows composing multiple callbacks (both continuous and discrete) with controlled ordering and execution priority. Essential for complex simulations with multiple events. + +## Why This Feature Matters + +- Manage multiple callbacks cleanly +- Control execution order +- Enable/disable callbacks dynamically +- Foundation for advanced callback patterns + +## Dependencies + +- Discrete callbacks (feature #5) +- Continuous callbacks (already implemented) + +## Implementation Approach + +```rust +pub struct CallbackSet<'a, const D: usize, P> { + continuous_callbacks: Vec>, + discrete_callbacks: Vec>, + // Optional: priority/ordering information +} +``` + +## Implementation Tasks + +- [ ] Define CallbackSet struct +- [ ] Builder pattern for adding callbacks +- [ ] Execution order management +- [ ] Integration with Problem solve loop +- [ ] Testing with multiple callbacks +- [ ] Documentation + +## Complexity Estimate + +**Effort**: Small (3-4 hours) +**Risk**: Low diff --git a/roadmap/features/07-saveat-functionality.md b/roadmap/features/07-saveat-functionality.md new file mode 100644 index 0000000..94b01cc --- /dev/null +++ b/roadmap/features/07-saveat-functionality.md @@ -0,0 +1,51 @@ +# Feature: Saveat Functionality + +## Overview + +Save solution at specific timepoints + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/10-fbdf-method.md b/roadmap/features/10-fbdf-method.md new file mode 100644 index 0000000..b648a6b --- /dev/null +++ b/roadmap/features/10-fbdf-method.md @@ -0,0 +1,51 @@ +# Feature: FBDF Method + +## Overview + +Fixed-leading-coefficient BDF for very stiff problems + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/11-rodas4-rodas5p-methods.md b/roadmap/features/11-rodas4-rodas5p-methods.md new file mode 100644 index 0000000..28e62be --- /dev/null +++ b/roadmap/features/11-rodas4-rodas5p-methods.md @@ -0,0 +1,51 @@ +# Feature: Rodas4/Rodas5P Methods + +## Overview + +Higher-order Rosenbrock methods + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/12-auto-switching.md b/roadmap/features/12-auto-switching.md new file mode 100644 index 0000000..0267d97 --- /dev/null +++ b/roadmap/features/12-auto-switching.md @@ -0,0 +1,306 @@ +# Feature: Auto-Switching & Stiffness Detection + +## Overview + +Automatic algorithm switching detects when a problem transitions between stiff and non-stiff regimes and switches to the appropriate solver automatically. This is one of the most powerful features for robust, user-friendly ODE solving. + +**Key Characteristics:** +- Automatic stiffness detection via eigenvalue estimation +- Seamless switching between non-stiff and stiff solvers +- CompositeAlgorithm infrastructure +- Configurable switching criteria +- Basis for DefaultODEAlgorithm (solve without specifying algorithm) + +## Why This Feature Matters + +- **User-friendly**: User doesn't need to know if problem is stiff +- **Robustness**: Handles problems with changing character +- **Efficiency**: Uses fast explicit methods when possible, switches to implicit when needed +- **Production-ready**: Essential for general-purpose library +- **Real problems**: Many problems are "mildly stiff" or transiently stiff + +## Dependencies + +### Required + +- [ ] At least one stiff solver (Rosenbrock23 or FBDF) +- [ ] At least two non-stiff solvers (have DP5, Tsit5) +- [ ] BS3 recommended for completeness + +### Recommended + +- [ ] Vern7 for high-accuracy non-stiff +- [ ] Rodas4 or Rodas5P for high-accuracy stiff +- [ ] Multiple controllers (PI, PID) + +## Implementation Approach + +### Stiffness Detection + +**Eigenvalue Estimation**: +``` +ρ = ||δf|| / ||δy|| +``` + +Where: +- δy = y_{n+1} - y_n +- δf = f(t_{n+1}, y_{n+1}) - f(t_n, y_n) +- ρ approximates spectral radius of Jacobian + +**Stiffness ratio**: +``` +S = |ρ * h| / stability_region_size +``` + +If S > tolerance (e.g., 1.0), problem is stiff. + +### Algorithm Switching Logic + +1. **Detect stiffness** every few steps +2. **Switch condition**: Stiffness detected for N consecutive steps +3. **Switch back**: Non-stiffness detected for M consecutive steps +4. **Hysteresis**: N < M to avoid chattering + +Typical values: +- N = 3-5 (switch to stiff solver) +- M = 25-50 (switch back to non-stiff) + +### CompositeAlgorithm Structure + +```rust +pub struct CompositeAlgorithm { + pub nonstiff_alg: NonStiff, + pub stiff_alg: Stiff, + pub choice_function: AutoSwitchCache, +} + +pub struct AutoSwitchCache { + pub current_algorithm: AlgorithmChoice, + pub consecutive_stiff: usize, + pub consecutive_nonstiff: usize, + pub switch_to_stiff_threshold: usize, + pub switch_to_nonstiff_threshold: usize, + pub stiffness_tolerance: f64, +} + +pub enum AlgorithmChoice { + NonStiff, + Stiff, +} +``` + +### Implementation Challenges + +1. **State transfer**: When switching, need to transfer state cleanly +2. **Controller state**: Each algorithm may have different controller state +3. **Interpolation**: Dense output from previous algorithm +4. **First step**: Which algorithm to start with? + +## Implementation Tasks + +### Core Infrastructure + +- [ ] Define `CompositeAlgorithm` struct + - [ ] Generic over two integrator types + - [ ] Store both algorithms + - [ ] Store switching logic state + +- [ ] Define `AutoSwitchCache` + - [ ] Current algorithm choice + - [ ] Consecutive step counters + - [ ] Thresholds + - [ ] Stiffness tolerance + +- [ ] Implement switching logic + - [ ] Eigenvalue estimation function + - [ ] Stiffness detection + - [ ] Decision to switch + - [ ] Reset counters appropriately + +### Integrator Changes + +- [ ] Modify `Problem` to work with composite algorithms + - [ ] May need `IntegratorEnum` or dynamic dispatch + - [ ] Or: make Problem generic and handle in solve loop + +- [ ] State transfer mechanism + - [ ] Transfer y, t from one integrator to other + - [ ] Transfer/reset controller state + - [ ] Clear interpolation data + +- [ ] Solve loop modifications + - [ ] Check for switch every N steps + - [ ] Perform switch if needed + - [ ] Continue with new algorithm + +### Eigenvalue Estimation + +- [ ] Implement basic estimator + - [ ] Track previous f evaluation + - [ ] Compute ρ = ||δf|| / ||δy|| + - [ ] Update estimate smoothly (exponential moving average) + +- [ ] Handle edge cases + - [ ] Very small ||δy|| + - [ ] First step (no history) + - [ ] After callback event + +### Default Algorithm + +- [ ] `AutoAlgSwitch` function/constructor + - [ ] Takes tuple of non-stiff algorithms + - [ ] Takes tuple of stiff algorithms + - [ ] Returns CompositeAlgorithm + - [ ] With default switching parameters + +- [ ] `DefaultODEAlgorithm` (future) + - [ ] Analyzes problem + - [ ] Selects algorithms based on size, tolerance + - [ ] Returns configured CompositeAlgorithm + +### Testing + +- [ ] **Transiently stiff problem** + - [ ] Starts non-stiff, becomes stiff, then non-stiff again + - [ ] Example: Van der Pol with time-varying μ + - [ ] Verify switches at right times + - [ ] Verify solution accuracy throughout + +- [ ] **Always non-stiff problem** + - [ ] Should never switch to stiff solver + - [ ] Verify minimal overhead + +- [ ] **Always stiff problem** + - [ ] Should switch to stiff early + - [ ] Stay on stiff solver + +- [ ] **Chattering prevention** + - [ ] Problem near stiffness boundary + - [ ] Verify doesn't switch back and forth rapidly + - [ ] Hysteresis should prevent chattering + +- [ ] **State transfer test** + - [ ] Switch mid-integration + - [ ] Verify no discontinuity in solution + - [ ] Interpolation works across switch + +- [ ] **Comparison test** + - [ ] Run transient stiff problem three ways: + - [ ] Auto-switching + - [ ] Non-stiff only (should fail or be very slow) + - [ ] Stiff only (should work but possibly slower) + - [ ] Auto-switching should be nearly optimal + +### Benchmarking + +- [ ] ROBER problem (chemistry, transiently stiff) +- [ ] HIRES problem (atmospheric chemistry) +- [ ] Compare to manual algorithm selection +- [ ] Measure switching overhead + +### Documentation + +- [ ] Explain stiffness detection +- [ ] Document switching thresholds +- [ ] When auto-switching helps vs hurts +- [ ] Examples with different problem types +- [ ] How to configure switching parameters + +## Testing Requirements + +### Transient Stiffness Test + +Van der Pol oscillator with time-varying stiffness: + +```rust +// μ(t) = 100 for t < 20 +// μ(t) = 1 for 20 <= t < 40 +// μ(t) = 100 for t >= 40 +``` + +Expected behavior: +- Start with non-stiff (or quickly switch to stiff) +- Switch to non-stiff around t=20 +- Switch back to stiff around t=40 +- Solution remains accurate throughout + +Track: +- When switches occur +- Number of switches +- Total steps with each algorithm + +### ROBER Problem + +Robertson chemical kinetics: +``` +y1' = -0.04*y1 + 1e4*y2*y3 +y2' = 0.04*y1 - 1e4*y2*y3 - 3e7*y2² +y3' = 3e7*y2² +``` + +Very stiff initially, becomes less stiff. + +Expected: Should start with (or quickly switch to) stiff solver. + +## References + +1. **Stiffness Detection**: + - Shampine, L.F. (1977) + - "Stiffness and Non-stiff Differential Equation Solvers, II" + - Applied Numerical Mathematics + +2. **Auto-switching Algorithms**: + - Hairer & Wanner (1996), "Solving ODEs II", Section IV.3 + - Discussion of when to use stiff solvers + +3. **Julia Implementation**: + - `OrdinaryDiffEq.jl/lib/OrdinaryDiffEqDefault/src/default_alg.jl` + - `AutoAlgSwitch` and `default_autoswitch` functions + +4. **MATLAB's ode45/ode15s switching**: + - MATLAB documentation on automatic solver selection + +## Complexity Estimate + +**Effort**: Large (15-25 hours) +- Composite algorithm infrastructure: 6-8 hours +- Stiffness detection: 4-6 hours +- Switching logic and state transfer: 5-8 hours +- Testing and tuning: 4-6 hours + +**Risk**: Medium-High +- Complexity in state transfer +- Getting switching criteria right requires tuning +- Interaction with controllers needs care +- Edge cases (callbacks during switch, etc.) + +## Success Criteria + +- [ ] Handles transiently stiff problems automatically +- [ ] Switches at appropriate times +- [ ] No chattering between algorithms +- [ ] Solution accuracy maintained across switches +- [ ] Overhead < 10% on problems that don't need switching +- [ ] Performance within 20% of manual optimal selection +- [ ] Documentation complete with examples +- [ ] Robust to edge cases + +## Future Enhancements + +- [ ] More sophisticated stiffness detection + - [ ] Multiple detection methods + - [ ] Learning from past behavior + +- [ ] Multi-algorithm selection + - [ ] More than 2 algorithms (low/medium/high accuracy) + - [ ] Tolerance-based selection + +- [ ] Automatic tolerance selection + +- [ ] Problem analysis at start + - [ ] Estimate problem size effect + - [ ] Sparsity detection + - [ ] Initial algorithm recommendation + +- [ ] DefaultODEAlgorithm with full analysis + - [ ] Based on problem size, tolerance, mass matrix, etc. diff --git a/roadmap/features/13-default-algorithm-selection.md b/roadmap/features/13-default-algorithm-selection.md new file mode 100644 index 0000000..0792ca6 --- /dev/null +++ b/roadmap/features/13-default-algorithm-selection.md @@ -0,0 +1,51 @@ +# Feature: Default Algorithm Selection + +## Overview + +Smart defaults based on problem characteristics + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/14-automatic-initial-step-size.md b/roadmap/features/14-automatic-initial-step-size.md new file mode 100644 index 0000000..569c0d0 --- /dev/null +++ b/roadmap/features/14-automatic-initial-step-size.md @@ -0,0 +1,51 @@ +# Feature: Automatic Initial Step Size + +## Overview + +Algorithm to determine good initial dt + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/15-presettimecallback.md b/roadmap/features/15-presettimecallback.md new file mode 100644 index 0000000..78769c1 --- /dev/null +++ b/roadmap/features/15-presettimecallback.md @@ -0,0 +1,51 @@ +# Feature: PresetTimeCallback + +## Overview + +Callbacks at predetermined times + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/16-terminatesteadystate.md b/roadmap/features/16-terminatesteadystate.md new file mode 100644 index 0000000..5c11008 --- /dev/null +++ b/roadmap/features/16-terminatesteadystate.md @@ -0,0 +1,51 @@ +# Feature: TerminateSteadyState + +## Overview + +Auto-detect steady state + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/17-savingcallback.md b/roadmap/features/17-savingcallback.md new file mode 100644 index 0000000..17e0dc9 --- /dev/null +++ b/roadmap/features/17-savingcallback.md @@ -0,0 +1,51 @@ +# Feature: SavingCallback + +## Overview + +Custom saving logic + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/18-linear-solver-infrastructure.md b/roadmap/features/18-linear-solver-infrastructure.md new file mode 100644 index 0000000..35ecce5 --- /dev/null +++ b/roadmap/features/18-linear-solver-infrastructure.md @@ -0,0 +1,51 @@ +# Feature: Linear Solver Infrastructure + +## Overview + +Generic linear solver interface and dense LU + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/19-jacobian-computation.md b/roadmap/features/19-jacobian-computation.md new file mode 100644 index 0000000..ff31370 --- /dev/null +++ b/roadmap/features/19-jacobian-computation.md @@ -0,0 +1,51 @@ +# Feature: Jacobian Computation + +## Overview + +Finite difference and auto-diff Jacobians + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/20-low-storage-runge-kutta.md b/roadmap/features/20-low-storage-runge-kutta.md new file mode 100644 index 0000000..8fa1d4c --- /dev/null +++ b/roadmap/features/20-low-storage-runge-kutta.md @@ -0,0 +1,51 @@ +# Feature: Low-Storage Runge-Kutta + +## Overview + +2N/3N storage variants for large systems + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/21-ssp-methods.md b/roadmap/features/21-ssp-methods.md new file mode 100644 index 0000000..7d1b4cb --- /dev/null +++ b/roadmap/features/21-ssp-methods.md @@ -0,0 +1,51 @@ +# Feature: SSP Methods + +## Overview + +Strong Stability Preserving methods + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/22-symplectic-integrators.md b/roadmap/features/22-symplectic-integrators.md new file mode 100644 index 0000000..19225db --- /dev/null +++ b/roadmap/features/22-symplectic-integrators.md @@ -0,0 +1,51 @@ +# Feature: Symplectic Integrators + +## Overview + +Verlet, Leapfrog, KahanLi for Hamiltonian systems + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/23-verner-methods-suite.md b/roadmap/features/23-verner-methods-suite.md new file mode 100644 index 0000000..22119b7 --- /dev/null +++ b/roadmap/features/23-verner-methods-suite.md @@ -0,0 +1,51 @@ +# Feature: Verner Methods Suite + +## Overview + +Vern6, Vern8, Vern9 + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/24-sdirk-methods.md b/roadmap/features/24-sdirk-methods.md new file mode 100644 index 0000000..a8d4547 --- /dev/null +++ b/roadmap/features/24-sdirk-methods.md @@ -0,0 +1,51 @@ +# Feature: SDIRK Methods + +## Overview + +KenCarp3/4/5 for stiff problems + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/25-exponential-integrators.md b/roadmap/features/25-exponential-integrators.md new file mode 100644 index 0000000..fc01896 --- /dev/null +++ b/roadmap/features/25-exponential-integrators.md @@ -0,0 +1,51 @@ +# Feature: Exponential Integrators + +## Overview + +Exp4, EPIRK4 for semi-linear problems + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/26-extrapolation-methods.md b/roadmap/features/26-extrapolation-methods.md new file mode 100644 index 0000000..b708689 --- /dev/null +++ b/roadmap/features/26-extrapolation-methods.md @@ -0,0 +1,51 @@ +# Feature: Extrapolation Methods + +## Overview + +Richardson extrapolation with adaptive order + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/27-stabilized-methods.md b/roadmap/features/27-stabilized-methods.md new file mode 100644 index 0000000..faccdbb --- /dev/null +++ b/roadmap/features/27-stabilized-methods.md @@ -0,0 +1,51 @@ +# Feature: Stabilized Methods + +## Overview + +ROCK2, ROCK4, RKC for mildly stiff + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/28-i-controller.md b/roadmap/features/28-i-controller.md new file mode 100644 index 0000000..bff5f89 --- /dev/null +++ b/roadmap/features/28-i-controller.md @@ -0,0 +1,51 @@ +# Feature: I Controller + +## Overview + +Basic integral controller + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/29-predictive-controller.md b/roadmap/features/29-predictive-controller.md new file mode 100644 index 0000000..21e0c65 --- /dev/null +++ b/roadmap/features/29-predictive-controller.md @@ -0,0 +1,51 @@ +# Feature: Predictive Controller + +## Overview + +Advanced predictive controller + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/30-vectorcontinuouscallback.md b/roadmap/features/30-vectorcontinuouscallback.md new file mode 100644 index 0000000..a68ef56 --- /dev/null +++ b/roadmap/features/30-vectorcontinuouscallback.md @@ -0,0 +1,51 @@ +# Feature: VectorContinuousCallback + +## Overview + +Multiple simultaneous events + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/31-positivedomain.md b/roadmap/features/31-positivedomain.md new file mode 100644 index 0000000..a584f4e --- /dev/null +++ b/roadmap/features/31-positivedomain.md @@ -0,0 +1,51 @@ +# Feature: PositiveDomain + +## Overview + +Enforce positivity constraints + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/32-manifoldprojection.md b/roadmap/features/32-manifoldprojection.md new file mode 100644 index 0000000..f267889 --- /dev/null +++ b/roadmap/features/32-manifoldprojection.md @@ -0,0 +1,51 @@ +# Feature: ManifoldProjection + +## Overview + +Project onto constraint manifolds + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Medium + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/33-nonlinear-solver-infrastructure.md b/roadmap/features/33-nonlinear-solver-infrastructure.md new file mode 100644 index 0000000..841a3a8 --- /dev/null +++ b/roadmap/features/33-nonlinear-solver-infrastructure.md @@ -0,0 +1,51 @@ +# Feature: Nonlinear Solver Infrastructure + +## Overview + +Newton and quasi-Newton methods + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/34-krylov-linear-solvers.md b/roadmap/features/34-krylov-linear-solvers.md new file mode 100644 index 0000000..b168a3b --- /dev/null +++ b/roadmap/features/34-krylov-linear-solvers.md @@ -0,0 +1,51 @@ +# Feature: Krylov Linear Solvers + +## Overview + +GMRES, BiCGStab for large sparse systems + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/35-preconditioners.md b/roadmap/features/35-preconditioners.md new file mode 100644 index 0000000..d06b4cd --- /dev/null +++ b/roadmap/features/35-preconditioners.md @@ -0,0 +1,51 @@ +# Feature: Preconditioners + +## Overview + +ILU, Jacobi preconditioners + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Large + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/36-fsal-optimization.md b/roadmap/features/36-fsal-optimization.md new file mode 100644 index 0000000..3432dc0 --- /dev/null +++ b/roadmap/features/36-fsal-optimization.md @@ -0,0 +1,51 @@ +# Feature: FSAL Optimization + +## Overview + +First-Same-As-Last function reuse + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/37-custom-norms.md b/roadmap/features/37-custom-norms.md new file mode 100644 index 0000000..12ef776 --- /dev/null +++ b/roadmap/features/37-custom-norms.md @@ -0,0 +1,51 @@ +# Feature: Custom Norms + +## Overview + +User-definable error norms + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified) diff --git a/roadmap/features/38-step-stage-limiting.md b/roadmap/features/38-step-stage-limiting.md new file mode 100644 index 0000000..e6f410a --- /dev/null +++ b/roadmap/features/38-step-stage-limiting.md @@ -0,0 +1,51 @@ +# Feature: Step/Stage Limiting + +## Overview + +Limit state values during integration + +## Why This Feature Matters + +(To be detailed) + +## Dependencies + +(To be detailed) + +## Implementation Approach + +(To be detailed) + +## Implementation Tasks + +- [ ] Core implementation +- [ ] Integration with existing code +- [ ] Testing +- [ ] Documentation +- [ ] Benchmarking + +## Testing Requirements + +(To be detailed) + +## References + +1. Julia implementation: OrdinaryDiffEq.jl +2. (Additional references to be added) + +## Complexity Estimate + +**Effort**: Small + +**Risk**: (To be assessed) + +## Success Criteria + +- [ ] Implementation complete +- [ ] Tests pass +- [ ] Documentation written +- [ ] Performance acceptable + +## Future Enhancements + +(To be identified)