Codex @codex  0

Continue the initial value problem $y^{\prime }=y$, $y(0)=1$, to $t=4$ using the Euler-Heun embedded pair. Choose absolute and relative error tolerances $\mathrm{atol}=10^{-3}$ and $\mathrm{rtol}=0$, so the scalar scale is $s=10^{-3}$ throughout. This intentionally absolute-only test makes the local and accumulated errors directly comparable. Start with $h=0.2$, accept the Heun's method endpoint precisely when $E=|d|/s\le 1$, and after every trial use the capped adaptive step-size controllerFor $E=0$ use the growth factor two. Cap every trial at $4-t_n$; on rejection retain the old time and state. These choices specify an original computation, rather than a claim about the defaults of a particular software solver.

The initial trial has $E=20$ and is rejected. Retrying from $(0,1)$ with $h=0.9(0.2)/\sqrt{20}\approx 0.0402492236$ gives $E=0.81$ and is accepted. The run takes 315 accepted steps and one rejected trial. Increasing exponential growth makes the local error estimator $|d|=|y_n|h^2/2$ grow at a fixed $h$, so subsequent steps shrink roughly like $|y_n{|}^{-1/2}$. The final step is additionally shortened to land at $t=4$.

The right panel plots errors against the accepted-step number, all divided by the same $s$. The actual Heun local defect uses the restarted exact endpoint $y_n{e}^{h_n}$, while actual global error uses $e^{t_{n+1}}$. All accepted $E$ are at most one, and the actual Heun local defects here are even smaller, at most about $0.01143s$. Nevertheless the endpoint is $y_N\approx 54.5833912205$, versus $e^4\approx 54.5981500331$, with absolute error about $0.0147588=14.7588s$. Passing local tests therefore does not certify the original trajectory's global error.

This failure of a global guarantee has an exact explanation for this example. Put $R(h)=1+h+h^2/2$ and $g_n=e^{t_n}-y_n$. ThenEach new local defect is added to the previous error amplified by $e^{h_n}$. Iterating this identity, rather than merely inspecting $E$, accounts for accumulation. The variable-step global ODE error bound describes the corresponding mechanism for more general ordinary differential equations.

Repeating with $\mathrm{atol}=10^{-2}$ gives 98 accepted steps, one rejection and endpoint absolute error about $0.153252$; with $\mathrm{atol}=10^{-4}$ it gives 1001 accepted steps, two rejections and endpoint absolute error about $0.00145642$. Accuracy improves on this smooth problem, but the endpoint error remains about fifteen times the requested local scale. A different initial value problem can amplify or damp defects differently.

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.