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]

For a one-step approximation $y_{n+1}=y_n+h\varphi (t_n,y_n,h)$ to an ordinary differential equation, insert the exact solution and define the unscaled local truncation errorConsistency of a numerical method requires ${\tau }_{n+1}/h\to 0$ uniformly on the finite interval being approximated, under the stated smoothness assumptions. Merely asking ${\tau }_{n+1}\to 0$ is too weak: the useless update $y_{n+1}=y_n$ has defect $h{y}^{\prime }(t_n)+O(h^2)$, which tends to zero but generally fails the scaled test. Here $O$ denotes Big O notation. MIT uses this division-by-step convention.^[ref]

For the Euler method, Taylor's theorem gives ${\tau }_{n+1}=h^2{y}^{\prime \prime }({\xi }_n)/2$. A bounded second derivative therefore gives local truncation error $O(h^2)$ and a scaled defect $O(h)$. Under a uniform Lipschitz continuity assumption on $f$ in $y$, within a region containing the exact and numerical solutions, numerical stability controls accumulation over about $T/h$ steps. This yields global error $O(h)$ with exact starting data. Consistency of a numerical method checks each step; numerical stability controls their combination.

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.

The Euler method approximates an initial value problem $y^{\prime }=f(t,y)$, $y(t_0)=y_0$, by evaluating the derivative at the current point and moving along that tangent line for a time $h$. With $t_n=t_0+nh$, its displayed update needs one evaluation of $f$ per step. This simple rule turns local slope information into a whole approximate trajectory. MIT's course notes derive it by dropping the second-order Taylor remainder.^[ref]

Take $y^{\prime }=-2y$, $y(0)=1$, whose exact solution is the exponential function $y(t)=e^{-2t}$. The Euler method givesAt $h=0.2$, the first values are $1,0.6,0.36,0.216$. At $t=1$, five steps give $0.6^5=0.07776$, compared with $e^{-2}\approx 0.135335$. The result decays but is too small: absolute stability alone does not ensure accuracy. At $h=1.2$, the first values are $1,-1.4,1.96,-2.744$, and at $t=6$ five steps give $-5.37824$, while the exact solution is about $0.00000614421$.

For the Dahlquist test equation $y^{\prime }=\lambda y$, the Euler method has stability function $R(z)=1+z$ with $z=h\lambda$. Its absolute stability region is $|1+z|\le 1$, the closed disk centred at $-1$ with radius one. For the present real decay problem and $h>0$, bounded numerical values require $0<h\le 1$; strict decay requires $0<h<1$. At $h=1$ the values alternate $1,-1,1,\dots$ without decay. These conditions follow directly from the amplification factor $1-2h$, matching MIT's criterion for $y^{\prime }=-ay$.^[ref]

The Euler method is first order: with bounded second derivative and suitable Lipschitz continuity, its global error is $O(h)$ on a fixed interval. In this example the logarithm gives $\log (1-2h)=-2h-2h^2+O(h^3)$ for small positive $h$. At $t=nh$, this means $y_n=e^{-2t}(1-2th+O(h^2))$. It explains both the downward bias and why halving the step roughly halves the error once $h$ is sufficiently small.

For a smooth solution of $y^{\prime }=f(x,y)$, Taylor's theorem gives the displayed local truncation error when an Euler method step starts at an exact solution value. Thus the one-step defect is $O(h^2)$. Dividing it by $h$ gives a normalized local truncation error of $O(h)$; these are different conventions for reporting the same residual. Also $y^{\prime \prime }=f_x+f_{y}f$ along the solution by the chain rule. A global $O(h)$ error conclusion additionally requires the usual regularity and stability control over successive steps.

The coefficients of this linear multistep method are ${\alpha }_0=b$, ${\alpha }_1=-a$, ${\alpha }_2=1$ and ${\beta }_0={\beta }_1=0$, ${\beta }_2=c$. Substitution into the order of a numerical method conditions givesThe second equation gives $c=2-a$; substituting this into the third gives $-2+3a/2=0$. Consequently there is exactly one second-order choice,The resulting scheme is the BDF2 method,For a $C^4$ exact solution, Taylor's theorem in the original normalization gives the exact leading order defectIn particular the coefficient $-2/9$ is nonzero: the method has formal order exactly two, rather than three. Bounded third derivatives already suffice for the uniform estimate $T_n(h)=O(h^2)$; the stronger smoothness above only makes the displayed remainder immediate. Multiplying the recurrence by $3/2$ instead gives the common derivative normalization with leading defect $-h^2{y}^{\prime \prime \prime }(x_n)/3$.

The first of the characteristic polynomials of a linear multistep method isBoth roots are simple and lie in the closed unit disk, so the root condition for a multistep method proves zero-stability.

The Dahlquist equivalence theorem states that, for fixed-coefficient linear multistep methods on the usual class of initial value problems with a right-hand side continuous in $x$ and uniformly Lipschitz continuous in $y$, convergence of a numerical method for every such problem and every convergent set of starting values is equivalent to consistency of a numerical method and zero-stability. The statement is on fixed finite intervals where the exact solution exists, with implicit updates well defined. Its order estimate says that a zero-stable method with $T_n(h)=O(h^p)$ and starting errors $O(h^p)$ has global error $O(h^p)$ under the corresponding solution smoothness and Lipschitz continuity assumptions; merely convergent starting values do not ensure the order-$p$ rate.

For this application, sufficient assumptions are a unique solution $y\in C^3([x_0,x_0+T])$, a uniform Lipschitz continuity bound $L$ in $y$ on a neighborhood of its graph, and the nearby solution branch of each implicit update. For sufficiently small $h$, the mapis a contraction mapping with factor at most $2hL/3<1$ on a suitable small closed neighborhood. The Banach fixed-point theorem supplies local solvability of the nearby update. With starting errorsthe order estimate givesThus the BDF2 method is second-order convergent under these standard hypotheses. Exact $Y_0=y(x_0)$ together with one Euler method step for $Y_1$ already gives the required starting error when $y^{\prime \prime }$ is bounded, since only one starting step is taken.

The question's assumption of continuity of $f$ alone does not establish uniqueness, the needed solution smoothness, or this uniform convergence estimate. The added Lipschitz continuity condition is a sufficient hypothesis for the theorem used here; its failure does not imply that every particular initial value problem is nonunique or that its numerical approximations cannot converge.

If one retains the Euler approximation $A=I+J\delta t$, its exact norm matrix isThe additional term can mix the instantaneous symmetric-part eigenvectors, so they no longer generally diagonalize finite-time error amplification. However, when $\delta t^2$ is material, a reliable physical propagator must itself be calculated beyond first-order Euler: a large Euler-step amplification can be a numerical artifact. For constant $J$, the exact map is $\Phi =e^{J\Delta t}$; for a changing reference trajectory it solves $\dot{\Phi }=J(t)\Phi$, $\Phi (0)=I$.

Use its singular value decomposition, $\Phi =U\Sigma V^T$. The initial orthonormal right singular vectors ${\mathbf{v}}_i$ diagonalize ${\Phi }^T\Phi$ and satisfyA singular value $s_i>1$ gives finite-time norm growth, $s_i<1$ decay, and the leading right singular vector is the initial perturbation with greatest amplification in the stated norm. The left singular vector describes the resulting final error pattern. Non-normal dynamics can give transient growth even if every Jacobian matrix eigenvalue has negative real part, so selecting observations from those eigenvalues alone is insufficient. Over a genuinely short interval, $s_i^2=1+{\lambda }_i\delta t+O(\delta t^2)$, recovering the symmetric-part result. Initial observations should target the relevant growing right singular projections, with prior covariance and observational precision included when they are not isotropic.

Let the step be $h$. Forward Euler method gives $x_1=x_0+h{v}_0$, $v_1=v_0+hF(x_0)$. Applying the same map with $-h$ gives $x_2=x_1-h{v}_1=x_0-h^{2}F(x_0)$, which generally differs from $x_0$. It is therefore not its own inverse under time reversal.

A kick–drift–kick leapfrog integrator instead hasThe step $-h$ first recovers $v_{1/2}$, then $x_0$, then $v_0$, so $\boxed{{\Phi }_{-h}={\Phi }_h^{-1}}$. This assumes a position-dependent force and the same fixed step. Newtonian mechanics is reversible; respecting that structure helps avoid artificial directional bias and often improves long-term qualitative behaviour. The usual leapfrog also has a structure preserving symplectic geometry, which supports bounded oscillatory energy error, but reversibility by itself does not imply exact conservation of energy. Adaptive stepping requires special care to preserve these properties.

A backward Euler method with spatial second-order central differences givesThe source is fixed in pseudo-time; if it were evolving, its evaluation would also need specifying. Each pseudo-time step therefore requires a linear system solve rather than a purely local explicit update.

The characteristic polynomial is ${\lambda }^2+2\gamma \lambda +k^2=0$. Since $0<\gamma <|k|$,Forward Euler method has amplification eigenvalues $g_{\pm }=1+h{\lambda }_{\pm }$. Their squared modulus isThusStrict inequality gives numerical damping; at the upper boundary the amplitude is constant despite physical decay, so that boundary is a poor accuracy choice. In the undamped limit $\gamma =0$, $|g{|}^2=1+k^2{h}^2>1$ for every $h>0$: forward Euler method is unstable for any nonzero step on an undamped harmonic oscillator.

Backward Euler method has update ${\mathbf{y}}^{n+1}=(I-hM{)}^{-1}{\mathbf{y}}^n$ and amplification eigenvalues $g_{\pm }=(1-h{\lambda }_{\pm }{)}^{-1}$. Hencefor every $h\ge 0$. It is therefore unconditionally stable for this linear system. Even with $\gamma =0$, it damps every nonzero numerical mode: stability is obtained at the cost of artificial damping and phase error. Large $h$ does not accurately resolve the oscillation period.

An explicit numerical method normally has a simple, inexpensive step and needs no simultaneous solve, but stability can force very small steps, especially for stiff differential equations. An implicit numerical method requires a linear or nonlinear solve, storage and suitable solver tolerances; its larger stability region can permit larger steps when stability is the main restriction. Accuracy, phase fidelity and preservation of physical structure still constrain the step. For this oscillator, backward Euler method is stable but dissipative, while forward Euler method is cheap but unstable without damping. Neither statement makes all implicit methods superior to all explicit ones; method order, problem structure and total cost determine the useful choice.

The line references below use the original listing's line numbers.

For this linear homogeneous problem, the exponential update is exact for the retained modes, apart from roundoff error. It is not a forward Euler method step and does not require the explicit finite-difference diffusion stability restriction. Sampling, periodic modelling and spatial resolution still matter.

Let $U\subset \mathbb{R}\times {\mathbb{R}}^d$ be an open set, $d<\infty$, and $f:U\to {\mathbb{R}}^d$ a continuous function. Through every $(t_0,y_0)\in U$ there is a local solution of a differential equation: a continuously differentiable function $y:I\to {\mathbb{R}}^d$ on an open interval containing $t_0$, satisfying the displayed initial value problem pointwise. No Lipschitz condition is assumed.^[ref]

Continuity of a proposed rate law therefore guarantees at least one trajectory for a short time. For a physical model, this is a first consistency check: can the stated rule and initial condition actually be followed? Prediction also needs uniqueness. The delayed solutions of a square-root differential equation below show how the same rule and initial condition can leave a departure time undetermined. The Picard-Lindelöf theorem supplies a sufficient additional hypothesis for uniqueness; failure of that hypothesis alone does not prove nonuniqueness.

For the initial value problem $y^{\prime }=2\sqrt{|y|}$, $y(0)=0$, a forward solution of a differential equation can stay zero until any $c\ge 0$, then follow $y=(t-c{)}^2$. The delayed solutions of a square-root differential equation verify this freedom exactly.

Local existence gives neither a global solution of a differential equation nor a bound on all future values. The continuous function $f(t,y)=y^2$ already permits finite-time blowup; the finite-dimensional continuation criterion explains the additional issue of extending a trajectory within its domain.

Proof outline: choose a compact set $K=[t_0-a,t_0+a]\times {\overline{B}}_r(y_0)\subset U$. By the extreme value theorem, $M={\max }_K\Vert f\Vert <\infty$. Take $0<h\le a$ with $Mh\le r$. Forward Euler method polygons remain in the closed ball and satisfy $\Vert y_n(t)-y_n(s)\Vert \le M|t-s|$. The Arzelà-Ascoli theorem yields a convergent subsequence with uniform convergence. Uniform continuity of $f$ on $K$ makes their integral equation defect tend uniformly to zero, so the limit satisfiesThe fundamental theorem of calculus gives a continuously differentiable function. Repeat backward and join at $t_0$, where both derivatives equal $f(t_0,y_0)$.^[ref]

The compact set here uses a finite-dimensional vector space; bounded closed balls in an arbitrary Banach space need not be compact. When $M=0$, the constant solution of a differential equation $y=y_0$ works on the entire time interval in $K$. In general the proof selects a subsequence, without asserting that all polygons converge to one unique trajectory.

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.

For a numerical evolution, numerical stability means that perturbations remain controlled over a fixed finite interval by a bound independent of sufficiently small discretization parameters. With perturbations ${\delta }_n$ to starting data and added step defects $d_n$, a typical estimate iswhere $C_T$ is independent of $h$. For a linear homogeneous recurrence, a uniform bound on powers of its update matrix supplies this control. von Neumann stability analysis and the energy method provide complementary tools.

Finite-interval numerical stability differs from absolute stability at a fixed step size. Applying the Euler method to $y^{\prime }=-2y$ multiplies a perturbation by $1-2h$ each step. At $h=0.2$ this is $0.6$, so perturbations shrink; at $h=1.2$ it is $-1.4$, so they alternate and grow. Nevertheless the Euler method converges as $h\to 0$ for this problem: failure at a coarse step does not contradict the convergence theorem. Absolute stability asks whether a selected step preserves boundedness of the Dahlquist test equation over arbitrarily many steps. MIT's analysis uses strict decay, $|1+h\lambda |<1$, rather than the boundary-inclusive boundedness convention.^[ref]