advanced33 min studyLesson 324 in Mathematics

Stiff equations and A-stability

Identify stability-limited time stepping and distinguish A-stability, L-stability, damping of fast modes, and temporal resolution.

On this page

What you’ll understand

  • Recognize stiffness as a mismatch between stability and accuracy step requirements.
  • Compare A-stable and L-stable amplification on rapid decay.
  • Explain why implicit stability does not permit arbitrary loss of physical resolution.

A simulation may contain a quickly disappearing transient together with a slowly evolving quantity that matters for the final answer. An explicit method can remain forced to follow the fast time scale even after that transient is visually negligible. This is a characteristic stiffness problem: numerical stability imposes much smaller steps than the desired slow-scale accuracy would suggest. Implicit methods can remove some of that restriction, while introducing algebraic-solver work and their own accuracy questions.

Stiffness depends on the problem and the task

A pair of decaying modes with rates minus one and minus one hundred contains a time-scale ratio of one hundred. Forward Euler must use steps below 0.02 to damp both modes, even if the slow mode changes acceptably little over a much larger interval. This is a useful model of stiffness, but a large eigenvalue ratio is not a universal definition. Initial amplitudes, forcing, requested tolerances, nonnormality, and which observables matter all affect the practical difficulty.

If the goal includes resolving the rapid initial transient, small steps are physically necessary regardless of the method. The computational advantage of a stiff integrator appears when fast stable modes can be followed accurately enough through their decay and then suppressed without destabilizing the slow evolution. A fast oscillation that remains important cannot simply be ignored because an implicit step is stable. Temporal resolution and stability are distinct constraints that may dominate at different stages of one calculation.

ys′=−ys,yf′=−100yf,RFE(−100h)=1−100hy_s'=-y_s,\qquad y_f'=-100y_f,\qquad R_{\mathrm{FE}}(-100h)=1-100h
The fast mode imposes h<0.02 for strict forward-Euler damping, although the slow mode has an order-one decay time.↔ Scroll wide equations horizontally.

A-stability covers the left half-plane

A one-step method is A-stable when its stability function has no poles in the open left half-plane and has modulus at most one throughout that half-plane. Backward Euler is A-stable: the real part of one minus z exceeds one there, so its reciprocal has modulus below one. No finite step instability occurs on a negative-real scalar mode. This is a scalar linear stability property, not a theorem that all nonlinear implicit steps have unique solutions or that every step is accurate.

The trapezoidal method is also A-stable, but its amplification approaches minus one along the negative real axis. Very rapid decays may therefore survive as alternating numerical components. L-stability adds the requirement that amplification tend to zero for increasingly stiff decay; backward Euler has that property. Strong suppression is valuable for unwanted transients, while excessive numerical damping is undesirable for physically important oscillations. Choosing a method requires deciding which behavior should be resolved and which should decay away.

RBE(z)=11−z⟶0,Rtrap(z)=1+z/21−z/2⟶−1(z→−∞)R_{\mathrm{BE}}(z)=\frac1{1-z}\longrightarrow0,\qquad R_{\mathrm{trap}}(z)=\frac{1+z/2}{1-z/2}\longrightarrow-1\quad(z\to-\infty)
Both methods are A-stable; backward Euler is L-stable, while trapezoidal integration does not strongly suppress very fast negative-real modes.↔ Scroll wide equations horizontally.

Spatial models can create very fast modes

Diffusion is a common source of stiffness after spatial discretization. Finer spatial scales decay more quickly: for a Fourier heat mode, the rate grows with the square of its wave number. A spatial mesh that resolves fine detail therefore introduces increasingly restrictive explicit time-step limits, even when the final observable is smooth. The exact Fourier-mode formula provides a clean way to see this scale separation without confusing it with a particular finite-difference implementation.

Implicit integration trades many small explicit steps for systems involving the Jacobian and the step size. Sparse factorization or preconditioned iterative solves can make that trade favorable, but their cost and convergence must be included in the comparison. Adaptive error control still matters, especially across initial layers, changes in forcing, or sharp transitions. An A-stable label is not a substitute for checking the solver's actual tolerances, nonlinear convergence, and resolution of the quantities the model is intended to predict.

u(x,t)=e−κn2tsin⁡(nx),ut=κuxxu(x,t)=e^{-\kappa n^2t}\sin(nx),\qquad u_t=\kappa u_{xx}
A periodic heat mode of wave number n has decay rate κn². This exact continuous solution illustrates scale separation, not a discretization-stability proof.↔ Scroll wide equations horizontally.

One step across two decay scales

For modes with rates −1 and −100, each initially one, compare forward Euler, backward Euler, and trapezoidal amplification at h=0.1.

  1. Forward Euler gives factors 0.9 and −9. The fast component grows ninefold and alternates sign, even though its exact factor is exp(−10)≈0.0000454.

  2. Backward Euler gives factors 1/1.1≈0.90909 and 1/11≈0.090909. Both damp, but the fast component is not accurately resolved by this single coarse step.

  3. Trapezoidal integration gives factors 0.95/1.05≈0.904762 and (1−5)/(1+5)=−2/3. It is stable, yet leaves a large alternating remnant of the fast mode.

RFE(−10)=−9,RBE(−10)=111,Rtrap(−10)=−23R_{\mathrm{FE}}(-10)=-9,\qquad R_{\mathrm{BE}}(-10)=\frac1{11},\qquad R_{\mathrm{trap}}(-10)=-\frac23
The exact fast-mode factor is exp(−10); numerical boundedness alone is a much weaker requirement.↔ Scroll wide equations horizontally.
The result

Stability distinguishes the explicit failure, while stiff damping distinguishes the two implicit methods. None of these coarse-step factors faithfully resolves the fast transient.

Explore the model

Interactive lab

Interactive lab · heat

Decay scales from spatial frequency

Pause and Reset, then set mode n=1 and diffusivity κ=1. Predict the amplitude after one Step, which advances physical time to 0.5. For a fresh comparison, Reset and set n=10 and κ=1, then Step once again. Compare the two analytic decay factors.

Decay scales from spatial frequencyThe exact periodic heat-equation mode is u(x,t) = exp(−κn²t) sin(nx), on [0,2π]. Playback advances physical time. The initial sine is shown separately, and the unnormalized L² norm is √π exp(−κn²t). This analytic instrument does not run a finite-difference solver or test discretization stability. Current readings: Current amplitude: 1; Spatial L² norm: 1.772; Decay rate · κn²: 3.6 /s. All parameters are adjustable with the labeled controls below.Periodic heat equation mode 3 at time 0: the amplitude is 1.0π/2π3π/22π−1−0.500.51x · periodic intervalTemperature perturbation uExact mode · n = 3
Current temperature modeInitial sineCurrent amplitude
Current amplitude
1
Spatial L² norm
1.772
Decay rate · κn²
3.6 /s
t = 0 s

The exact periodic heat-equation mode is u(x,t) = exp(−κn²t) sin(nx), on [0,2π]. Playback advances physical time. The initial sine is shown separately, and the unnormalized L² norm is √π exp(−κn²t). This analytic instrument does not run a finite-difference solver or test discretization stability. Parameter changes pause playback and recompute the model at the current elapsed time. Reset restores both parameters and returns time to zero. Playback starts paused. Adjust either parameter, then choose Play or Step.

Look forThe amplitudes are exp(−0.5)≈0.60653 and exp(−50), effectively zero at the readout scale. The lab evaluates an exact heat mode; it does not run an implicit solver or demonstrate the stability of a spatial discretization.

Check your understanding

Which statement distinguishes L-stability from A-stability for the methods discussed here?

Try it yourself

Evaluate backward-Euler and trapezoidal factors at z=−20. Which leaves the smaller fast-mode magnitude, and which exact factor are they approximating?

Take a moment to make a prediction or work it out on paper. Open a hint when you need one.

Give me a hint

Insert −20 into both rational stability functions before comparing absolute values.

Compare with the solution

Backward Euler gives 1/21≈0.047619. Trapezoidal integration gives (1−10)/(1+10)=−9/11≈−0.818182. The backward-Euler magnitude is much smaller; the trapezoidal remnant also changes sign. The exact factor is exp(−20)≈2.06115×10⁻⁹, so even the smaller numerical factor is a coarse approximation to the fast transient, although it damps it strongly relative to trapezoidal integration.

Take it with you

  • Stiffness often makes explicit steps stability-limited rather than slow-scale accuracy-limited.
  • A-stability and stiff-limit damping are separate properties.
  • Important fast dynamics still require resolution even when the integrator is stable.

Read further

Original sources and deeper explanations to follow at your own pace.

Optional progress, saved in this browser when available. Revisit any lesson at any time.

Lesson not marked complete.