Implicit Methods and Stiffness - Executable Engineering: Interactive Numerical Methods in Python

Up to this point, we have examined explicit solver schemes, in which the next state $y_{i+1}$ is computed directly from known information at $(t_i, y_i)$ using $f(t_i, y_i)$ .

Implicit methods evaluate the derivative at the new time step $f(t_{i+1}, y_{i+1})$ in order to determine $y_{i+1}$ ! Because $y_{i+1}$ appears on both sides of the update equation, we generally must solve a root-finding problem (or linear system) at every time step. While this adds computational expense per step, implicit methods are dramatically more resilient to numerical instability.

import numpy as np

try:
    import executable_engineering as exe
except ImportError:
    %pip install -q executable_engineering
    import executable_engineering as exe

The Backward Euler Method¶

The simplest implicit scheme is the Backward Euler (implicit Euler) method. Like Forward Euler, it approximates the integral over $[t_i, t_{i+1}]$ using a constant slope, but it selects the slope at the end of the time step:

y_{i+1} = y_i + h f(t_{i+1}, y_{i+1})

(1)

which implicitly defines $y_{i+1}$ as the root of $g(y_{i+1}) = y_{i+1} - y_i - h f(t_{i+1}, y_{i+1}) = 0$ .

differential_equations_module.py

def backward_euler(f, t0, y0, t_end, h):
    from scipy.optimize import root

    t_points = np.arange(t0, t_end + h / 2, h)
    y_points = np.zeros_like(t_points, dtype=float)
    y_points[0] = y0

    for i in range(len(t_points) - 1):
        t_next = t_points[i + 1]
        y_curr = y_points[i]
        sol = root(lambda y_val: y_val - (y_curr + h * f(t_next, y_val)), y_curr)
        y_points[i + 1] = sol.x[0]

    return {
        "t_points": t_points,
        "y_points": y_points,
        "h": h,
        "method": "Backward Euler",
    }

exe.display_source(exe.backward_euler)

f = lambda t, y: -y
res_be = exe.backward_euler(f, t0=0.0, y0=2.0, t_end=4.0, h=1.0)
print("Time nodes:     ", res_be["t_points"])
print("Backward Euler: ", np.round(res_be["y_points"], 4))

Numerical Instability¶

An algorithm is numerically unstable when small errors introduced during computation (from truncation or floating-point round-off) are amplified from step to step rather than damped. Numerical instability is a property of the algorithm and step size, not of the underlying differential equation.

Example: Integrate $\frac{\mathrm{d}y}{\mathrm{d}t} = -y$ with $y(0) = 2$ from $t = 0$ to $t = 20$ using Forward and Backward Euler and compare against the exact decay $y(t) = 2e^{-t}$ across increasing step sizes $h$ .

exe.plot_euler_stability(t_end=20.0, y0=2.0)

Analytical Stability Analysis¶

Gadzooks! Why does Forward Euler oscillate and explode when $h$ becomes large? Consider the linear test equation $\frac{\mathrm{d}y}{\mathrm{d}t} = -ay$ for a decay constant $a > 0$ . Applying Forward Euler yields:

\begin{aligned} y_{i+1} &= y_i + h \frac{\mathrm{d}y}{\mathrm{d}t}(y_i) = y_i - a h y_i \\ y_{i+1} &= (1 - a h) y_i \end{aligned}

(2)

For smooth monotonic decay, we require $1 - ah > 0 \implies h < \frac{1}{a}$ .
If $\frac{1}{a} < h < \frac{2}{a}$ , the amplification factor is negative with $|1 - ah| < 1$ , causing damped sign-alternating overshoots.
If $h > \frac{2}{a}$ , then $|1 - ah| > 1$ and the numerical solution catastrophically diverges!

Compare this with Backward Euler on the same equation:

\begin{aligned} y_{i+1} &= y_i + h \frac{\mathrm{d}y}{\mathrm{d}t}(y_{i+1}) = y_i - a h y_{i+1} \\ y_{i+1} &= \frac{y_i}{1 + a h} \end{aligned}

(3)

Because $0 < \frac{1}{1 + ah} < 1$ for all $h > 0$ , Backward Euler is monotonically decreasing and unconditionally stable! Focusing on the slope at the end of the time step naturally damps out errors from large time steps.

Stiffness¶

One might hope that adaptive explicit time steppers can simply cruise past stability limits once a fast transient dies out, but stiffness makes that impossible.

Consider the initial value problem:

\frac{\mathrm{d}y}{\mathrm{d}t} = -1000 y + 3000 - 2000 e^{-t}, \qquad y(0) = 0

(4)

Its exact analytical solution is:

y(t) = 3 - 0.998 e^{-1000 t} - 2.002 e^{-t}

(5)

which features an extremely rapid boundary-layer transient $e^{-1000t}$ (time scale $\tau_{\text{fast}} = 0.001$ ) followed by a gentle decay $e^{-t}$ (time scale $\tau_{\text{slow}} = 1$ ). Use the slider below to zoom from the initial boundary layer out to the long-time progression:

exe.plot_stiff_ode()

The Tyranny of the Fast Timescale¶

Stability of the $-1000y$ term under Forward Euler requires $h < \frac{2}{1000} = 0.002$ to avoid catastrophic blow-up. Surely that tiny step restriction only applies during the first few milliseconds ( $t < 0.01$ ) while $e^{-1000t}$ is active, right?

Wrong! Even on the interval $t \in [0.01, 0.05]$ , where $e^{-1000t} < 4.5 \times 10^{-5}$ has completely vanished from the exact curve, the $-1000y$ term in the differential equation still governs error amplification. Increase $h$ past 0.0020 in the interactive plot below to watch Explicit Euler tear itself apart while Implicit Euler remains rock solid:

exe.plot_stiff_euler_comparison(t0=0.01, t_end=0.05)

Definition of Stiffness¶

While a single universal definition of stiffness is elusive, a practical engineering definition is:

Stiff Equation: A differential equation in which the maximum allowable step size $h$ for explicit methods is dictated by numerical stability rather than by the smoothness/accuracy of the solution.

Stiffness arises constantly in engineering physics whenever coupled phenomena evolve on widely separated timescales (e.g., fast acoustic/elastic vibrations coupled with slow thermal diffusion, or chemical kinetics with fast and slow reactions). Even after the fast transients have decayed to equilibrium, explicit solvers remain hostage to the fastest timescale.

Engineers handle stiffness in two ways:

Change the physical model to eliminate the fast timescale analytically (e.g., quasi-static elasticity or incompressible flow).
Use implicit solvers (such as Backward Euler, Crank-Nicolson, or BDF methods), paying the cost of a nonlinear/linear solve per step in exchange for taking large, physically meaningful time steps.