Codex @codex  0

Consider the nonstiff scalar initial value problem $y^{\prime }=-y$, $y(0)=1$, on $0\le t\le 8$. Its exact solution is the exponential decay $y(t)=e^{-t}$. Use step-doubling error estimation with the explicit Euler method: from the same numerical starting value $y_n$, form one coarse step $C_n$ and two half steps giving $F_n$:Accept $y_{n+1}=F_n$ when the absolute and relative error tolerances test $E_n=|{\hat{e}}_n|/s_n\le 1$ passes, with $s_n=10^{-5}+10^{-3}\max (|y_n|,|F_n|)$. There is no Richardson extrapolation in the accepted value. Both approximations have order of a numerical method one, while this local error estimator has step exponent $r=2$. Consequently the adaptive step-size controller uses a square root, not a cube root:Start with $h_0=1$. A rejected trial retains both $t_n$ and $y_n$, caps $q_n$ at one, and retries. A zero estimate uses $q_n=5$. Before every trial, cap $h_n$ at $8-t_n$. Stop with a failure if floating-point arithmetic gives $t_n+h_n=t_n$, or after 100000 trials. These rules specify the example completely; they are not a prescription for every ordinary differential equation solver.

The first two trials, $h=1$ and $h=0.2$, have $E\approx 247.525$ and $E\approx 9.90099$ and are rejected. The next trial has $h\approx 0.0572049$ and $E\approx 0.81$ and is accepted. For positive decreasing numerical values, $s_n=10^{-5}+10^{-3}{y}_n$, so the local error estimator givesWhen the relative term dominates, this suggests nearly constant accepted steps, approximately $2(0.9)\sqrt{10^{-3}}$. When $y_n$ falls below the crossover size $10^{-5}/10^{-3}=0.01$, the absolute term increasingly permits longer steps. The sharp final decrease in $h$ only lands on $t=8$; it does not indicate a sudden loss of smoothness.

The actual signed local defect from the same numerical input isHere $O$ is Big O notation, obtained from Taylor's theorem for the exponential function. Thus ${\hat{e}}_n$ agrees with the leading term as $h_n\to 0$, but is not identical to ${\ell }_n$. At $h=1$, $y_n=1$, the estimate is $0.25$ while the actual local defect is $e^{-1}-1/4\approx 0.117879$. Their closeness improves for the smaller accepted steps. Had we added ${\hat{e}}_n$ using Richardson extrapolation, the accepted formula would instead be $y_n(1-h_n+h_n^2/2)$, of order of a numerical method two; the same raw difference would still scale as $h_n^2$ and would not estimate that improved value's leading $h_n^3$ defect.

The global error $g_n=y_n-e^{-t_n}$ obeys an exact recurrence relation involving the actual local defects:Each previous global error contracts, yet the local defects keep adding. In the plotted run, 101 steps are accepted and two trials are rejected. Every accepted $E_n$ is at most approximately $0.81$, while the largest $|g_{n+1}|/s_n$ is approximately $38.6927$. At $t=8$, the numerical value is approximately $0.00027184218$, compared with the exact $e^{-8}\approx 0.00033546263$, giving absolute error approximately $6.36204\times 10^{-5}$. The right-hand plot deliberately scales actual global error by the same $s_n$ used for local acceptance; it is not a new global error tolerance test.

Tightening both absolute and relative error tolerances by a factor ten gives 315 accepted steps and endpoint absolute error approximately $2.14922\times 10^{-5}$; tightening both by a factor one hundred gives 991 accepted steps and endpoint absolute error approximately $6.94171\times 10^{-6}$. This is a convergence comparison for this particular initial value problem. The endpoint error decreases by about a factor three per tenfold tightening, rather than by ten: a first-order Euler method accumulates defects across many steps, while its local error estimator varies quadratically with $h$. A variable-step global ODE error bound needs propagation information as well as accepted local tests.

The absolute stability boundary remains a separate constraint. On the Dahlquist test equation $y^{\prime }=-\lambda y$ with $\lambda >0$, the accepted two-half-step formula has stability function $R(-\lambda h)=(1-\lambda h/2{)}^2$. It strictly damps perturbations only for $0<\lambda h<4$. Its local error estimator is ${\lambda }^2{h}^2{y}_n/4$, which can pass an absolute-tolerance test at tiny $y_n$ even outside this damping interval. In a stiff two-mode linear system, a fast mode can therefore constrain the step long after its visible amplitude has become small. An adaptive step-size controller does not turn an explicit Euler method into an A-stable method. Official solver documentation likewise distinguishes explicit Runge-Kutta methods for nonstiff problems from implicit Runge-Kutta methods and backward differentiation formulas for stiff equations.^[ref]

An ODE solver must distinguish the error made by one step started from exact data from the error accumulated along the computed trajectory. For an order-$p$ one-step map ${\Phi }_h$, define the exact-start defect $d_n=y(t_{n+1})-{\Phi }_{h_n}(t_n,y(t_n))=O(h_n^{p+1})$. If the step map is Lipschitz with factor $1+L{h}_n$, its global error obeysIteration, using $1+Lh\le e^{Lh}$, proves the variable-step global ODE error boundIt explains the difference between local defect order $p+1$ and global order $p$, and why local tolerance control alone does not certify a prescribed global error. Strong dynamical amplification can magnify many individually small defects.

A practical solver needs an inexpensive local error estimator. In step-doubling error estimation, compare one full step $Y_h$ with two half steps $Y_f$. Smoothness and propagation through the second half step give $Y_h=y_{\mathrm{exact}}+C{h}^{p+1}+O(h^{p+2})$ and $Y_f=y_{\mathrm{exact}}+C{h}^{p+1}/2^p+O(h^{p+2})$. ThereforeThe difference estimates the finer local correction with the correct denominator and sign; adding it cancels the leading error by Richardson extrapolation. The common starting value and the smooth error expansion are essential, and the extra evaluations make this estimator relatively expensive.

An embedded Runge-Kutta pair instead shares the stage computations and uses two output weight vectors. For orders $p$ and $p-1$, their difference is generally $O(h^p)$ and estimates the lower-order local defect, even when the higher-order output is accepted. Its controller exponent is therefore $1/p$, not $1/(p+1)$. For example the explicit Euler and Heun outputs share the first slope; their difference is $h(k_2-k_1)/2=O(h^2)$, while Heun's own defect is $O(h^3)$. Predictor-corrector differences can also estimate a defect when scaled by the appropriate leading error constants; an uncalibrated raw difference need not be the corrected method's error.

For a vector equation, choose positive component scalesThe absolute and relative error tolerances allow meaningful tests both near zero and for large components with different physical scales. The max-norm criterion $E\le 1$ controls every scaled estimated component. A root-mean-square criterion is another common choice but has a weaker individual-component implication.

If the estimate scales like $E\simeq C{h}^q$, an adaptive step-size controller useswith lower and upper bounds on the change factor. A rejected step leaves the old state and time unchanged and recomputes with a smaller step; an accepted step advances with the chosen output and adjusts the next step. A zero estimate receives capped finite growth rather than division by zero. The safety factor allows for variation in the error constant. Step doubling has $q=p+1$, whereas a $p/(p-1)$ embedded pair usually has $q=p$. If a solver instead controls defect per unit time, division by $h$ changes the scaling exponent and must be included in the controller design.

Using previous estimates can smooth the sequence of accepted steps. A proportional-integral controller can take $h_{n+1}/h_n=\eta E_n^{-a}{E}_{n-1}^b$ with $0<b<a$, appropriate scaling by the estimator order and bounded factors. It avoids overreacting to a single noisy estimate, but needs safeguards after rejected steps, changes of order and abrupt coefficient changes. Event locations and discontinuities require locating the event and restarting the error expansion; blindly applying a smooth-step estimator across a jump can be misleading.

For a multistep method, adaptive step sizes change the interpolation history and coefficients. Constant-step coefficients cannot simply be reused with arbitrary step ratios. Variable-step formulas, compatible starting/restarting procedures and bounded ratios preserve the needed stability; local error and sometimes order are selected together to balance evaluation cost against accuracy.

Stiffness imposes an additional requirement. On a rapidly decaying mode $y^{\prime }=-\kappa y$, explicit Euler requires $h\kappa \le 2$ for absolute stability. Once $y$ is small, a local absolute-tolerance test can propose a large step even though roundoff or a small perturbation would then be amplified. A stiff solver therefore also uses a suitable implicit A-stable or L-stable method and solves its algebraic equations accurately enough for the estimator to be meaningful. A-stability alone does not ensure stiff decay, as the Lobatto stability function in question 4 tends to one for large negative arguments.

Finally, over-tightening tolerance eventually encounters roundoff and algebraic-solve error rather than further truncation-error improvement. Inexact Newton or linear solves must have errors below the local budget; a tiny time step with an inaccurate stage solve does not deliver the formal order. Global error assessment can use controlled refinements, defect propagation or adjoint weighting for a chosen output. Sensible error and step-size control therefore combines asymptotic estimators, tolerance scaling, stability, event handling and implementation error, rather than treating the accepted local estimate as a complete guarantee about the trajectory.

For a discretized ordinary differential equation, error control must distinguish the local truncation error of one exact-start step from accumulated global error. For a smooth order-$p$ one-step method, these are usually $O(h^{p+1})$ and $O(h^p)$ respectively, provided perturbations propagate stably. A small local defect is useful only when the underlying problem and the numerical update do not amplify it excessively.

There are several practical local error estimators. In step-doubling error estimation, compute one step $Y_h$ and two half steps $Y_{h/2,h/2}$ from the same starting value. If the leading method error is $C{h}^{p+1}$, the two-half-step error is $C{h}^{p+1}/2^p$, up to higher terms. Their difference $D=Y_{h/2,h/2}-Y_h$ therefore givesfor exact minus fine-step error. Adding that estimate is Richardson extrapolation. This costs extra evaluations but gives an independent check on a single-method implementation; nonsmooth events can invalidate the assumed expansion.

An embedded Runge-Kutta pair shares stages and uses two output weight vectors to obtain orders $p$ and $q<p$. Their difference generally estimates the lower-order local error $O(h^{q+1})$, not the higher-order error directly. Stage reuse makes this much cheaper than independent integrations. Predictor-corrector error estimation uses the appropriately calibrated difference of two formulas; raw predictor/corrector differences require their known error constants. For multistep schemes, changing the step size also changes interpolation coefficients and requires controlled step ratios and proper history updates, rather than reusing constant-step coefficients.

Component scales should reflect the requested physical accuracy. An example of absolute and relative error tolerances isPositive absolute tolerances avoid singular scaling when a component vanishes, while relative tolerances track widely differing magnitudes. Accept the step if $E\le 1$; otherwise reject it without committing the new state or multistep history. If the estimated defect scales as $h^r$, a basic adaptive step-size controller proposesGrowth limits, a safety factor, and proportional-integral smoothing reduce erratic changes. The exponent must match the actual estimator: $r=p+1$ for a genuine order-$p$ local defect, but $r=p$ for the usual $p/(p-1)$ embedded difference. It is not always $1/p$ or always $1/(p+1)$. A zero estimate uses a finite maximum growth factor, not division by zero.

Local acceptance does not itself certify a global error tolerance. If a dense approximation $\eta (t)$ has residual $d(t)={\eta }^{\prime }(t)-f(t,\eta (t))$ and $f$ is Lipschitz continuous with constant $L$ on the relevant trajectory region, Gronwall's inequality gives the residual bound for global ODE errorThe discrete analogue propagates and sums local defects. A posteriori defect correction, repeat integrations with tighter tolerances, and adjoint-weighted estimates for a chosen final functional provide stronger global checks. A smaller local target may be needed in an unstable flow; trajectory sensitivity should not be mislabeled as a code defect.

For a stiff differential equation, absolute stability can constrain explicit steps far more tightly than accuracy alone. Use an appropriate implicit method and keep nonlinear stage/linear-solve errors below the discretization-error budget. A-stability, L-stability and B-stability address different linear damping and nonlinear contractivity properties. An error controller cannot repair a method outside its stability domain. Tolerance-driven adaptation also does not automatically preserve positivity, conservation laws or a symplectic structure; these may require a suitable method or an additional acceptance test.

Finally, event detection and dense output have their own errors. Locate an event with a consistent interpolant and account for its state sensitivity; a discontinuous right side may require restarting the integrator. Refining indefinitely eventually exposes roundoff error and conditioning limits. Reliable error control balances truncation, solution sensitivity, algebraic-solve error and floating-point error, while reporting failure to meet a target instead of endlessly shrinking a step.

Richardson extrapolation improves an approximation by combining two resolutions so that their leading truncation errors cancel. The useful surprise is that the size of that leading error need not be known: its scaling under refinement is enough. This lets an existing finite difference method or quadrature rule produce a more accurate answer from two runs, without first deriving a new method from scratch.

Suppose $A(h)$ approximates a fixed limit $L$ as the positive resolution parameter $h$ tends to zero, with an asymptotic expansionHere $p$ is known, $c$ is independent of $h$, and $O$ denotes Big O notation. The same asymptotic expansion, with the same $c$, must describe both $h$ and $h/2$. Nested grids make it natural to reuse samples, although nesting itself does not establish this error model. Define the linear combinationThe weights sum to one, preserving $L$, while $2^{p}c(h/2{)}^p-c{h}^p=0$. Thus Richardson extrapolation gives $R(h)=L+O(h^{p+q})$: the next available term in the asymptotic expansion determines the improved order. Justin Solomon's numerical analysis text develops this cancellation for finite differences and applies it to numerical integration.^[ref]

Geometrically, regard the approximation as a function of $t=h^p$. Its first two terms are a straight line $L+ct$. The two computed points lie at positive $t$; extending the line through them to $t=0$ gives $R(h)$. This is why the operation is extrapolation: it estimates a value outside the sampled parameter interval. The remaining curvature corresponds to the uncancelled terms of the asymptotic expansion.

Worked example: a centered derivative of sine. For the smooth function $f(x)=\sin x$, the exact derivative at $x=0$ is $L=\cos 0=1$. The centered difference quotient, a central finite difference, isThe Taylor series of sine gives the signed truncation errorThus $p=2$, $q=2$ and $c=-1/6$. Halving $h$ divides the leading truncation error by four. Combining the two centered difference quotients with Richardson extrapolation givesIndeed, the $h^4$ coefficient becomes $(4/16-1)/(3\cdot 120)=-1/480$; the $h^2$ coefficient vanishes. This is the even-power case of Richardson extrapolation discussed in Persson's authored lecture notes.^[ref]

For $h=0.4$ and $h/2=0.2$, the actual computations areThe absolute error drops from about $2.6454\times 10^{-2}$ for $A(0.4)$ to $5.3080\times 10^{-5}$ for $R(0.4)$, an improvement by roughly a factor of 498. Refining again gives $R(0.2)\approx 0.999996670632607$, whose absolute error is about $3.3294\times 10^{-6}$. The refined truncation errors decrease by nearly four for $A$ and nearly sixteen for $R$, reflecting their second and fourth orders.

The logarithmic scale on each axis turns power laws into straight lines. The dashed reference curves are $h^2/6$ and $h^4/480$, the magnitudes of the leading truncation errors. The plotted errors are computed with 70-digit decimal arithmetic, rather than clipped to a display floor; this isolates the asymptotic expansion from roundoff error.

Worked example: integrating an exponential. Approximate the definite integral of the exponential function on $[0,1]$:For numerical integration, the composite trapezoidal rule adds the areas of trapezoids between successive samples; this is the quadrature use of the endpoint average that also occurs in the ODE trapezoidal rule. On $n$ equal intervals, with $h=1/n$, defineSumming the geometric series gives $T_n=(e-1)\frac{h}{2}\frac{e^h+1}{e^h-1}$. Expanding the exponential function near zero yieldsThe leading truncation error again has $p=2$, $q=2$. On the uniform nested grids with $n$ and $2n$ intervals, Richardson extrapolation therefore givesThis linear combination is exactly composite Simpson's rule on the finer grid, whose spacing is $h/2$. On each coarse panel, an old endpoint receives weight $h/6$ and the new midpoint weight $2h/3$; shared old endpoints receive weight $h/3$. These are the usual Simpson's rule weights $1,4,2,\dots ,4,1$ multiplied by $h/6$. The equivalence is explained directly in LeMesurier's university notes.^[ref]

For two coarse panels and four fine panels,The absolute errors are about $0.0356493$, $0.00894008$ and $0.0000370135$, respectively. Reusing the fine-grid samples makes Richardson extrapolation especially attractive here: the resulting Simpson's rule needs no samples beyond those already used by $T_{2n}$.

When the improvement is justified. Richardson extrapolation requires the stated asymptotic expansion and the correct positive $p$. If the actual leading exponent is $r\ne p$, its coefficient after combination is $c(2^{p-r}-1)/(2^p-1)$, so the leading truncation error generally remains. Merely knowing $A(h)-L=O(h^p)$ does not supply a common coefficient $c$ or a higher-order remainder. A discontinuity of a function, a singular endpoint, or unrelated changes between grids can invalidate the assumed asymptotic expansion; refinement alone does not guarantee improvement for every numerical problem.

Roundoff error and errors in the sampled values survive the cancellation of truncation error. If the coarse and fine values have absolute perturbations bounded by $\epsilon$, the perturbation in $R(h)$ is bounded by $(2^p+1)\epsilon /(2^p-1)$; for $p=2$ this is $5\epsilon /3$. In floating-point arithmetic, a centered difference quotient can also amplify sample errors through division by $h$, so an indefinitely smaller $h$ need not be more accurate. Persson's notes describe this balance between truncation error and roundoff error.^[ref]

The same cancellation idea appears in step-doubling error estimation. For an order-$p$ one-step method for an ordinary differential equation, two half steps each contribute local truncation error of order $h^{p+1}$, so their accumulated leading local error has ratio $2^{-p}$ to one full step. This explains the denominator $2^p-1$ there; the appropriate exponent follows from that local error model, rather than automatically from a generic use of Richardson extrapolation.

For an order-$p$ one-step method with leading numerical error $C{h}^{p+1}$, two half steps have leading error $C{h}^{p+1}/2^p$. The displayed difference estimates exact minus the finer approximation. Adding it implements Richardson extrapolation. It requires a common starting value, the smooth local expansion, and accuracy of all algebraic solves; it can fail across a nonsmooth event.