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.

The Euler-Heun 1(2) pair combines the explicit Euler method and Heun's method. Here 1(2) names their respective order of a numerical method; the accepted endpoint in this article is the second-order one. Its shared Runge-Kutta stages and endpoints areThe first output row in this Butcher tableau is Heun's method; the second is the explicit Euler method. Both outputs cost two evaluations of $f$ together, rather than one for Euler method plus two for a separately computed Heun's method.

For the initial value problem $y^{\prime }=y$, $y(0)=1$, take $h=0.2$. The first slope is $k_1=1$; the stage state is $1+0.2k_1=1.2$, giving $k_2=1.2$. ThereforeThe exact exponential function gives signed local defects approximately $0.02140275816$ for the explicit Euler method and $0.00140275816$ for Heun's method. Their difference is exactly $0.02$. Thus the pair difference estimates the larger Euler defect; calling it the actual Heun defect would miss by a factor about fourteen here.

More generally, a step from $y_n$ on $y^{\prime }=y$ hasThese original Taylor expansions separate three quantities: the accepted Heun's method has local truncation error $O(h^3)$ and global error $O(h^2)$ under suitable numerical stability, while the tested local error estimator is $O(h^2)$. For $E=|d|/s$ with $s$ held approximately fixed, use $h_{\mathrm{new}}=\eta h{E}^{-1/2}$. A cube-root update would model a different quantity, the actual higher-order local defect.

A German mathematician who studied linear differential equations and later held a chair in technical mechanics at Karlsruhe. His work on second-order equations with four singular points led to the Heun equation, extending the family associated with the hypergeometric differential equation.

His name is also attached to Heun's method, a predictor-corrector Runge-Kutta method for numerical integration. The predictor estimates an endpoint slope and the correction averages it with the initial slope. The special-function equation and the numerical method are distinct contributions linked by the study of differential equations, rather than two names for the same construction. ^[ref]

Homepage: mathshistory.st-andrews.ac.uk/Biographies/Heun/

Local errors of the Euler schemes. Let $h=\Delta t>0$, and evaluate $f$ and its partial derivatives at an exact starting point $(y(t),t)$. DefineThe exact Taylor expansion is $y(t+h)=y+hf+\frac{1}{2}{h}^{2}A+O(h^3)$. The forward Euler method retains only $y+hf$, so its one-step defect, exact minus computed value, isFor the Backward Euler method, substitution of the exact solution into the scheme gives the residualbecause $f(y(t+h),t+h)=f+hA+O(h^2)$. Under local smoothness, solving the implicit step changes this residual by a factor $1+O(h)$, so the one-step defect has the same leading term:These are unnormalized local truncation errors, as in the question. A convention dividing by $h$ instead gives $O(h)$. Both schemes have first-order global error under the usual smoothness and stability hypotheses.

Euler stability. For the Dahlquist test equation $y^{\prime }=-\Lambda y$, put $z=h\Lambda$, with $\mathrm{Re}z>0$. The numerical amplification factors areStrict decay requires $|R|<1$; the boundary $|R|=1$ gives bounded but generally nondecaying values. For Forward Euler method,For positive real $\Lambda$, this reduces to $0<h\Lambda <2$. The absolute stability region in the $z$ plane is the disk $|z-1|\le 1$. For Backward Euler method, $|1+z|>1$ whenever $\mathrm{Re}z>0$, so it decays for every step size. It is A-stable and also L-stable, since $R_{BE}(z)\to 0$ for large $z$ in this half-plane. Here $z$ has the opposite sign from the alternative convention $z=h\lambda$ for $y^{\prime }=\lambda y$.

Predictor-corrector accuracy. The printed two-stage construction is Heun's method, equivalentlyAlthough the second stage resembles a Backward Euler method evaluation, it uses the predicted value, so the combined method is explicit. Expand that stage around the exact starting point, puttingThenand the numerical step is $y+hf+\frac{1}{2}{h}^{2}A+\frac{1}{4}{h}^{3}B+O(h^4)$. Along the exact solution, $y^{\prime \prime \prime }=B+f_{y}A$, so comparison with the exact Taylor expansion givesThus the local defect improves by one power of $h$, giving second-order global error under the usual hypotheses. For the linear test equation the defect reduces to $-h^3{\Lambda }^{3}y/6+O(h^4)$.

Predictor-corrector stability. Applied to the test equation, the predictor is $y^{(1)}=(1-z)y_n$ and the second value is $y^{(2)}=(1-z+z^2)y_n$. Their average givesFor real positive $z$, $R_H(z)=\frac{1}{2}[(z-1{)}^2+1]$, so strict decay holds exactly for $0<z<2$, the same real-axis step bound as Forward Euler method. For complex $z=x+iy$ the precise condition isThe region is bounded and does not contain the full right half-plane: for example $z=3$ gives $R_H=5/2$. It is therefore not A-stable. Its imaginary-axis values satisfy $|R_H(iy){|}^2=1+y^4/4>1$ for $y\ne 0$. Better accuracy has not inherited the unconditional stability of an actual implicit Backward Euler method 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.