Stiff Systems

Simulating a combustion reaction involves thousands of chemical species across time scales from nanoseconds to seconds. Explicit RK4 needs about 10⁹ steps. BDF or Radau solve the same problem in hundreds: orders of magnitude less compute.

• **Chemical kinetics:** combustion, polymer synthesis-stiff systems spanning 10⁻¹² to 10³ seconds; LSODA/BDF are the standard
• **Power grid simulation:** fast transients and slow dynamics in the same system; implicit methods are mandatory
• **Stiff Neural ODEs:** some control and physics problems produce stiff dynamics; adjoint+BDF in torchdiffeq

Предварительные знания

• Numerical Methods for ODEs

What is Stiffness

A stiff ODE system contains components that change at vastly different rates. Formally: the ratio of the largest to smallest eigenvalue magnitudes of the Jacobian is large. For explicit methods, this forces an extremely small step for stability, even when accuracy does not require it.

**Stiffness ratio:** stiffness ratio = max|Re(λᵢ)| / min|Re(λᵢ)| where λᵢ are eigenvalues of ∂f/∂y. **Explicit Euler stability:** stable when h·|λ| ≤ 2 For a stiff system with λ_max = 10⁶: - Explicit Euler needs h < 2·10⁻⁶ (millions of steps) - Implicit Euler: stable for any h (A-stability) **Example:** dy/dt = −1000y + cos(t)

Signs of stiffness: solve_ivp with RK45 takes thousands of tiny steps, integration is unexpectedly slow. Fix: use method='Radau' or method='BDF' in solve_ivp.

Implicit Euler and A-Stability

Implicit Euler defines yₙ₊₁ through an equation at the **new time level** rather than from the known yₙ alone. Each step requires solving a nonlinear system (via Newton's method), but the payoff is A-stability: the method is stable for any h > 0.

**Implicit Euler:** yₙ₊₁ = yₙ + h·f(tₙ₊₁, yₙ₊₁) Each step solves: G(yₙ₊₁) = yₙ₊₁ − yₙ − h·f(tₙ₊₁, yₙ₊₁) = 0 → Newton's method: O(n²) or O(n³) per step **A-stability:** stable for all λ with Re(λ) < 0 at any h Explicit Euler stability region: |1 + hλ| ≤ 1 Implicit Euler: |1/(1 − hλ)| ≤ 1 → entire left half-plane!

BDF Methods and LSODA

BDF (Backward Differentiation Formula) is a family of implicit multistep methods for stiff ODEs. They use several previous values yₙ, yₙ₋₁, ... for higher accuracy. LSODA automatically switches between stiff (BDF) and non-stiff (Adams) strategies.

**BDF of order k:** Σⱼ₌₀ᵏ aⱼ yₙ₊₁₋ⱼ = h·β₀·f(tₙ₊₁, yₙ₊₁) BDF1 = Implicit Euler: O(h) BDF2: O(h²), BDF4: O(h⁴), BDF6: O(h⁶) **Stability cutoff:** BDF7 and higher are not zero-stable (Cryer, 1971), so order ≤ 6 in practice. **LSODA:** auto-switches between Adams (non-stiff) and BDF (stiff), detecting stiffness on the fly. `solve_ivp(f, ..., method='LSODA')` or `method='BDF'`

Key Ideas

• **Stiffness:** large max|λ|/min|λ| ratio forces explicit methods to take h → 0 for stability
• **Implicit Euler:** A-stable for any h; each step is expensive (Newton for nonlinear); O(h) accuracy
• **BDF methods:** implicit multistep; orders 1-6; BDF7 and higher fail zero-stability (Cryer 1971)
• **Practice:** method='Radau' for stiff problems; method='LSODA' when stiffness is unknown

Related Topics

Stiffness connects ODE methods to linear algebra and PDE applications:

• Numerical Methods for ODEs — Explicit methods (RK45) are the starting point; stiffness explains why they are inefficient for some problems
• Finite Difference Methods for PDEs — PDE discretization produces stiff ODE systems; the CFL condition is the stability condition for explicit methods

Вопросы для размышления

• How can stiffness in an ODE system be detected without computing the Jacobian eigenvalues explicitly?
• Radau IIA is an implicit Runge-Kutta method: why is it not bound by Dahlquist's second barrier in the way linear multistep methods are?
• What happens to Neural ODE training when the dynamics are stiff, and how does that change the adjoint method?

Связанные уроки

• la-01-vectors-intro