Implicit vs Explicit Methods: Stability, Accuracy, and When to Use Each

Choosing a time-integration method is one of the most important decisions in a scientific simulation. The method determines how the numerical solution moves from one time level to the next, how small the time step must be, how much each step costs, and whether fast physical processes are resolved or suppressed.

The common distinction between explicit and implicit methods is useful, but it does not provide a complete selection rule. An implicit method is not automatically more accurate, and an explicit method is not automatically unsuitable for serious simulations. Stability, accuracy, computational cost, stiffness, numerical damping, and the physical time scales of the problem must be considered together.

The central principle is simple: numerical stability only tells us whether errors remain controlled. It does not tell us whether the computed solution is close to the true physical solution.

What Is Time Integration?

After a PDE has been discretized in space, it often becomes a system of ordinary differential equations:

du/dt = F(u, t)

A time integrator approximates how the vector u changes over a finite step:

tⁿ → tⁿ⁺¹ = tⁿ + Δt

The exact evolution is usually unavailable, so the algorithm constructs an approximation from known values, derivative evaluations, or a system involving the unknown future state.

An explicit method calculates the new state directly from information already available. An implicit method defines the new state through an equation that must be solved.

Explicit Time Integration

Forward Euler is the simplest explicit method:

uⁿ⁺¹ = uⁿ + Δt F(uⁿ, tⁿ)

Everything on the right side is known. No linear or nonlinear system is required. This makes each step inexpensive and easy to parallelize.

Higher-order explicit Runge-Kutta methods calculate several intermediate stages. The classical fourth-order method uses four derivative evaluations:

k₁ = F(uⁿ, tⁿ)

k₂ = F(
    uⁿ + 0.5 Δt k₁,
    tⁿ + 0.5 Δt
)

k₃ = F(
    uⁿ + 0.5 Δt k₂,
    tⁿ + 0.5 Δt
)

k₄ = F(
    uⁿ + Δt k₃,
    tⁿ + Δt
)

uⁿ⁺¹ = uⁿ
      + Δt(k₁ + 2k₂ + 2k₃ + k₄) / 6

Explicit methods are attractive when each derivative evaluation is affordable and the stability limit does not force an excessive number of steps.

Implicit Time Integration

Backward Euler evaluates the derivative at the unknown future state:

uⁿ⁺¹ = uⁿ + Δt F(uⁿ⁺¹, tⁿ⁺¹)

The new value appears on both sides. A linear problem may require a matrix solve, while a nonlinear problem may require Newton iterations or another nonlinear algorithm.

Crank–Nicolson averages the derivative between the current and future states:

uⁿ⁺¹ = uⁿ
      + 0.5 Δt [
          F(uⁿ, tⁿ)
          + F(uⁿ⁺¹, tⁿ⁺¹)
        ]

Implicit methods cost more per step, but suitable schemes can remain stable for time steps that would make an explicit method diverge. This is especially valuable for stiff systems and fine spatial grids.

Stability Is Not Accuracy

A stable computation does not necessarily provide an accurate physical trajectory. The distinction can be studied with the linear test equation:

dy/dt = λy

The exact solution after one time step is:

y(t + Δt) = exp(λΔt)y(t)

A numerical method instead produces:

yⁿ⁺¹ = R(z)yⁿ

z = λΔt

The function R(z) is the amplification factor. Absolute stability requires:

|R(z)| ≤ 1

This condition prevents unbounded numerical growth for a decaying test problem. It does not guarantee that R(z) closely approximates exp(z).

An implicit method may remain bounded with a very large time step while reproducing the decay rate, phase, or transient response poorly. The NASA report Explicit, Implicit, and Hybrid Methods discusses the need to consider accuracy rather than using stability alone to justify a method.

Understanding Stability Regions

The stability region is the set of values of z = λΔt for which the amplification factor remains bounded.

Forward Euler

Forward Euler has:

R(z) = 1 + z

Its stability region satisfies:

|1 + z| ≤ 1

This forms a disk centered at −1 with radius one. Along the negative real axis, the stable interval is:

−2 ≤ z ≤ 0

Classical RK4

The classical fourth-order Runge-Kutta method has a larger but still bounded stability region. Along the negative real axis, it remains stable approximately until:

z ≈ −2.785

This is considerably larger than the Forward Euler interval, but no explicit Runge-Kutta method can include the entire left half of the complex plane.

A practical introduction to absolute stability, A-stability, and L-stability is available in the Crash Course on Numerical ODEs.

Backward Euler

Backward Euler has:

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

Its stability region contains the complete left half-plane. It is therefore A-stable.

As z becomes increasingly negative, the amplification factor approaches zero. Strongly decaying modes are rapidly suppressed. This makes Backward Euler L-stable, although it is only first-order accurate.

Crank–Nicolson

Crank–Nicolson has:

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

It is also A-stable because its stability region includes the left half-plane. However, as z → −∞:

R(z) → −1

Highly stiff modes are not driven to zero. Instead, they can alternate in sign while retaining nearly constant magnitude. Crank–Nicolson is therefore not L-stable and may produce nonphysical temporal oscillations when very large steps are applied to stiff systems.

A-Stability and L-Stability

An A-stable method is stable for every test-equation eigenvalue with a non-positive real part, regardless of the time-step size.

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

R(z) → 0 as z → −∞

This distinction is important for stiff systems. A-stability prevents explosive growth, while L-stability ensures that unresolved fast decaying modes are strongly damped.

Not every implicit method is A-stable, and not every A-stable method is L-stable. The properties belong to the individual scheme rather than to the entire implicit category.

Order of Accuracy

The order of a method determines how quickly its error decreases as the time step becomes smaller.

For a method of order p:

Local truncation error = O(Δt^(p+1))
Global error           = O(Δt^p)

Forward and Backward Euler have the same formal order even though their stability properties differ greatly. RK4 can be much more accurate than Backward Euler at the same step size when stability permits its use.

The explicit or implicit label primarily describes how a step is calculated. It does not determine the formal order.

What Is Stiffness?

A system is stiff when it contains strongly separated time scales and explicit stability requirements force steps much smaller than those needed to resolve the behavior of interest.

Consider:

dy/dt = -1000(y - cos(t)) - sin(t)

The desired solution may vary on a time scale of order one, but a rapidly decaying component has a time scale near 0.001. An explicit method may need to resolve the fast mode for stability even after that mode has become physically unimportant.

A suitable implicit method can step over the fast decay and follow the slower solution. This is the main reason implicit methods are used for stiff reaction systems, diffusion equations, electrical circuits, and tightly coupled multiphysics models.

Numerical Damping

Implicitness does not automatically imply strong damping. Damping is controlled by the method’s amplification factor.

Backward Euler strongly suppresses modes when |λΔt| is large. This can be desirable when those modes represent unresolved stiffness. It can be undesirable when they represent waves or transients that should be measured.

Crank–Nicolson introduces much less high-frequency damping. This preserves some oscillatory behavior, but it can also allow unwanted stiff numerical oscillations to remain.

The FLOW-3D discussion of implicit and explicit numerical methods illustrates how large implicit steps can distort transient behavior. The effect should not be interpreted as a universal fixed under-relaxation factor. Its magnitude depends on the integration scheme, time step, equation, and iterative solver.

Under-Relaxation Is a Separate Issue

Under-relaxation is often used inside nonlinear or coupled iterative solvers:

u(updated) =
    u(old)
    + α [
        u(computed)
        - u(old)
      ]

The parameter α is usually between zero and one. Smaller values can stabilize an iterative solve but slow its convergence and alter the apparent transient when iterations are stopped before full convergence.

Under-relaxation is not an unavoidable property of every implicit time integrator. It is an additional algorithmic choice that may appear inside the nonlinear solution process.

CFL Restrictions for Explicit Methods

For an advection equation, explicit methods commonly follow a Courant condition:

Δt ≤ C Δx / |v|

The constant C depends on the spatial method and time integrator.

For an explicitly integrated diffusion equation, the limit typically scales as:

Δt ≤ C Δx² / D

This quadratic dependence can become expensive on fine meshes. Halving the cell size may require approximately four times as many time steps for a diffusion-controlled explicit scheme.

These restrictions do not mean explicit methods are inaccurate. They define a stability range. In hyperbolic simulations, the physical need to resolve wave travel may already require a step similar to the CFL limit.

Cost per Step

Explicit methods usually require function evaluations, flux calculations, or sparse matrix-vector products. Their steps are relatively inexpensive and often scale well on parallel hardware.

Implicit methods may require:

• Matrix assembly
• Jacobian construction
• Linear system solution
• Preconditioner setup
• Nonlinear Newton iterations
• Convergence checks

An implicit method is efficient only when the larger usable step compensates for the additional cost of each solve.

The comparison should therefore use total cost at a fixed error level rather than step count alone.

Adaptive Explicit Integration

Embedded Runge-Kutta pairs estimate error without completing two entirely independent integrations.

Dormand–Prince 5(4), often called RK45, shares a set of intermediate stages to construct both a fifth-order approximation and a lower-order error estimate.

The normalized error may be evaluated as:

error_ratio =
    estimated_error
    / (
        absolute_tolerance
        + relative_tolerance
          * solution_scale
      )

If the ratio is below one, the step can be accepted. If it exceeds one, the step is rejected and repeated with a smaller Δt.

A typical update has the form:

Δt(new) =
    safety
    * Δt(old)
    * error_ratio^(-1/(p+1))

Practical implementations also limit how rapidly the step may grow or shrink.

Adaptive Implicit Integration

Implicit solvers may estimate error through embedded formulas, variable-order BDF methods, defect estimates, or step doubling.

Step doubling compares:

• One step of length Δt
• Two steps of length Δt/2

The difference estimates temporal error. This may require several implicit solves, although matrix factorizations or preconditioners can sometimes be reused when the operator remains similar.

Large production models may combine several independent restrictions. The PISM time-stepping documentation demonstrates how CFL, diffusivity, output, and model-specific limits interact in a real simulation code.

IMEX Methods

Implicit-explicit methods split the right-hand side into stiff and non-stiff components:

du/dt = Fexplicit(u) + Fimplicit(u)

The inexpensive non-stiff term is evaluated explicitly, while the stiff term is treated implicitly.

For a convection-diffusion problem:

∂u/∂t
+ v · ∇u
= D∇²u

the advection term may be explicit and the diffusion term implicit. This avoids a global nonlinear solve for the complete equation while removing the severe explicit diffusion restriction.

IMEX schemes require compatible explicit and implicit formulas. Their order and stability depend on the complete paired method, not only on each component in isolation.

Explicit, Implicit, and IMEX Comparison

Selecting a Method by Physics

Method Order Still Matters

A first-order implicit method may require a small time step for accuracy even when stability permits a large one. A fourth- or fifth-order explicit method may be far more efficient for a smooth non-stiff problem.

Conversely, a high-order explicit method cannot overcome severe stiffness if its stability region excludes the relevant eigenvalues.

Selection therefore involves two separate questions:

1. Is the stability region suitable for the system and intended step size?
2. Is the order high enough to meet the required error at an affordable cost?

Verification Through Time-Step Refinement

A time integrator should be tested by repeating the simulation with smaller steps. Compare physically relevant quantities such as:

• Peak amplitude
• Wave arrival time
• Oscillation phase
• Total mass or energy
• Reaction yield
• Interface position
• Steady-state value

If the result changes significantly after the time step is halved, the original step was not temporally converged.

Temporal testing should be separated from mesh convergence. Refining both space and time simultaneously makes it difficult to determine which error source caused the change.

Common Selection Mistakes

One common mistake is taking a very large implicit step simply because the method remains stable.

Other frequent errors include:

• Confusing local and global order
• Assuming every implicit method is A-stable
• Assuming every A-stable method strongly damps stiff modes
• Using explicit diffusion on a fine mesh without estimating its stability limit
• Using Backward Euler when phase accuracy is important
• Using Crank–Nicolson for severe stiffness without checking temporal oscillations
• Ignoring nonlinear-solver tolerances in an implicit method
• Comparing algorithms at different accuracy levels
• Reporting stability without performing time-step convergence tests
• Applying a generic explicit-versus-implicit rule to every PDE

A Practical Selection Workflow

1. Identify the important physical time scales.
2. Determine whether the semi-discrete system is stiff.
3. Estimate advection, diffusion, reaction, and wave-based restrictions.
4. Decide whether fast modes must be resolved or may be damped.
5. Select a method with a suitable stability region.
6. Select an order that can meet the accuracy target.
7. Include the cost of matrix and nonlinear solves.
8. Use adaptive stepping where appropriate.
9. Repeat the simulation with tighter tolerances or smaller steps.
10. Compare total runtime at the same measured error.

Related Guides

• Time Integration Methods for PDE Solvers — A foundational introduction to explicit, implicit, and hybrid formulations.
• Operator Splitting Methods: Strang Splitting and IMEX Schemes — Learn how coupled physical processes can use different solvers.
• When to Use FEM, FVM, or FDM — Compare the main choices for spatial discretization.
• Mesh Quality and Convergence Studies — Test whether spatial refinement changes the computed result.
• Verification vs. Validation in Scientific Simulations — Place temporal convergence inside the complete credibility process.

Further Reading

• Belytschko, T. et al. Explicit, Implicit, and Hybrid Methods. NASA Technical Report.
• George, J. D. et al. Explicit Complex Time Integrators for Stiff Problems.
• Kim, J. Crash Course on Numerical ODEs: Absolute Stability, A-Stability, and L-Stability.
• FLOW-3D CFD-101. Implicit Versus Explicit Numerical Methods.
• PISM Documentation. Understanding Adaptive Time Stepping.

Conclusion

Explicit and implicit time integrators solve different numerical problems. Explicit methods provide inexpensive steps and are effective for non-stiff equations, waves, and problems whose physical resolution already requires small time increments. Implicit methods can avoid restrictive stability limits and are often necessary for stiff reactions, diffusion-dominated models, and tightly coupled systems.

Stability does not guarantee accuracy. An implicit calculation can remain bounded while missing fast transients, introducing phase error, or using a time step too large to reproduce the physical trajectory.

The individual method matters more than its broad label. Backward Euler is strongly damping and first-order accurate. Crank–Nicolson is second-order and A-stable but does not suppress extremely stiff modes. RK4 offers high accuracy for non-stiff systems but has a bounded stability region. IMEX methods combine explicit and implicit treatments when the operators can be separated.

The right method is the one that meets the required error at the lowest credible computational cost. That decision should be demonstrated through stability analysis, adaptive error control, and time-step refinement rather than assumed from the words “explicit” or “implicit.”