Reading Time: 3 minutes

Adaptive time stepping changes the time-step size during a simulation. The solver takes smaller steps when the solution changes rapidly and larger steps when the evolution becomes smooth.

A fixed-step calculation must use one value of Δt throughout the entire simulation. That value is often selected according to the most demanding event, transient, or stability restriction. The same small step is then used during long intervals where it is unnecessary.

An adaptive method estimates the error of every proposed step. It accepts the step when the normalized error is below the requested tolerance and rejects it when the error is too large. A controller then selects a new step size.

Error control is only one part of the process. Explicit methods must still respect CFL and diffusion limits. Implicit methods may need restrictions for accuracy, nonlinear convergence, wave resolution, events, or physical constraints. A reliable solver combines these requirements rather than allowing one error estimate to control everything.

What Is Adaptive Time Stepping?

After spatial discretization, a time-dependent PDE commonly becomes a system of ordinary differential equations:

du/dt = F(t, u)

A numerical method advances the solution from tⁿ to:

tⁿ⁺¹ = tⁿ + Δt

An adaptive solver follows an accept-or-reject cycle:

  1. Propose a time step.
  2. Calculate a trial solution.
  3. Estimate the local temporal error.
  4. Normalize the error using absolute and relative tolerances.
  5. Accept or reject the trial step.
  6. Select the next step size.

Rejected steps do not advance simulation time. The solver restores the previous accepted state and repeats the calculation with a smaller value of Δt.

Local Error and Global Error

Local truncation error is the error introduced by one time step when the step begins from the exact solution.

Global error is the accumulated difference between the numerical and exact solutions after many steps.

For a method of order p:

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

Global error = O(Δt^p)

Adaptive methods normally estimate local error. Keeping every local estimate below a tolerance does not provide a rigorous universal bound on global error, but it offers an effective practical control mechanism.

Absolute and Relative Tolerances

A raw error value has little meaning without the scale of each solution component. Adaptive solvers therefore calculate a normalized error.

A common component scale is:

scaleᵢ =
atolᵢ
+
rtol × max(
    |uᵢⁿ|,
    |uᵢ(trial)|
)

The normalized weighted root mean square error is:

E =
sqrt(
    1 / m
    × Σ [
        errorᵢ / scaleᵢ
      ]²
)

The step is normally accepted when:

E ≤ 1

The absolute tolerance controls components near zero. The relative tolerance controls components whose magnitudes are large compared with atol.

Different fields may require different absolute tolerances. A temperature measured in hundreds of kelvin should not necessarily use the same absolute tolerance as a concentration near 10⁻⁹.

How the Next Step Is Selected

If the estimated local error behaves like:

error ∝ Δt^(p+1)

a basic controller can use:

Δt(new) =
safety
× Δt(old)
× E^(-1/(p+1))

The safety factor is smaller than one. It reduces the chance that the next step will sit directly on the rejection boundary.

Production implementations also limit step-size changes:

growth_min
≤
Δt(new) / Δt(old)
≤
growth_max

The correct exponent depends on the order represented by the error estimator and the controller convention. It should not be copied blindly between unrelated methods.

I, PI, and PID Controllers

A simple integral controller uses only the current normalized error. This can cause the time-step history to alternate between aggressive increases and rejections.

PI and PID controllers also use errors from previous accepted steps. This produces smoother step-size changes and can reduce rejection rates.

A simplified PI controller has the form:

Δt(n+1) =
safety
× Δt(n)
× E(n)^(-k₁)
× E(n-1)^(k₂)

The coefficients depend on the method order and controller design. Sophisticated control is useful, but it cannot repair an inaccurate error estimator or an unsuitable time integrator.

Embedded Runge–Kutta Pairs

An embedded Runge–Kutta pair calculates two approximations from the same set of intermediate stages.

The two formulas have different orders:

u(high) =
uⁿ
+
Δt Σ bᵢ kᵢ

u(low) =
uⁿ
+
Δt Σ b̂ᵢ kᵢ

The difference provides an error estimate:

e =
u(high)
-
u(low)

Because both formulas share their stage evaluations, the estimator is much cheaper than completing two independent integrations.

Embedded explicit Runge–Kutta methods are best suited to non-stiff systems. If stiffness controls the step size, increasing or decreasing Δt

Dormand–Prince 5(4)

Dormand–Prince 5(4) is a seven-stage embedded pair. It produces a fifth-order solution and a fourth-order error estimate.

The pair uses local extrapolation: the fifth-order approximation becomes the accepted solution, while the lower-order formula controls the step.

It also has an FSAL property in ordinary use. The last derivative evaluation of one accepted step can become the first evaluation of the next step, reducing the number of new function evaluations after startup.

Dormand–Prince is used by SciPy’s RK45 and MATLAB’s ode45. It is a strong general-purpose choice for smooth non-stiff initial-value problems.

Tsit5

Tsit5 is the Tsitouras 5(4) embedded Runge–Kutta method. It is widely used in the SciML ecosystem as a default method for non-stiff problems at ordinary tolerances.

The method is designed to provide a favorable balance between error-estimation quality, stage cost, and practical efficiency.

It remains an explicit method. If a semi-discrete PDE is stiff because of diffusion, reactions, or strongly separated time scales, an implicit, Rosenbrock, BDF, or IMEX method may be more efficient.

Bogacki–Shampine 3(2)

Bogacki–Shampine 3(2) is a four-stage embedded pair. It provides a third-order accepted solution with a second-order error estimator.

It is implemented as RK23 in SciPy and underlies MATLAB’s ode23.

The pair can be useful when moderate accuracy is sufficient or when function evaluations are inexpensive. Its lower order may require more accepted steps than a fifth-order method when strict tolerances are requested.

Step-Doubling

Step-doubling estimates temporal error by comparing the same numerical method at two resolutions.

The procedure calculates:

  1. One full step of length Δt.
  2. Two half-steps of length Δt/2.
  3. The difference between the two results.

For a method of order p, an estimate of the error in the refined two-half-step result is:

e(refined) ≈
[
    u(half)
    -
    u(full)
]
/
(
    2^p - 1
)

The accepted solution should normally be the result of the two half-steps because it is more accurate.

Step-doubling requires three applications of the underlying method: one full step and two half-steps. An implicit implementation may therefore require three nonlinear or linear solves unless some setup work can be reused.

Residual and Defect Estimation

Residual-based methods assess how closely a reconstructed numerical trajectory satisfies the differential equation between accepted points.

A continuous reconstruction ũ(t) has a defect:

d(t) =
dũ/dt
-
F(t, ũ)

The defect can support an error estimate when it is combined with stability information or an appropriate norm.

The algebraic residual of an implicit solve is not automatically a temporal error estimate. Newton’s method can reduce the discrete equation residual almost to zero even when the time step is too large to represent the physical solution accurately.

Algebraic convergence and temporal truncation error therefore require separate controls.

Adaptive BDF Methods

Backward differentiation formulas are implicit multistep methods used for stiff systems. Variable-order implementations commonly switch among orders one through five.

A BDF solver stores information from previous accepted states. It constructs a predictor for the new state, solves the implicit equation, and estimates the local truncation error using divided differences or a Nordsieck-like history representation.

The solver can change:

  • The time-step size
  • The BDF order
  • The frequency of Jacobian updates
  • The linear solver or preconditioner setup

A smooth solution may justify increasing the order. An event, startup phase, discontinuity, or failed error test may force the solver to reduce its order and step size.

BDF2 is a backward differentiation formula. It is not Crank–Nicolson, although both methods are second order.

Error Control in COMSOL BDF

COMSOL’s BDF controller evaluates a weighted root mean square local error estimate. A step passes when the normalized error is below one.

The solver may use BDF orders from one to five. It examines solution smoothness, selects an order, attempts a step, and repeats the step with a smaller size when the local error test fails.

Its solver log includes:

  • Order for the current BDF order
  • Tfail for failed time-step error tests
  • NLfail for nonlinear solver failures

A high Tfail count suggests that the controller frequently proposed steps that were too large for the requested accuracy. A high NLfail count indicates that the algebraic nonlinear solve struggled.

These two failures should not be interpreted as the same problem.

The detailed mechanism is described in COMSOL’s Automatic Time Step and Order Selection in Time-Dependent Problems.

Algebraic Error and Temporal Error

An implicit method solves a nonlinear system at every time step. If Newton or the linear solver terminates too early, algebraic error can become comparable to the temporal error.

The time-step controller may then react to an error that comes from incomplete algebraic convergence rather than from the time discretization.

A credible solver configuration should make the algebraic tolerance sufficiently tighter than the temporal target. Making it unnecessarily strict also wastes work.

Some solvers coordinate nonlinear convergence and time-step selection. A nonlinear controller may reduce the step when Newton stability requires a smaller value than the BDF accuracy estimator proposes.

CFL Limits and Adaptive Stepping

Explicit error control cannot override an explicit stability restriction.

For advection, a typical CFL limit is:

Δt(CFL) =
C
× h
/
|v|max

For explicitly integrated diffusion:

Δt(diffusion) =
C
× h²
/
α

The constants depend on the spatial discretization, dimension, mesh geometry, and time integrator.

The usable step may be:

Δt =
min(
    Δt(error),
    Δt(CFL),
    Δt(diffusion),
    Δt(event),
    Δt(maximum)
)

An adaptive Runge–Kutta method may report a very small error and still become unstable if the CFL bound is violated.

Do Implicit Methods Need a CFL Limit?

An A-stable implicit discretization of a linear model may not have the same explicit CFL stability limit. That does not mean the time step can be arbitrarily large.

Large implicit steps may:

  • Miss waves and fast transients
  • Introduce excessive numerical damping
  • Lose monotonicity
  • Make Newton iterations difficult
  • Cross an event without resolving it
  • Reduce the accuracy of coupled nonlinear processes

These restrictions should be called accuracy, nonlinear-convergence, or physical-resolution limits rather than a universal CFL condition.

Multiple Time-Step Constraints

Large simulation codes often calculate several candidate step sizes. PISM, for example, considers model-specific restrictions connected with advection, diffusion, reporting intervals, and other physical processes.

Recording the reason that limits each accepted step is useful. A calculation dominated by an error controller requires different optimization from one limited by the fastest velocity or the finest mesh cell.

The PISM time-stepping documentation provides a practical example of multi-constraint control.

Output Times and Internal Steps

Requested output times do not always need to become integration steps. Many solvers use dense output or interpolation to calculate stored results between accepted internal time levels.

Forcing every output time to become an exact solver step can reduce efficiency when frequent output is requested.

Exact step alignment may still be appropriate when:

  • A discontinuous load changes at a known time.
  • A boundary condition switches abruptly.
  • An event changes the governing equations.
  • A control action must occur at a precise instant.
  • The method lacks sufficiently accurate dense output.

After a discontinuity, the multistep history or adaptive controller may need to restart.

Adaptive Time Stepping in FiPy

FiPy solves spatially discretized finite volume equations but does not turn a call such as:

equation.solve(
    var=phi,
    dt=dt,
    rtol=...,
    atol=...
)

into an embedded temporal error controller. Solver tolerances in FiPy primarily control the algebraic linear solve.

Current FiPy documentation recommends the external steppyngstounes package for adaptive step management. FiPy examples use controller errors derived from equation residuals, conservation checks, or other application-specific criteria.

See the official FiPy adaptive stepping documentation.

A FiPy Step-Doubling Example

The following educational example wraps FiPy’s first-order implicit transient discretization in a step-doubling controller.

It performs one full Backward Euler step and two half-steps. The two-half-step result is accepted when its normalized error is below one.

import numpy as np

from fipy import (
    Grid1D,
    CellVariable,
    TransientTerm,
    DiffusionTerm
)

# Mesh and physical model
cell_count = 100
domain_length = 1.0
cell_size = domain_length / cell_count
diffusivity = 0.01

mesh = Grid1D(
    nx=cell_count,
    dx=cell_size
)

phi = CellVariable(
    name="concentration",
    mesh=mesh,
    value=0.0
)

x = mesh.cellCenters[0]

phi.setValue(
    1.0,
    where=x < 0.25
)

equation = (
    TransientTerm()
    ==
    DiffusionTerm(
        coeff=diffusivity
    )
)

def advance_from_state(
    state,
    total_step,
    substeps
):
    trial = CellVariable(
        mesh=mesh,
        value=state,
        hasOld=True
    )

    substep_size = (
        total_step / substeps
    )

    for _ in range(substeps):
        trial.updateOld()

        equation.solve(
            var=trial,
            dt=substep_size
        )

    return np.asarray(
        trial.value
    ).copy()


# Controller settings
current_time = 0.0
final_time = 0.5

dt = 1e-3
minimum_dt = 1e-10
maximum_dt = 0.05

relative_tolerance = 1e-4
absolute_tolerance = 1e-8

safety = 0.9
minimum_factor = 0.2
maximum_factor = 3.0

# Backward Euler has order p = 1
method_order = 1

accepted_steps = 0
rejected_steps = 0

while current_time < final_time:
    dt = min(
        dt,
        final_time - current_time,
        maximum_dt
    )

    state = np.asarray(
        phi.value
    ).copy()

    # One full step
    full_step = advance_from_state(
        state=state,
        total_step=dt,
        substeps=1
    )

    # Two half-steps
    half_steps = advance_from_state(
        state=state,
        total_step=dt,
        substeps=2
    )

    # Error estimate for the refined solution
    error_vector = (
        half_steps - full_step
    ) / (
        2 ** method_order - 1
    )

    scale = (
        absolute_tolerance
        + relative_tolerance
        * np.maximum(
            np.abs(state),
            np.abs(half_steps)
        )
    )

    normalized_error = np.sqrt(
        np.mean(
            (error_vector / scale) ** 2
        )
    )

    if normalized_error <= 1.0:
        # Accept the more accurate half-step result
        phi.setValue(half_steps)
        current_time += dt
        accepted_steps += 1

        if normalized_error == 0.0:
            factor = maximum_factor
        else:
            factor = (
                safety
                * normalized_error
                ** (
                    -1.0
                    / (method_order + 1)
                )
            )

        factor = np.clip(
            factor,
            minimum_factor,
            maximum_factor
        )

        dt *= factor

    else:
        # Reject and retry from the saved state
        rejected_steps += 1

        factor = (
            safety
            * normalized_error
            ** (
                -1.0
                / (method_order + 1)
            )
        )

        factor = np.clip(
            factor,
            minimum_factor,
            1.0
        )

        dt *= factor

        if dt < minimum_dt:
            raise RuntimeError(
                "Adaptive step fell below minimum_dt"
            )

print("Accepted steps:", accepted_steps)
print("Rejected steps:", rejected_steps)
print("Final time:", current_time)

This implementation is intentionally simple. Each attempted interval requires three implicit solves, so it may be expensive. A production BDF or embedded method normally provides a more efficient estimator and a smoother PI or PID controller.

Adding an Explicit CFL Constraint

For an explicitly integrated advection term, the error-selected step can be limited separately:

maximum_speed = 2.0
cfl_number = 0.5

dt_cfl = (
    cfl_number
    * cell_size
    / maximum_speed
)

dt = min(
    dt_error,
    dt_cfl,
    maximum_dt,
    final_time - current_time
)

For non-uniform meshes, use the most restrictive local cell scale or a face-based estimate appropriate to the spatial method.

Choosing Tolerances

There is no universal pair of tolerances suitable for every PDE.

A practical process is:

  1. Choose absolute tolerances according to the physical scale of every field.
  2. Select a moderate relative tolerance.
  3. Repeat the simulation with tolerances reduced by a factor of ten.
  4. Compare the quantities of scientific interest.
  5. Continue tightening until the result changes less than the accepted temporal error.

Reducing rtol makes the requirement stricter and normally increases computational cost. Increasing rtol loosens the requirement and normally reduces accuracy.

An absolute tolerance that is too small can force excessive work when a component passes through zero. One that is too large can hide meaningful changes in small variables.

Step Rejection

A rejected explicit step wastes stage evaluations. A rejected implicit step may waste nonlinear iterations, Jacobian calculations, factorizations, and preconditioner setup.

A small number of rejected steps is normal. A high rejection rate may indicate:

  • An initial step that is too large
  • An overly aggressive controller
  • A discontinuity or event
  • Poorly scaled tolerances
  • An unreliable error estimator
  • Stiffness in an explicit method
  • Incomplete nonlinear convergence
  • An incorrect or ill-conditioned model

Controller tuning should not be used to hide model or solver failures.

When Adaptive Stepping Can Miss Important Dynamics

An error estimator measures the numerical solution generated by the chosen method. It does not know automatically which physical features matter.

A solver may take steps that are too large for:

  • Short external pulses
  • High-frequency waves
  • Threshold events
  • Contact activation
  • Rapid phase changes
  • Discontinuous source terms

Use maximum-step constraints, event handling, or forced time levels when such features must be resolved.

When Fixed Steps Are Useful

Fixed time stepping remains appropriate when:

  • A structure-preserving method requires a constant step.
  • Identical discrete trajectories are needed for debugging.
  • The solution has a nearly uniform time scale.
  • The simulation is synchronized with measured data.
  • A GPU or parallel implementation performs better with predictable workloads.
  • The time-step convergence study requires controlled refinement.

Adaptive solvers can still return results at common output times through interpolation. Fixed stepping is therefore not required solely to compare runs at specified timestamps.

Adaptive Stepping and Reproducibility

Small changes in tolerances, solver versions, hardware arithmetic, or nonlinear convergence can alter the accepted step sequence.

A reproducible report should record:

  • Time integrator and version
  • Relative and absolute tolerances
  • Initial, minimum, and maximum steps
  • Error norm
  • Controller type
  • CFL or physics-based restrictions
  • Event times
  • Accepted and rejected step counts
  • Nonlinear and linear solver tolerances

Common Mistakes

  • Using algebraic residual tolerance as a temporal error tolerance
  • Assuming FiPy changes dt through rtol and atol arguments
  • Accepting the less accurate full-step result in step-doubling
  • Forgetting the factor 2^p - 1 in Richardson estimation
  • Using a fixed controller exponent for every estimator
  • Ignoring the CFL limit of an explicit scheme
  • Calling every implicit limitation a CFL condition
  • Using Dormand–Prince for a strongly stiff diffusion problem
  • Loosening rtol when greater accuracy is required
  • Allowing a step to cross a discontinuous load without an event
  • Confusing Tfail with nonlinear solver failure
  • Comparing adaptive solvers only by accepted step count

A Practical Workflow

  1. Verify the fixed-step implementation first.
  2. Determine whether the semi-discrete system is stiff.
  3. Select an error estimator appropriate to the method.
  4. Scale every solution component with justified tolerances.
  5. Add CFL and physics-based maximum-step constraints.
  6. Handle discontinuities and events explicitly.
  7. Record accepted steps, rejections, and limiting reasons.
  8. Repeat the simulation with tighter tolerances.
  9. Compare physical quantities rather than only solver logs.
  10. Benchmark total runtime at the same measured error.

Related Guides

Further Reading

Conclusion

Adaptive time stepping improves efficiency by matching computational effort to the evolving solution. It reduces the step near rapid transients and increases it when the solution becomes smooth.

Embedded Runge–Kutta pairs estimate error from shared stages. Step-doubling compares one full step with two half-steps. Variable-order BDF methods use solution history, prediction, and divided differences to control both step size and order.

Error control does not replace stability analysis. Explicit methods must respect CFL and diffusion restrictions. Implicit methods still need steps small enough to resolve the physics and allow reliable nonlinear convergence.

FiPy does not provide automatic temporal adaptation through ordinary rtol and atol arguments to solve(). Adaptive FiPy workflows require an external step controller or custom error logic.

A credible adaptive calculation documents its tolerances, error norm, controller, rejected steps, stability constraints, event handling, and convergence under tighter settings. The objective is not to minimize the number of steps. It is to reach the required physical accuracy with the lowest reliable computational cost.