differential-equations/roadmap/features/12-auto-switching.md

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

pub struct CompositeAlgorithm<NonStiff, Stiff> {
    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:

// μ(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.