Choosing the Right NDSolve Method for Stiff Differential Equations

The Problem with Stiff Equations

When solving differential equations numerically, not all problems are created equal. A stiff equation is one where the solution changes rapidly in some regions but slowly in others, causing explicit methods to require impractically small step sizes. Consider a simple chemical reaction system where concentrations change dramatically over microseconds but need to be tracked for seconds.

Wolfram's NDSolve handles this through adaptive methods, but choosing the right one matters. By default, it uses LSODA—a method that automatically switches between non‑stiff and stiff solvers—but for certain problems, explicit control yields better results.

Understanding NDSolve’s Method Options

The Method option in NDSolve lets you specify how the integration proceeds. For stiff problems, two key approaches stand out:

• LSODA: Automatically detects stiffness and switches methods (default for many ODEs)
• BDF (Backward Differentiation Formula): A multi‑step implicit method designed specifically for stiff systems

For hybrid systems with discrete events, WhenEvent integrates seamlessly with any method, allowing you to model impacts, threshold crossings, or state resets.

Worked Example: Solving a Stiff Chemical Kinetics Problem

Let’s model a simple irreversible reaction A → B with a fast initial transient:

(* Define the ODE: dA/dt = -k*A *)
k = 1000;
eqn = {a'[t] == -k*a[t], a[0] == 1};

(* Default method – may struggle *)
sol1 = NDSolveValue[eqn, a, {t, 0, 1}];

(* Explicit BDF method – better for stiff problems *)
sol2 = NDSolveValue[eqn, a, {t, 0, 1}, Method -> "BDF"];

(* Compare accuracy at t = 1 *)
{sol1[1], sol2[1]}

The BDF method typically produces more accurate results with fewer steps for this stiff problem. You can verify by checking sol["Method"] on the solution to see which algorithm was used.

Trade‑offs and Limitations

While BDF excels at stiff ODEs, it has limitations:

• Computational cost: Implicit methods require solving nonlinear equations at each step
• Higher dimensions: Large systems may need manual Jacobian specification for efficiency
• Chaotic systems: Default tolerances may miss subtle dynamics; consider AccuracyGoal and PrecisionGoal

For high‑dimensional PDEs, NDSolve may require manual discretization before applying ODE solvers to the resulting system.

Actionable Takeaway

Always test stiffness sensitivity in your problem. Start with the default, then explicitly specify Method -> "BDF" for stiff systems. Use StiffnessSwitching if you want automatic detection with fallback to BDF. For verification, compare against known analytical solutions or refine AccuracyGoal until results stabilize.