Adaptive ODE solvers and error estimation

Let $d\ge 1$ and $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ be open, and let $f:D\to {\mathbb{R}}^d$ be continuous and locally Lipschitz in the state variable. There is no maximal solution $x:I\to {\mathbb{R}}^d$ of $x^{\mathrm{\prime }}=f(t,x)$ whose interval has a finite right endpoint $b$ and for which some $t_0\in I$, $t_0<b$, and compact set $C\subseteq D$ satisfy $(t,x(t))\in C$ for every $t\in [t_0,b)$. The compact tube condition supplies boundedness of $f$ along the graph and a limit point inside the open domain.

Assume the Axioms of Countable Choice and Dependent Choice. Let $d\ge 1$ and $a<b$, let $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ be open, let $f:D\to {\mathbb{R}}^d$ be continuous, and let $y:[a,b]\to {\mathbb{R}}^d$ be an exact solution of $y^{\mathrm{\prime }}=f(t,y)$ whose graph lies in $D$. Let $U:[a,b]\to {\mathbb{R}}^d$ be componentwise absolutely continuous with graph in $D$. Suppose $\mathrm{\ell },\rho \in L^1[a,b]$, $\rho \ge 0$ almost everywhere, and, for almost every $t$, the points $U(t),y(t)$ satisfy

Define $r(t):=U^{\mathrm{\prime }}(t)-f(t,U(t))$ wherever $U^{\mathrm{\prime }}$ exists, and assume $\parallel r(t){\parallel }_2\le \rho (t)$ almost everywhere. Then, for every $t\in [a,b]$,

The bound is a posteriori: it uses the residual of the actual reconstruction. An embedded Runge–Kutta difference alone does not supply the hypothesis $\parallel r(t){\parallel }_2\le \rho (t)$.

Let $d\ge 1$ and $a<T$, let $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ be open, let $f:D\to {\mathbb{R}}^d$ be continuous, and let $y$ be an exact solution of $y^{\mathrm{\prime }}=f(t,y)$ on $[a,T]$. Let $U:[a,T]\to {\mathbb{R}}^d$ be continuous and $C^1$ on each member of a finite partition (continuously differentiable on each closed panel, with one-sided endpoint derivatives), with graph in $D$. Let $\mathrm{\ell },\rho :[a,T]\to [0,\mathrm{\infty })$ be continuous on each partition interval, with one-sided endpoint values, and suppose on every open panel $⟨U-y,f(t,U)-f(t,y)⟩\le \mathrm{\ell }(t)\mathrm{\parallel }U-y{\mathrm{\parallel }}_2^2,\mathrm{\parallel }{U}^{\mathrm{\prime }}(t)-f(t,U(t)){\mathrm{\parallel }}_2\le \rho (t)\mathrm{.}$ Then for every $t\in [a,T]$, $\mathrm{\parallel }U(t)-y(t){\mathrm{\parallel }}_2\le e^{{\int }_a^t\mathrm{\ell }(v)dv}\mathrm{\parallel }U(a)-y(a){\mathrm{\parallel }}_2+{\int }_a^t{e}^{{\int }_s^t\mathrm{\ell }(v)dv}\rho (s)ds\mathrm{.}$ The proof uses only finite-piece ordinary calculus and no choice axiom. The general absolutely continuous case is treated separately.

For an ODE $y^{\mathrm{\prime }}=f(t,y)$ on states in ${\mathbb{R}}^d$ with $d\ge 1$, two explicit Runge–Kutta weight vectors $b=(b_i{)}_{i=1}^s$ and $\widehat{b}=({\widehat{b}}_i{)}_{i=1}^s$ are an embedded pair when they use the same explicit stage recursion $K_i=f\text{\hspace{0.17em}⁣}(t_n+c_{i}h,y_n+h{\sum }_{j<i}{a}_{ij}{K}_j),1\le i\le s,$ but produce two one-step values $Y=y_n+h{\sum }_{i=1}^s{b}_i{K}_i,\widehat{Y}=y_n+h{\sum }_{i=1}^s{\widehat{b}}_i{K}_i\mathrm{.}$ The two methods are the high and low members when their local orders are respectively at least $p+1$ and $p$, for an integer $p\ge 1$. Their embedded difference estimator is the nonnegative scalar $E_n(h):=\mathrm{\parallel }\widehat{Y}-Y{\mathrm{\parallel }}_2,$ with the Euclidean norm. This definition asserts only how the estimator is computed; its reliability as a bound for either local defect is a separate claim and is not part of the definition.

The stage recursion and the local defect sign convention are those of Explicit Runge–Kutta tableaus, stage maps, and failure status and One-step ODE methods, local defect, uniform order, and step stability.

Inputs are an IVP in the sense of First-order systems, initial value problems, and solutions on intervals, $y^{\mathrm{\prime }}=f(t,y)$ on an open domain $D\subseteq \mathbb{R}\times {\mathbb{R}}^d$ with initial data $(t_0,y_0)\in D$ and requested interval $[t_0,T]$, where $t_0\le T$ and $d\ge 1$; an initial proposal $0<h_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}\le h_{\max }$; an embedded pair with low order $p\ge 1$ as in Embedded Runge–Kutta pair and difference estimator; tolerances $\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{l}>0$ and $\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{l}\ge 0$; safety $0<\sigma <1$; scale bounds $0<{\gamma }_{\min }\le 1\le {\gamma }_{\max }$; step bounds $0<h_{\min }\le h_{\max }$; and a positive integer attempt budget $N$.

Output. Return a record consisting of one of success, failure, minimum-step-rejected, or attempt-budget-exhausted, together with the accepted mesh $\mathcal{A}=((t_0,y_0),\dots ,(t_m,y_m))$. The mesh always contains the initial pair and every accepted endpoint; on success, $t_m=T$ and $y_m$ is the returned endpoint approximation. On any other status, $\mathcal{A}$ is the accepted prefix and $y_m$ is the last accepted approximation.

Every maximum and minimum below is over a displayed nonempty finite list of real numbers; existence follows from Every nonempty finite set of reals has a maximum and a minimum.

Procedure. Initialize $t:=t_0$, $y:=y_0$, $\mathcal{A}:=((t_0,y_0))$, and proposal $h_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}}:=h_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}$. If $t_0=T$, return success with $\mathcal{A}$. For each attempt index $k=1,\dots ,N$, perform exactly the following steps:

1. Set $h:=\min (h_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}},h_{\max },T-t)$. The last term permits a final clipped step shorter than $h_{\min }$. If $h\le 0$ before reaching $T$, return failure with $\mathcal{A}$.
2. Evaluate the shared stages in tableau order, checking before each evaluation that its time-state input lies in $D$. If a stage input is outside $D$, its evaluation is undefined, or its value is nonfinite, return failure with $\mathcal{A}$. Compute both pair outputs; if either output or any arithmetic needed to compute them is undefined, fails, or is nonfinite in the rounded realization, return failure with $\mathcal{A}$.
3. Define $s_i:=\mathrm{a}\mathrm{t}\mathrm{o}\mathrm{l}+\mathrm{r}\mathrm{t}\mathrm{o}\mathrm{l}\max (\mathrm{\mid }(y{)}_i\mathrm{\mid },\mathrm{\mid }(Y{)}_i\mathrm{\mid })$ and $Q:={\max }_i\mathrm{\mid }(\widehat{Y}-Y{)}_i\mathrm{\mid }\mathrm{/}{s}_i$, where $Y$ is the high-order value and $\widehat{Y}$ the low-order value. If any scale or $Q$ operation fails or is nonfinite in the rounded realization, return failure with $\mathcal{A}$. If $Q\le 1$, first check $(t+h,Y)\in D$; if not, return failure with $\mathcal{A}$. Otherwise accept: set $t:=t+h$, $y:=Y$, and append $(t,y)$ to $\mathcal{A}$. If $Q>1$, reject and retain $t,y,\mathcal{A}$.
4. If the attempt was accepted and $t=T$, return success with $\mathcal{A}$. Otherwise, if it was rejected and $h\le h_{\min }$, return minimum-step-rejected with $\mathcal{A}$. These terminal checks precede the budget check, so a final accepted step that reaches $T$ succeeds and a rejection at the minimum step reports that safeguard.
5. If $k=N$, return attempt-budget-exhausted with $\mathcal{A}$. Otherwise compute $g=\{\begin{array}{ll}{\gamma }_{\max },& Q=0,\\ \min ({\gamma }_{\max },\max ({\gamma }_{\min },\sigma \mathrm{/}\sqrt[p+1]{Q})),& Q>0.\end{array}$ Set $h_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}}:=\min (h_{\max },\max (h_{\min },gh))$ and continue to attempt $k+1$. A failed or nonfinite rounded computation of $g$ or $h_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{p}}$ returns failure with $\mathcal{A}$.

For $Q>0$, $\sqrt[p+1]{Q}$ is the unique positive $(p+1)$st root supplied by Existence and uniqueness of $n$-th roots: a unique $a^{1\mathrm{/}n}\ge 0$ with $(a^{1\mathrm{/}n}{)}^n=a$. The exponent $1\mathrm{/}(p+1)$ is a leading-error model for step selection; neither acceptance nor this controller is a global error certificate or a stiffness test.

The mathematical specification evaluates each displayed operation, including maximum, minimum, absolute value, and positive-root evaluation, exactly. A rounded implementation instead applies its fixed finite-format rounding rule to every stage, output, scale, and controller-factor arithmetic operation and uses the resulting stored high and low values in the displayed difference and acceptance test; nonfinite or failed operations return failure. The finite termination theorem applies to either realization when each attempt has finite work. The rounded realization has no a posteriori guarantee: a tolerance below its resolution can be accepted because both stored values coincide.

For $y^{\mathrm{\prime }}=f(t,y)$, take the shared explicit stages $K_1=f(t_n,y_n),K_2=f(t_n+h,y_n+h{K}_1)\mathrm{.}$ The Euler member is $Y_E=y_n+h{K}_1$ and the Heun explicit trapezoidal member is $Y_H=y_n+\frac{h}{2}(K_1+K_2)$. Their embedded difference is $Y_H-Y_E=\frac{h}{2}(K_2-K_1)$. The local orders are proved separately. In Embedded Runge–Kutta pair and difference estimator notation, the common tableau has $c=(0,1)$, $A_{21}=1$, high weights $(1\mathrm{/}2,1\mathrm{/}2)$ and low weights $(1,0)$.

Fix one initial value problem and initial state. Let $D_p(h):=\widehat{Y}(h)-y(t+h)$ be the low member's one-step defect and $D_{p+1}(h):=Y(h)-y(t+h)$ the high member's defect, using the same exact initial state and step start. Suppose, as $h\downarrow 0$, $D_p(h)=a{h}^{p+1}+O(h^{p+2}),D_{p+1}(h)=O(h^{p+2})$ for a fixed vector $a\in {\mathbb{R}}^d$. Here $O(h^r)$ for a vector means a remainder whose Euclidean norm is at most $C{h}^r$ for all sufficiently small positive $h$. Then $\widehat{Y}(h)-Y(h)=a{h}^{p+1}+O(h^{p+2})$. If $a\mathrm{\ne }0$, the estimator satisfies $E(h)\mathrm{/}\mathrm{\parallel }{D}_p(h){\mathrm{\parallel }}_2\to 1$. If $a=0$, these hypotheses do not imply that this ratio tends to one or that the estimator is nonzero; the Euler–Heun pair on $y\mathrm{\prime }=t^2$, $y(0)=0$, has ratio $3\mathrm{/}2$ as computed below.

For one Heun step from $(t_n,y_n)$ of size $h$ with stages $K_1,K_2$ as in Euler–Heun embedded Runge–Kutta pair, define its quadratic continuous extension for $0\le \theta \le 1$ by $P_n(\theta ):=y_n+h(\theta K_1+\frac{1}{2}{\theta }^2(K_2-K_1))\mathrm{.}$ It has $P_n(0)=y_n$ and $P_n(1)=Y_H$, the accepted Heun endpoint. It is a separate within-step approximation; its order is analyzed independently of the endpoint order.

Fix a dimension $d\ge 1$, rational times $a<T$, a rational initial vector $y_0$, and an ODE field $f$ represented by a finite acyclic expression DAG with rational constants, time/state inputs, and only $+,-,\times ,÷$ nodes. Its domain consists of points where every division is defined. Write $D_{y}f$ for the matrix of partial derivatives with respect to the state coordinates. A finite rational residual-tube certificate consists of the following exact rational data:

1. A finite partition $a=t_0<t_1<\cdots <t_m=T$. On each panel $[t_j,t_{j+1}]$, a vector polynomial $U_j(t)$ with dyadic (hence rational) coefficients, whose neighboring endpoint values agree exactly. These stored coefficients define the actual reconstruction, including any endpoint correction or rounded stage value used by an implementation.
2. For each panel, a closed rational state box $X_j\subset {\mathbb{R}}^d$, rational bounds $L_j\ge 0$, ${\rho }_j\ge 0$, a rational $d_j\ge 0$, and a polynomial radius $R_j(t)$ with rational coefficients. The radius polynomials agree at neighboring endpoints. For a panel's finite rational subdivision, the intended inequalities are $R_j(t)>0$, $R_j^{\mathrm{\prime }}(t)-L_j{R}_j(t)-{\rho }_j>0$, and $U_j(t)+[-R_j(t),R_j(t){]}^d\subseteq X_j$.
3. A finite rational subdivision of every panel, with exact endpoint coverage, no gaps, and pairwise-disjoint interiors. It supplies cells on which each required range test below is performed. Rational cell endpoints and constants are enclosed by directed outward rounding to finite dyadic endpoints before interval evaluation; an unrepresentable or overflowing conversion cannot validate.
4. An initial check target $d_0\ge 0$ with ${\sum }_i(U_0(a{)}_i-(y_0{)}_i{)}^2\le d_0^2$ and $d_0<R_0(a)$.

The certificate is validated only if outward-rounded real-interval evaluation on every supplied cell proves all three radius/tube inequalities, the range inclusion in $X_j$, that every division denominator in the field and its state derivatives excludes zero on the time-state box $[t_j,t_{j+1}]\times X_j$, and ${\sum }_{i=1}^d\sup \mathrm{\mid }{U}_{j,i}^{\mathrm{\prime }}-f_i(t,U_j(t)){\mathrm{\mid }}^2\le {\rho }_j^2,{\sum }_{i,k=1}^d\sup \mathrm{\mid }{\mathrm{\partial }}_{y_k}{f}_i(t,y){\mathrm{\mid }}^2\le L_j^2$ there. Here each supremum is replaced in the check by the maximum absolute endpoint of a finite outward interval enclosure; the residual expression is evaluated at the actual polynomial $U_j$, so interpolation, stage, and storage errors already present in that reconstruction are included. Strict inequalities must have a certified positive interval lower endpoint; an inconclusive, nonfinite, unsupported, or failed test is not validation.

This format certifies only rational expression fields and polynomial reconstructions using the stated finite checks. It uses real intervals, not complex rectangles or balls, and does not claim an evaluator for transcendental functions or general library floating-point operations. The page proves the soundness of a validated certificate separately.

The finite DAG and real-interval conventions are those of Real intervals, complex rectangles and ball enclosures and Outward-rounded interval, rectangle and ball operations.

Global error from uniform one-step defect and stability assumes a prescribed uniform mesh $t_n=a+nh$, a one-step defect bound on that mesh, and a stability estimate for propagation on those steps. Its conclusion therefore applies when those hypotheses hold. The controller in Finite-budget adaptive Runge–Kutta controller accepts variable, state-dependent step sizes based on a pair difference. That difference is not the exact one-step defect hypothesis of the fixed-mesh theorem. No uniform mesh, defect reliability estimate, mesh-ratio control, or variable-step propagation theorem is supplied by that theorem. It may be used on a particular adaptive run only after those distinct hypotheses have independently been proved; its fixed-mesh statement alone yields no adaptive global error certificate.

For finite-dimensional states and the finite-budget procedure in Finite-budget adaptive Runge–Kutta controller, the controller returns one of success, failure, minimum-step-rejected, or attempt-budget-exhausted after at most $N$ trial attempts. Each individual attempt uses finitely many stage evaluations and arithmetic operations.

Let $d\ge 1$ and let $f$ be $C^2$ on an open neighborhood of a compact convex time-state set containing the initial point, the short exact solution segment, and the Euler stage points. Assume the first and second derivatives of $f$ are bounded there. For the exact solution $y^{\mathrm{\prime }}=f(t,y)$ starting from $y(t)=y_0$, put $g:={\mathrm{\partial }}_{t}f(t,y_0)+D_{y}f(t,y_0)f(t,y_0)\mathrm{.}$ Here ${\mathrm{\partial }}_{t}f$ is the time derivative and $D_{y}f$ is the Jacobian in the state variables. The notation $O(h^k)$ means a norm bounded by a constant times $h^k$ for sufficiently small positive $h$, uniformly over the compact starting set when stated. For all sufficiently small $h>0$, the Euler and Heun one-step values in Euler–Heun embedded Runge–Kutta pair satisfy $Y_E-y(t+h)=-\frac{1}{2}{h}^{2}g+O(h^3),Y_H-y(t+h)=O(h^3),Y_H-Y_E=\frac{1}{2}{h}^{2}g+O(h^3)\mathrm{.}$ Moreover, uniformly for $0\le \theta \le 1$, $y(t+\theta h)=y_0+\theta h{K}_1+\frac{1}{2}{\theta }^2{h}^{2}g+O(h^3),K_1=f(t,y_0)\mathrm{.}$ For a family of starting points, assume one common $h_0>0$ for which all exact segments and stage points stay in the same neighborhood. Under this condition the estimates are uniform over the compact starting set because the derivative bounds are common. Thus Euler has local order one and Heun has local order at least two on that set. If $g\mathrm{\ne }0$, the embedded norm estimates Euler's defect asymptotically with ratio one; when $g=0$, that ratio conclusion does not follow from these orders.