Stiff Problems and Stiff Solvers: When Explicit Methods Fail

A time-dependent simulation may look well formulated and still force an explicit solver to take thousands of extremely small steps. Increasing the order from Forward Euler to RK4 may delay the problem, but it does not always remove it. The system may be stiff.

Stiffness appears when numerical stability requires a much smaller time step than the step needed to resolve the behavior of scientific interest. Fast components may decay almost immediately, while slower components control the long-term solution. An explicit integrator must remain stable for both.

Implicit stiff solvers can step over rapidly decaying modes because their stability regions extend much farther into the left half of the complex plane. The trade-off is that every accepted step may require Jacobian evaluation, matrix factorization, iterative linear algebra, or nonlinear solution.

The correct question is therefore not whether implicit methods are always better. It is whether the computational cost of an implicit step is lower than the cost of the many small explicit steps needed for stability.

The broader guide Time Integration Methods for PDE Solvers compares explicit, implicit, and IMEX approaches. This article focuses specifically on stiffness, stability regions, Jacobians, and stiff solver families.

What Stiffness Means

There is no single definition of stiffness that applies perfectly to every nonlinear system, numerical method, tolerance, and integration interval.

A useful practical description is:

A problem is stiff when a numerical method must take steps much smaller than the scale required for accuracy because of stability restrictions.

This definition emphasizes that stiffness depends partly on the solver. A problem that is stiff for an explicit Runge–Kutta method may be manageable for an L-stable implicit method.

It also depends on the requested output. A fast chemical transient may be scientifically important in one study and irrelevant in another. If it must be resolved, a small step is an accuracy requirement rather than an unnecessary stiffness restriction.

A Simple Stiff System

Consider two decaying components:

dy₁/dt = -y₁

dy₂/dt = -1000y₂

The exact solutions are:

y₁(t) = y₁(0)e⁻ᵗ

y₂(t) = y₂(0)e⁻¹⁰⁰⁰ᵗ

The second component disappears much faster than the first. After its initial transient, the scientifically relevant solution may be controlled almost entirely by y₁.

Forward Euler applied to the second equation is stable only when:

|1 - 1000Δt| ≤ 1

This gives:

0 ≤ Δt ≤ 0.002

The solver must continue using a step near this scale even after y₂ has become negligible. The slow component could be represented accurately with much larger steps, but explicit stability is controlled by the fast mode.

Why the Stiffness Ratio Is Incomplete

For a linear system:

dy/dt = Ay

a commonly quoted stiffness ratio is:

κ =
max |Re(λᵢ)|
/
min |Re(λᵢ)|

The ratio can be informative when all relevant eigenvalues have negative real parts and the matrix behaves regularly. A large value suggests widely separated decay rates.

However, it is not a universal test. It does not include the integration interval, the requested accuracy, the initial condition, or the direction in which the solution evolves. It can also become undefined or misleading when the Jacobian contains zero eigenvalues, conservation modes, or strongly non-normal behavior.

In a nonlinear system, the Jacobian changes along the trajectory:

J(t) =
∂f
/
∂y
evaluated at y(t)

A system may therefore move between stiff and non-stiff regions during one simulation.

Jacobian-Based Indicators

Eigenvalue estimates of the local Jacobian are useful because fast stable modes are associated with eigenvalues that have large negative real parts.

A cheaper indicator sometimes uses the scaled trace:

τ(J) =
trace(J) / d
=
div(f) / d

Over an integration interval with characteristic length T, a strongly negative value of:

T trace(J) / d

can indicate rapid phase-volume contraction and possible explicit step restrictions.

This is a heuristic indicator, not a necessary and sufficient definition of stiffness. A trace averages the eigenvalue real parts and can hide one strongly stiff direction when other directions compensate for it.

Non-normal systems can also exhibit transient growth that is not predicted reliably by eigenvalues or trace alone. Logarithmic norms, direct solver diagnostics, and observed step-size restrictions may provide additional information.

The mathematical background is summarized in the overview of stiff equations. A deeper analysis of why stiffness resists a simple eigenvalue-based definition is provided in Stiffness of ODEs.

Practical Signs of Stiffness

A model may be stiff when:

• An adaptive explicit method repeatedly reduces its step despite a smooth solution.
• The explicit solver takes many steps after a fast transient has already decayed.
• A small change in step size switches the calculation from stable to explosive.
• The Jacobian contains strongly negative modes with very different scales.
• An implicit solver completes the same calculation with far fewer accepted steps.
• Chemical reactions, diffusion, relaxation, or penalty terms introduce very short time scales.

These symptoms should not be used to dismiss possible coding or modeling errors. Incorrect signs, units, boundary conditions, or Jacobians can create behavior that resembles stiffness.

The Dahlquist Test Equation

Stability is commonly studied with:

dy/dt = λy

The exact solution after one step is:

y(t + Δt) =
e^(λΔt)y(t)

A numerical method produces:

yⁿ⁺¹ =
R(z)yⁿ

z = λΔt

The function R(z) is the stability function. The method is absolutely stable when:

|R(z)| ≤ 1

The set of all such values forms the absolute stability region.

A clear introduction to the Dahlquist equation and A-stability is available in the IIT numerical analysis notes.

Why Explicit Methods Fail

Forward Euler has:

R(z) = 1 + z

Its stability region is a disk centered at −1 with radius one. Along the negative real axis:

−2 ≤ z ≤ 0

Classical RK4 has a larger region, extending to approximately:

z ≈ −2.785

It is still bounded. No explicit Runge–Kutta method can be A-stable.

When a Jacobian eigenvalue is approximately −1000, RK4 requires a step smaller than roughly:

Δt ≲ 2.785 / 1000

A higher-order explicit method improves accuracy and may enlarge the stable interval, but it cannot extend stability arbitrarily far along the negative real axis.

A-Stability

An A-stable method contains the entire left half-plane in its stability region:

Re(z) ≤ 0
⇒
|R(z)| ≤ 1

Backward Euler is A-stable because:

R(z) =
1 / (1 - z)

For modes with negative real parts, the method remains stable regardless of the step size.

This does not mean an arbitrarily large step is accurate. A-stability removes one stability restriction but does not remove truncation error.

L-Stability

An L-stable method is A-stable and also satisfies:

R(z) → 0
as
z → -∞

This property strongly damps unresolved fast modes. Backward Euler is L-stable.

Crank–Nicolson is A-stable but not L-stable. Its stability function is:

R(z) =
(1 + z/2)
/
(1 - z/2)

As z → -∞, the factor approaches −1. A very stiff component can alternate in sign instead of decaying rapidly. This may create nonphysical oscillations even though the method remains bounded.

L-stability is especially valuable when fast transients should disappear rather than remain visible in the computed solution.

The Prothero–Robinson Problem

The Prothero–Robinson test problem has the general form:

dy/dt =
λ[y - g(t)]
+ g'(t)

Its exact solution is:

y(t) =
g(t)
+
[y(0) - g(0)]e^(λt)

Choosing:

g(t) = sin(t)

gives:

dy/dt =
λ[y - sin(t)]
+ cos(t)

and:

y(t) =
sin(t)
+
y(0)e^(λt)

when g(0)=0.

For a large negative λ, the exponential transient decays rapidly, while the smooth forcing sin(t) remains. An explicit method may still be forced to resolve the rapidly decaying mode for stability.

The problem is also used to study stiff order reduction. Some methods fail to achieve their nominal convergence order when stiffness, forcing, and boundary or constraint effects interact.

Python Comparison With SciPy

The following example compares an explicit method with two stiff solvers:

import numpy as np

from scipy.integrate import solve_ivp

stiffness = -1000.0
time_interval = (0.0, 10.0)
initial_state = [1.0]

def right_hand_side(time, state):
    value = state[0]

    return [
        stiffness
        * (value - np.sin(time))
        + np.cos(time)
    ]

def jacobian(time, state):
    return np.array([
        [stiffness]
    ])

def exact_solution(time):
    return (
        np.sin(time)
        + np.exp(stiffness * time)
    )

methods = [
    ("RK45", {}),
    ("Radau", {"jac": jacobian}),
    ("BDF", {"jac": jacobian}),
    ("LSODA", {})
]

for method, options in methods:
    solution = solve_ivp(
        fun=right_hand_side,
        t_span=time_interval,
        y0=initial_state,
        method=method,
        rtol=1e-7,
        atol=1e-9,
        **options
    )

    reference = exact_solution(
        solution.t
    )

    maximum_error = np.max(
        np.abs(
            solution.y[0]
            - reference
        )
    )

    print(method)
    print("  Success:", solution.success)
    print("  Accepted points:", solution.t.size)
    print("  Function evaluations:", solution.nfev)
    print("  Jacobian evaluations:", solution.njev)
    print("  LU factorizations:", solution.nlu)
    print("  Maximum error:", maximum_error)

The exact counts depend on tolerances, software version, initial state, Jacobian information, and solver settings. The important comparison is total computational cost at the same error level.

Families of Stiff Solvers

No single solver is best for every stiff system. Important considerations include:

• System dimension
• Jacobian sparsity and bandwidth
• Cost of evaluating the right-hand side
• Availability of analytic derivatives
• Requested tolerance
• Presence of constraints or a mass matrix
• Whether the problem changes between stiff and non-stiff phases

Backward Differentiation Formulas

BDF methods are implicit multistep formulas. They approximate the derivative using several previous solution values and evaluate the governing equation at the new time.

A first-order BDF formula is Backward Euler:

(
    yⁿ⁺¹ - yⁿ
) / Δt
=
f(tⁿ⁺¹, yⁿ⁺¹)

BDF2 uses two previous values:

(
    3yⁿ⁺¹
    - 4yⁿ
    + yⁿ⁻¹
) / (2Δt)
=
f(tⁿ⁺¹, yⁿ⁺¹)

Variable-order implementations may use several BDF orders and adapt both order and step size. Orders one and two are A-stable. Higher-order BDF formulas sacrifice full A-stability but remain useful for many stiff systems.

BDF methods are often effective for large systems when the Jacobian is sparse and a suitable sparse linear solver or preconditioner is available.

Examples include SciPy’s BDF, MATLAB’s ode15s, SUNDIALS CVODE BDF, and SciML methods such as QNDF and FBDF.

Rosenbrock Methods

Rosenbrock methods are linearly implicit Runge–Kutta methods. They linearize the stage equations and solve systems of the form:

(
    I - γΔtJ
)kᵢ
=
known stage terms

This avoids a complete nonlinear solve at every stage. The method still requires Jacobian information or a useful approximation.

Rosenbrock methods are often effective for small and medium stiff systems, especially at strict tolerances. The same matrix structure may be reused across several stages.

Examples in the SciML ecosystem include:

• Rosenbrock23
• Rodas4
• Rodas5P

The current SciML Rosenbrock documentation provides method-specific recommendations and implementation details.

SDIRK Methods

SDIRK means singly diagonally implicit Runge–Kutta. Every implicit stage uses the same diagonal coefficient in the Runge–Kutta matrix.

This produces stage systems with a repeated structure:

(
    I - γΔtJ
)kᵢ
=
stage residual

A matrix factorization or preconditioner may therefore be reused across stages when the Jacobian is held fixed.

SDIRK methods are single-step methods, which makes them easier to restart after events or discontinuities than multistep BDF formulas.

Examples include ESDIRK and KenCarp methods. Some provide L-stability and embedded error estimates.

TR-BDF2

TR-BDF2 is a composite one-step method. It combines a trapezoidal-rule stage with a second-order backward differentiation stage.

It should not be grouped with ordinary variable-order multistep BDF solvers. Its special construction gives second-order accuracy and strong stability properties, including L-stability for the standard parameter choice.

TR-BDF2 is often a practical general-purpose stiff method when moderate order and robust damping are more important than very high asymptotic order.

Fully Implicit Runge–Kutta Methods

FIRK means fully implicit Runge–Kutta. The stages are coupled and normally require solving a larger nonlinear system.

Important families include:

• Radau IIA methods
• Gauss–Legendre methods
• Lobatto methods

Radau IIA methods can provide high order, A-stability, and L-stability. They are attractive for high-accuracy stiff problems but may have greater per-step cost than Rosenbrock, SDIRK, or BDF methods.

SciPy’s Radau solver is based on a fifth-order Radau IIA method.

Stabilized Explicit Methods

Stabilized explicit methods use specially constructed stability polynomials to extend the stable interval along the negative real axis.

Examples include:

• RKC methods
• ROCK2
• ROCK4

For a method with s stages, the stable interval can grow approximately like s². This makes the approach useful for large diffusion-dominated systems where the Jacobian spectrum lies mainly near the negative real axis.

These methods are not A-stable and do not solve every type of stiffness. They may be unsuitable when the system has substantial imaginary eigenvalues, highly nonlinear reaction stiffness, differential-algebraic constraints, or difficult coupling.

The EPFL material on stabilized explicit methods illustrates how Chebyshev-based methods extend the stability domain.

IMEX Methods

When stiffness can be separated from non-stiff dynamics, an implicit-explicit method may be more efficient than treating the complete system implicitly:

dy/dt =
f_nonstiff(y)
+
f_stiff(y)

The non-stiff term is integrated explicitly, while the stiff term is integrated implicitly.

Examples include advection treated explicitly and diffusion treated implicitly, or transport treated explicitly and chemical reactions treated implicitly.

IMEX methods require a physically and numerically meaningful operator split. Their stability depends on the complete explicit-implicit pair.

Automatic Stiffness Switching

Some algorithms estimate stiffness while integrating and switch between non-stiff and stiff formulas.

LSODA switches between an Adams method and a BDF method. It is available through ODEPACK, SciPy, and several language interfaces.

In the SciML ecosystem, an automatic method can combine an explicit algorithm with a stiff solver, for example:

AutoTsit5(
    Rosenbrock23()
)

Automatic switching is useful when stiffness is unknown or appears only during part of the trajectory. Detection is not perfect, so difficult problems may still benefit from selecting a stiff solver directly.

The current SciML ODE solver guide recommends auto-switching methods when stiffness is uncertain.

The Jacobian Matters

For an implicit step, the nonlinear residual may be:

G(yⁿ⁺¹) =
yⁿ⁺¹
- yⁿ
- Δt f(tⁿ⁺¹, yⁿ⁺¹)

Newton’s method requires:

[
    I
    - ΔtJ
]δ
=
-G

where:

J =
∂f / ∂y

The linear solution can dominate runtime. A stiff solver will perform poorly if it treats a large sparse Jacobian as a dense matrix.

Performance can improve through:

• An analytic Jacobian
• Automatic differentiation
• A declared sparsity pattern
• Banded matrix storage
• Sparse direct factorization
• Krylov methods
• Physics-based preconditioning
• Matrix-free Jacobian-vector products

The SciML tutorial on solving large stiff equations demonstrates how Jacobian structure and linear solver selection affect performance.

Solver Selection Guide

Why Fixed Dimension Thresholds Are Misleading

It is tempting to recommend Rosenbrock methods below 1,000 equations, SDIRK methods up to 10,000, and BDF methods above that. Such thresholds are not reliable.

A dense 500-equation problem can be more expensive than a sparse system containing hundreds of thousands of equations. The decisive factors include:

• Number of nonzero Jacobian entries
• Fill-in during factorization
• Preconditioner quality
• Right-hand-side cost
• Required tolerance
• Frequency of Jacobian updates
• Availability of matrix-free products
• Parallel scalability

Solver choice should be tested through work-precision comparisons on the actual problem.

The Flame Propagation Example

Cleve Moler’s flame propagation example uses the nonlinear scalar equation:

dy/dt =
y²(1 - y)

The initial value is very small. The solution remains near zero for a long induction period, transitions rapidly, and then approaches one.

In the published MATLAB comparison, the non-stiff ode45 solver took 3,040 steps and 20,179 function evaluations. The stiff Rosenbrock-based ode23s solver took 99 steps and 412 function evaluations.

This is approximately thirty times fewer accepted steps, not one percent of the original step count.

The example and its Lambert W analytical representation are described in Stiff Differential Equations by Cleve Moler.

The result should not be interpreted as a universal thirtyfold advantage for implicit methods. It demonstrates how strongly solver efficiency can change when stability rather than accuracy controls an explicit integrator.

Comparing Solvers Fairly

A fair benchmark should compare methods at equivalent accuracy. Record:

• Global error or error in a quantity of interest
• Accepted and rejected steps
• Right-hand-side evaluations
• Jacobian evaluations
• Linear solves and factorizations
• Nonlinear iterations
• Memory consumption
• Total wall-clock time

A solver taking fewer steps may still be slower if every step requires an expensive factorization. A method with more function evaluations may win when the function is cheap and the competing linear algebra is costly.

Common Mistakes

• Assuming every diverging explicit solution proves stiffness
• Treating the stiffness ratio as a definitive diagnosis
• Using the Jacobian trace as a universal stiffness definition
• Assuming every implicit method is A-stable
• Assuming every A-stable method is L-stable
• Calling TR-BDF2 an ordinary multistep BDF method
• Expanding FIRK incorrectly as an explicit Runge–Kutta family
• Using a dense Jacobian for a large sparse system
• Failing to provide a Jacobian or sparsity pattern when the solver can use one
• Selecting a solver only from the number of equations
• Using a stabilized explicit method for stiffness outside its spectral design
• Comparing step counts without comparing error and runtime
• Taking very large implicit steps without checking temporal accuracy

A Practical Workflow

1. Verify the equations, units, signs, and initial conditions.
2. Run an adaptive explicit method as a diagnostic.
3. Inspect accepted steps, rejected steps, and minimum step size.
4. Estimate or analyze the Jacobian spectrum where practical.
5. Check whether the fast components must be resolved physically.
6. Try an L-stable stiff solver at the same tolerances.
7. Supply Jacobian information and sparsity.
8. Compare work and accuracy, not only accepted steps.
9. Repeat with tighter tolerances.
10. Document solver, tolerances, linear algebra, and failure criteria.

Related Guides

• Time Integration Methods for PDE Solvers — Compare explicit, implicit, and hybrid approaches.
• Operator Splitting, Strang Splitting, and IMEX Schemes — Separate stiff and non-stiff physical processes.
• Mesh Quality and Convergence Studies — Determine whether spatial refinement increases stiffness and changes the solution.
• Verification vs. Validation in Scientific Simulations — Place solver convergence inside the broader credibility workflow.
• Simulation & Modeling Projects — Browse related MatForge resources.

Further Reading

• Stiff Equation: Definitions, Indicators, and Historical Context
• Higham and Trefethen: Stiffness of ODEs
• IIT Notes on Absolute Stability and A-Stability
• EPFL Material on Stabilized Explicit Methods
• SciML ODE Solver Recommendations
• MathWorks: Stiff Differential Equations

Conclusion

Stiffness occurs when numerical stability forces a method to resolve fast components more closely than the desired output requires. Explicit Runge–Kutta methods have bounded stability regions, so large negative Jacobian modes can impose extremely small steps.

A-stable implicit methods remove this left-half-plane stability restriction. L-stable methods go further by strongly damping unresolved fast modes. These properties explain why Backward Euler, BDF, Rosenbrock, SDIRK, TR-BDF2, and Radau methods can outperform explicit solvers on stiff systems.

No single diagnostic or solver works for every problem. Eigenvalues, stiffness ratios, scaled divergence, and observed step restrictions all provide partial information. The structure of the Jacobian and the linear solver may matter more than the nominal number of equations.

Stabilized explicit methods are useful when stiffness is moderate and dominated by negative-real-axis eigenvalues. IMEX methods work when stiff and non-stiff operators can be separated. Auto-switching methods help when the regime changes during integration.

The final choice should be based on work-precision testing. The best stiff solver is the method that reaches the required scientific accuracy with the lowest reliable computational cost.