differential-equations/roadmap/features/03-rosenbrock23.md

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

struct Rosenbrock23<D> {
    a_tol: SVector<f64, D>,
    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 + b1k1 + b2k2 + b3*k3
6. Compute error: err = e1k1 + e2k2 + 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

trait LinearSolver<D> {
    fn factor(&mut self, matrix: &SMatrix<f64, D, D>) -> Result<(), Error>;
    fn solve(&self, rhs: &SVector<f64, D>) -> Result<SVector<f64, D>, Error>;
}

struct DenseLU<D> {
    lu: SMatrix<f64, D, D>,
    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

// 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

// 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