Codex @codex  0

If a scaled local error estimator behaves as $E\propto h^r$ with approximately unchanged error scale, the displayed update targets a fixed local level. A safety factor $0<\eta <1$ aims below the acceptance boundary. For an embedded Runge-Kutta pair with order of a numerical method $p$ for the accepted output and $p-1$ for the companion formula, the difference generally has $r=p$; for the actual unscaled local truncation error of an order-$p$ formula, $r=p+1$. The exponent belongs to the quantity being controlled. Dividing the defect by $h$ changes that exponent. This follows by solving $E(h_{\mathrm{new}})\approx E(h)(h_{\mathrm{new}}/h{)}^r$ for the next step.^[ref]

Rejecting a trial means retrying from the same accepted time and state, with a smaller step. Advancing time or retaining the rejected output would change the initial value problem being approximated. Minimum and maximum change factors prevent extreme updates; a zero estimate requires finite capped growth. After acceptance, cap the next step at the remaining interval length. Detect a step that cannot advance time in floating-point arithmetic, and report failure rather than an endless retry. An adaptive step-size controller can use earlier estimates to smooth step sequences, but the memoryless rule here is sufficient to explain acceptance and rejection.^[ref]

A Butcher tableau records a Runge-Kutta method's stage matrix $A$, time nodes $c$, and output weights $b$ in the block arrayRead row $i$ above the line as the recipe for Runge-Kutta stage $i$: the left entry gives $t_n+c_{i}h$, and the right entries give the coefficients in $Y_i=y_n+h{\sum }_j{a}_{ij}{k}_j$. Below the line, $b_i$ weights slope $k_i$ in the final update. Zeros on and above the diagonal identify an Explicit Runge-Kutta method. For the internally consistent convention, row sums of $A$ equal $c$; first-order consistency of a numerical method also requires ${\sum }_i{b}_i=1$. An embedded Runge-Kutta pair adds another weight row while retaining $A,c$.

Two Runge-Kutta methods share a stage matrix and nodes but use different output weights. An order-$p$ and order-$q$ pair with $q<p$ gives an output difference generally of order $h^{q+1}$, estimating the lower-order defect. Shared stage evaluations reduce the cost. The adaptive step-size controller exponent must reflect this estimator order, rather than blindly the higher reported solution order.

Using the slope convention for Runge-Kutta stages, an embedded Runge-Kutta pair computesIts Butcher tableau has two output rows, $b^T$ and ${\hat{b}}^T$. Evaluating $f$ is often the expensive part of solving an ordinary differential equation; after the shared slopes have been computed, the extra row needs only another linear combination.^[ref]

Write $\Phi$ for the exact endpoint from the current numerical state, and define signed local defects ${\ell }_q=\Phi -Y^{[q]}$ and ${\ell }_p=\Phi -Y^{[p]}$. Then $d={\ell }_q-{\ell }_p$. If ${\ell }_q=C_q{h}^{q+1}+O(h^{q+2})$ and ${\ell }_p=O(h^{p+1})$, with $C_q\ne 0$, the difference reproduces the leading lower-order defect. It does not reproduce the leading higher-order defect. Here the powers use big O notation and a smooth Taylor expansion; a vanishing leading coefficient or accidental cancellation can make the estimate misleading.

A common choice accepts $Y^{[p]}$ after testing $d$ with absolute and relative error tolerances. With the scale approximately unchanged on a retry, the adaptive step-size controller uses exponent $1/(q+1)$. Accepting a better endpoint does not change the power in the quantity being tested. The Euler-Heun embedded pair makes this visible without a large coefficient table.^[ref]

Error control selects discretization parameters and checks approximations against requested accuracy. Local error estimators, absolute and relative error tolerances and an adaptive step-size controller provide practical local control. Propagation of defects determines the global error, while stiffness, roundoff error and algebraic-solve accuracy impose additional limits.

For an ordinary differential equation, the practical loop is to propose a step, form a local error estimator, scale it by absolute and relative error tolerances, and accept or retry before choosing the next step. An embedded Runge-Kutta pair supplies two endpoints from shared Runge-Kutta stages, making this extra accuracy check inexpensive. The Euler-Heun embedded pair below exposes the distinction between the accepted formula's order of a numerical method and the estimator's step exponent; its Adaptive Euler-Heun integration of exponential growth shows accepted local tests alongside accumulated global error.

A local error estimator approximates a numerical step's defect from the exact solution starting at the same input. Differences between formulas estimate error only after their asymptotic orders and constants are accounted for. Step-doubling error estimation, an embedded Runge-Kutta pair and predictor-corrector error estimation are common constructions. They do not by themselves bound accumulated global error.

Let ${\Phi }_h(t_n,y_n)$ be the exact endpoint of the ordinary differential equation started at the numerical value $y_n$, and let $Y_{n+1}$ be the candidate numerical endpoint. The signed local defect being estimated is ${\ell }_{n+1}={\Phi }_h(t_n,y_n)-Y_{n+1}$. In contrast, actual global error compares $Y_{n+1}$ with the solution of the original initial value problem, including all previous steps. The local truncation error convention inserts original exact data instead. These starting-value distinctions matter even when the same formula is used for every step. A practical local error estimator also relies on a smooth Taylor expansion; a small difference can reflect accidental cancellation rather than a rigorous upper bound.

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.

A Runge-Kutta method advances an initial value problem $y^{\prime }=f(t,y)$, $y(t_0)=y_0$, by sampling its vector field at several trial states within one time step. These samples are Runge-Kutta stages; a weighted combination of their slopes produces the next state. The construction improves on the single tangent line used by the Euler method while requiring evaluations of $f$, rather than explicit higher derivatives of the solution. It applies to scalar equations and systems of ordinary differential equations alike. The coefficient choices distinguish methods such as Heun's method, the classical fourth-order Runge-Kutta method and implicit Runge-Kutta methods.

For $s$ stages, write $A=(a_{ij})$ for the stage matrix, $c_i$ for relative time nodes and $b_i$ for output weights. A time step of length $h$ computesHere $k_i$ is a slope, with units of state per time; $h{k}_i$ is a state increment. Another common convention writes $K_i=h{k}_i$, so its stage states are $Y_i=y_n+{\sum }_j{a}_{ij}{K}_j$ and its output is $y_n+{\sum }_i{b}_i{K}_i$. These are equivalent conventions: including $h$ both inside $K_i$ and again in the output would count the time step twice. The textbook formulation uses the increment convention.^[ref]

Usually the Butcher tableau obeys $c_i={\sum }_j{a}_{ij}$, called internal consistency. Stage $i$ then combines slopes over an effective time $c_{i}h$. The second-order conditions with independent Runge-Kutta abscissae explain what changes when time nodes are specified independently. A Runge-Kutta stage is generally an auxiliary approximation, not a separate accepted solution value at its node. Several stages can even have the same time node but different state arguments, as in the classical fourth-order Runge-Kutta method.

The useful idea is to approximate the integral $y(t_n+h)-y(t_n)={\int }_{t_n}^{t_n+h}f(t,y(t))dt$ without knowing the intermediate exact trajectory. Trial states provide places to sample $f$, and Butcher order conditions make the resulting combination agree with the exact Taylor expansion. More stages alone do not guarantee a higher order of a numerical method: their coefficients must satisfy the appropriate conditions. Embedded Runge-Kutta pairs reuse these same samples for a second output and a practical local error estimator.