differential-equations/roadmap/features/01-bs3-method.md

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<D> 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 + ha[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