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]

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]

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.

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 finite difference operator maps grid values to weighted combinations of nearby grid values. Scaled finite differences approximate derivatives; their weights and grid spacing determine the order of a numerical method.

For the advection equation $u_t+c{u}_x=0$ on a uniform grid, the Lax-Wendroff update isIts amplification factor is $G=1-i\nu \sin \theta +{\nu }^2(\cos \theta -1)$, with $|G{|}^2=1-4{\nu }^2(1-{\nu }^2){\sin }^4(\theta /2)$. Thus von Neumann stability analysis on an infinite or periodic grid gives stability for $|\nu |\le 1$. The scheme has second order of a numerical method for sufficiently smooth solutions; finite-interval boundary conditions still need separate treatment.

For a time update with one-step defect ${\delta }_{h,k}$ after exact solution values are inserted, the normalized defect divides by the time step $k$. If ${\delta }_{h,k}=O(k^{p+1})$, the normalized local truncation error is $O(k^p)$. A spatial mesh contributes its own powers: for a diffusion update, ${\mathcal{T}}_{h,k}=O(k+h^2)$ means ${\delta }_{h,k}=O(k^2+k{h}^2)$. Under $k=O(h^2)$ these are respectively second and fourth powers of $h$. Neither normalization changes the actual order of a numerical method; it changes the power used to express its residual.

A numerical flux estimates the rate of transfer across a cell face from nearby states. In the finite volume method update, $F_{j+1/2}^n$ approximates the conservation law flux averaged over the time step; $k{F}_{j+1/2}^n$ is the transferred amount. For a two-state numerical flux $F(L,R)$ approximating a scalar conservation law with conservation law flux $f$, consistency of a numerical method requires $F(U,U)=f(U)$. Thus equal states recover the physical transfer rate. With constant equal-width cells, this also preserves a constant state because its incoming and outgoing numerical fluxes agree.^[ref]

Consistency and shared face flux cancellation answer different questions: matching $f$ at equal states approximates the intended equation, while sharing one oriented face value preserves the total. The example periodic upwind finite volume advection uses $F(L,R)=aL$ for $a>0$, so the arriving state comes from the left. The Godunov numerical flux uses an entropy solution of a local Riemann problem; the Engquist-Osher flux provides a different construction for a scalar conservation law. A monotone conservative scheme uses a numerical flux nondecreasing in its left state and nonincreasing in its right state, together with a time step that makes the complete update monotone. Neither two-state consistency nor conservation by itself supplies that time-step restriction or a general stability or order of a numerical method guarantee.

Suppose an $s$-stage Runge-Kutta method has order $2s-1$, an invertible matrix $A$, and $b^T{A}^{-1}1=1$. Its stability function is $R(z)=1+z{b}^T(I-zA{)}^{-1}1$. The adjugate matrix formula gives a rational function with denominator degree $s$ and numerator degree at most $s$; its zero limit at infinity lowers the latter degree to at most $s-1$. The order of a numerical method gives $R(z)-e^z=O(z^{2s})$, identifying the $[s-1/s]$ Padé approximant. Uniqueness follows because the cross-multiplied difference of two candidates has degree at most $2s-1$ but a zero of order at least $2s$, so is identically zero. The A-stability of near-diagonal exponential Padé approximants then proves A-stability, and the zero limit proves L-stability. These conditions identify the stability function; they do not alone claim uniqueness of all entries of the Runge-Kutta method tableau.

For an initial value problem, convergence of a numerical method means that on each fixed interval $0\le t\le T$, the largest error at mesh points tends to zero as the step size tends to zero and the starting errors vanish. Order of a numerical method describes how rapidly an exact-start local defect decreases; stability of a numerical method controls propagation of that defect and of errors already present. These are related but different properties. A high-order expansion alone cannot ensure that accumulated error remains small.

For a one-step map $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the exact-start defect $d_{n+1}=y(t_{n+1})-{\Phi }_h(t_n,y(t_n))$. Order $p$ means $d_{n+1}=O(h^{p+1})$ for sufficiently smooth solutions. The normalized local truncation error divides this defect by $h$ and is $O(h^p)$; explicitly specifying the normalization prevents an apparent one-order discrepancy. Suppose the map satisfies the local Lipschitz boundFor $e_n=Y_n-y(t_n)$ the error recurrence isIteration, or the discrete Gronwall inequality, yields $\Vert e_n\Vert \le e^{LT}(\Vert e_0\Vert +CT{h}^p)$ for $nh\le T$. Stable propagation converts an order-$p$ local formula into global error $O(h^p)$, provided starting data are equally accurate. Consistency is the weaker requirement that the normalized local defect tend to zero.

For a linear multistep method, zero-stability concerns the recurrence at $h=0$. Its roots must have modulus at most one, and unit-modulus roots must be simple. Roots outside the circle amplify errors exponentially in the number of steps; repeated unit roots introduce polynomial growth. The requirement is stronger than bounded propagation for one fixed number of steps: as $h\to 0$ on a fixed physical interval there are $O(1/h)$ steps. Under the usual Lipschitz assumptions and consistent starting values, the Dahlquist equivalence theorem says consistency plus zero-stability is equivalent to convergence. Its sufficiency combines a bounded discrete recurrence propagator with the sum of local defects, while necessity follows by exciting an unstable homogeneous root or by observing the persistent error of an inconsistent formula.

For example, $Y_{n+2}-3Y_{n+1}+2Y_n=-hf(Y_n)$ satisfies the first-order consistency identities, but its root polynomial is $(\zeta -1)(\zeta -2)$. For $f=0$, a small starting error can excite a $2^n$ mode and destroy convergence as the number of steps increases. In contrast, simple parasitic unit roots such as the root minus one of a two-step centered method need not destroy zero-stability: they remain bounded, though they can create unwanted oscillations and affect practical behavior. Accurate starting procedures and suppression of accumulated parasitic error remain relevant.

Absolute stability asks a different question: for $y^{\prime }=\lambda y$, does the discrete solution decay or stay bounded for a specified finite $z=h\lambda$? The Forward Euler method has $R(z)=1+z$ and, for real $\lambda <0$, requires $h\le 2/|\lambda |$. It converges as $h\to 0$, but a stiff decay mode can impose a severe practical step restriction unrelated to the accuracy needed for slowly varying components. The Backward Euler method has $R(z)=1/(1-z)$, hence is A-stable and L-stable. The trapezoidal rule has $R(z)=(1+z/2)/(1-z/2)$, is A-stable, but tends to minus one on a very stiff negative mode and can retain alternating transients rather than damping them rapidly.

Zero-stability therefore does not imply unrestricted absolute stability; the three-step family analyzed earlier in this paper is a particularly strong illustration. Conversely a stable large step need not be accurate. This distinction matters for a stiff equation. The Second Dahlquist barrier limits A-stable linear multistep methods to order at most two, whereas implicit Runge-Kutta methods can have much higher order and still be A-stable. That barrier concerns the specified method class, not every numerical time integrator.

Finally, the mathematical order depends on smoothness: nonsmooth solutions or forcing can lower observed rates. Roundoff, imperfect implicit solves and inaccurate derivative evaluations introduce additional defects. Smaller steps increase the number of such perturbations; truncation error alone does not describe arbitrarily fine computation. The useful organizing principle is consistency supplies vanishing defects, stability bounds their accumulation, convergence follows, and order quantifies its rate under the stated regularity and starting-error assumptions.

Apply the Runge-Kutta method to the Dahlquist test equation and put $z=h\lambda$. Its stability function isSince $A$ is an invertible matrix, $Q$ has polynomial degree $s$. The adjugate matrix formula gives $\mathrm{deg}P\le s$, but the additional condition givesso actually $\mathrm{deg}P\le s-1$. The specified order of a numerical method implies $R(z)=e^z+O(z^{2s})$ near zero. Therefore the Padé identification of a Radau stability function makes $R$ the $[s-1/s]$ Padé approximant to the exponential function: its numerator and denominator have the required degrees, and its Taylor series agrees through degree $2s-1$.

In the permitted A-stability of near-diagonal exponential Padé approximants, take $m=s-1$, $n=s$. Then $n-2\le m\le n$, so the scheme is A-stable. The same stability function also tends to zero at infinity, giving L-stability.

Expanding the symmetric spatial stencil by the Taylor theorem givesConsistency of a numerical method and cancellation of the $h^2$ term require $\alpha +2\beta +2\gamma =0$, $\beta +4\gamma =1$, and $\beta +16\gamma =0$. Solving this linear system givesThe resulting fourth-order centered second derivative hasThe nonzero next coefficient proves that the highest spatial order of a numerical method is four for this three-parameter stencil.

A constant-step linear multistep method has the formwith ${\alpha }_k\ne 0$. The characteristic polynomials of a linear multistep method encode both accuracy and error propagation. The method has order of a numerical method $p$ when insertion of a sufficiently smooth exact solution leaves an unnormalized local truncation error $O(h^{p+1})$. The exponential-symbol order criterion for a multistep method expresses this as $\rho (e^z)-z\sigma (e^z)=O(z^{p+1})$. In particular, first-order consistency of a numerical method requires $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$.

Accuracy alone is insufficient. For $f=0$, perturbations obey the recurrence $\rho (E)e_n=0$. The root condition for a multistep method requires every root of a polynomial $\rho$ to lie in the closed unit disk, with each unit-circle root simple. This is zero-stability. A root outside the unit disk causes geometric amplification as $n\sim T/h$ increases; a repeated unit-circle root permits polynomial growth in $n$. Either can destroy convergence even for starting perturbations tending to zero.

The Dahlquist equivalence theorem states that, for a fixed-coefficient method and a well-posed ordinary differential equation with a suitably Lipschitz continuous right-hand side, consistency of a numerical method and zero-stability are equivalent to convergence for all convergent starting data. If the starting errors are $O(h^p)$ and the exact solution is sufficiently smooth, the global error is $O(h^p)$ on a fixed time interval. A stability estimate bounds the propagated sum of $O(h^{p+1})$ local defects over $O(1/h)$ steps; the discrete Gronwall inequality controls the dependence of $f$ on the accumulated error. Implicit methods also require the locally consistent, uniquely solvable update branch for sufficiently small $h$.

For example, the explicit Euler method has order one and $\rho (\zeta )=\zeta -1$, satisfying the root condition for a multistep method. The two-step Adams-Bashforth method has order two and $\rho (\zeta )=\zeta (\zeta -1)$, again satisfying it. The trapezoidal rule has order two with one step. The two-step Simpson multistep method, $Y_{n+2}-Y_n=h(f_{n+2}+4f_{n+1}+f_n)/3$, has order four and the simple roots $1,-1$: it is zero-stable, although its parasitic amplification root can spoil long-time decay. These examples distinguish finite-time numerical convergence from absolute stability.

The first Dahlquist barrier gives the maximum possible order of a convergent ordinary real $k$-step method:These bounds assume only first-derivative evaluations, fixed coefficients and zero-stability; the sixth-order multiderivative multistep method in solution is outside that class.

To obtain arbitrarily high convergent order within the ordinary class, increase the number of steps. The Adams-Bashforth method integrates a polynomial interpolation of $k$ past derivative values and has order $k$. The Adams–Moulton method includes the new endpoint and has order $k+1$ for $k$ steps. In both families $\rho (\zeta )={\zeta }^{k-1}(\zeta -1)$, so the root condition for a multistep method holds for every fixed $k$. Starting values can be generated by a Runge-Kutta method of sufficient accuracy. The backward differentiation formula instead differentiates an interpolating polynomial; its order-$k$ members are zero-stable only for $1\le k\le 6$.

Finally, high order and favorable stiff stability are separate design goals. By the Second Dahlquist barrier, an ordinary A-stable linear multistep method has order at most two. One may accept a smaller linear stability domain, use suitable backward differentiation formula members for stiff problems, or switch to high-order implicit Runge-Kutta methods or multiderivative multistep methods. A convergent high-order formula still needs time steps appropriate to its stability properties and the problem's eigenvalues.

Let $g(y)=f_y(y)f(y)$. The chain rule gives $y^{\prime \prime }=g(y)$ along a smooth solution. Expanding every term about $t_n$, the unscaled exact-solution residual of this multiderivative multistep method isFor example, at derivative degree $q\ge 2$ its coefficient isthe coefficients at degrees zero, one and two vanish. Thus the formal order of a numerical method is two for $a\ne 11/5$, and three for $a=11/5$, where the fourth-degree coefficient is $-4/15$. This order statement describes the exact-solution defect; convergence also needs zero-stability.

At $h=0$ the characteristic polynomial isThe root condition for a multistep method holds exactly when $|a|\le 1$ and $a\ne 1$. In particular $a=-1$ has two distinct unit roots and is admissible, whereas $a=1$ has a repeated unit root. Applying the permitted Dahlquist equivalence theorem for this multiderivative setting, with smooth $f$, the nearby implicit solution branch and convergent starting values, yieldsAll those methods have global order two, provided the starting errors are $O(h^2)$. The formally third-order choice $a=11/5$ is not zero-stable and therefore is not convergent as a method for general initial-value problems.

One can see the failure without any nonlinear difficulty. For $y^{\prime }=0$, an error mode $e_n=a^n{e}_0$ grows exponentially if $|a|>1$. At $a=1$, $e_n=C+Dn$; choosing $e_0=0$, $e_1=h$ gives vanishing starting errors but $e_{T/h}=T$. This also explains why a small residual alone does not rescue that degenerate choice: at $a=1$ the first-derivative term disappears and the formula has no zero-stable first-order evolution interpretation.

Expand the exact solution about $t_{n+1}=t$. The unscaled local truncation error isOnly odd powers occur in this centered expansion. Unless $a=1/3$, the first nonzero coefficient is the $h^3$ coefficient, giving order of a numerical method two. For $a=1/3$, that term vanishes but the next coefficient is $-1/90$, giving order of a numerical method four. ThusThese are exact orders for general smooth ordinary differential equations, since the respective first surviving derivatives need not vanish. With starting errors $O(h^p)$, the zero-stability established in part (a) makes the global error $O(h^p)$.

For an evolutionary partial differential equation, stability of a numerical method means that errors in the starting data and forcing stay controlled on each fixed interval $0\le t\le T$, with constants independent of the mesh. Fourier stability analysis is especially effective for a constant-coefficient finite difference method on a uniform infinite or periodic grid: translation invariance makes different Fourier modes evolve independently.

For a scalar one-step scheme on the integer lattice, insert $U_m^n={\hat{U}}^n(\theta )e^{im\theta }$, with $-\pi \le \theta \le \pi$. The Fourier symbol of each shift is $e^{ij\theta }$, and the update reduces toThe amplification factor $G$ must be defined for every relevant frequency; for an implicit scheme this includes checking that its denominator is nonzero. The Parseval identity turns a bound ${\sup }_{\theta }|G(\theta ){|}^n\le C_T$ into a discrete L2 norm bound. Thus $|G|\le 1$ gives contractivity, while the more general bound $|G|\le 1+Ck$ gives $\Vert U^n{\Vert }_h\le e^{CT}\Vert U^0{\Vert }_h$ for $nk\le T$. Conversely, frequencies with amplification uniformly greater than one produce unstable wave packets or periodic Fourier modes under refinement.

For a system, $G(\theta )$ is an amplification matrix; for a multilevel scheme, a companion matrix evolves the vector of time levels. The actual requirement is a uniform bound on matrix powers, not just their spectral radius. A nontrivial Jordan block at a unit-modulus eigenvalue creates polynomial growth in the time index. Even simple eigenvalues inside the unit disk can fail to give a uniform bound if the eigenvector matrices become ill-conditioned as the mesh changes. Uniformly controlled diagonalization of a matrix, or a suitable quadratic energy estimate, resolves this issue. The power boundedness of a two-level Fourier scheme illustrates why root multiplicities and conditioning matter.

For the heat equation, put $\mu =k/h_x^2$. Centered space with the Forward Euler method hasThe von Neumann stability analysis condition $|G|\le 1$ holds exactly for $0\le \mu \le 1/2$ on the full frequency interval. The Backward Euler diffusion scheme instead has $G=(1+4\mu {\sin }^2(\theta /2){)}^{-1}$ and is stable for every $\mu \ge 0$. The Crank-Nicolson diffusion scheme hasIt too is unconditionally stable, but poorly resolved high-frequency modes have $G\approx -1$ for large $\mu$, giving oscillatory numerical transients. A-stability therefore does not guarantee strong damping; L-stability distinguishes the damping of the Backward Euler method.

For the advection equation $u_t+c{u}_x=0$, let $\nu =ck/h_x$. Forward time and centered space giveFor nonzero fixed $\nu$, repeated steps amplify some modes by a fixed factor greater than one, so the method is unstable under the usual refinement $k\propto h_x$. This calculation also shows the refinement qualification: if $k=O(h_x^2)$, its finite-time growth can be bounded by $\exp (C{c}^{2}Tk/h_x^2)$, though the restrictive scaling defeats the usual hyperbolic time-step choice. For $c>0$, the upwind finite difference scheme hasso it is stable for $0\le \nu \le 1$. The Lax-Wendroff advection scheme instead hasgiving $|\nu |\le 1$. These examples separate the effects of the spatial stencil, temporal approximation and Courant number; a higher order of a numerical method by itself does not establish stability.

A periodic boundary condition is ideal for this analysis. The discrete Fourier transform diagonalizes the circulant update, with only the discrete frequencies $2\pi j/J$ needed on a grid of $J$ points. A bound on all $\theta$ is a convenient guarantee over every such grid. On the whole line, the Fourier transform uses a continuous frequency interval and the Parseval identity proves the corresponding square-summable-data estimate.

Homogeneous Dirichlet boundary conditions do not admit arbitrary complex exponential Fourier modes. For the standard centered second difference on $M$ interior points, a discrete sine transform diagonalizes the actual boundary-value matrix:The same heat amplification formulas then apply at these sine frequencies. For example, the exact forward-Euler contractivity limit on that fixed grid is $\mu \le [2{\cos }^2(\pi /(2(M+1))){]}^{-1}$; the mesh-independent sufficient limit $1/2$ follows by including all frequencies. Homogeneous Neumann boundary conditions can similarly permit a discrete cosine transform, provided the endpoint discretization and its weighted inner product are chosen consistently. The constant mode is then present, reflecting conservation of the heat equation's spatial mean. For inhomogeneous boundary conditions, subtract a suitable lifting of the prescribed boundary data and estimate the induced source using the homogeneous evolution's bound and a Duhamel principle estimate. Bounds for the lifting and source must themselves be uniform in the mesh.

General boundary conditions require separate boundary stability of a finite-difference method. For an advection equation, incoming data are prescribed at the inflow boundary, while an outflow closure must respect the outgoing characteristics. An interior periodic Fourier symbol cannot detect a growing mode confined near a boundary. A normal-mode test $U_m^n=z^n{\kappa }^m$ on a half-line seeks modes with $|z|>1$ and $|\kappa |<1$ that satisfy both the interior recurrence and boundary closure. Their existence proves instability, and uniform control requires more than merely excluding isolated growing roots; the Uniform Kreiss--Lopatinskii condition addresses boundary resolvent bounds. Alternatively, a direct energy method using summation by parts can include the boundary terms and establish the required estimate.

Variable coefficients and nonuniform meshes usually destroy exact Fourier transform diagonalization. Frozen-coefficient Fourier stability analysis is then a useful diagnostic, but it is not automatically a proof for the full variable-coefficient boundary problem. The mesh-uniform energy method in question 5 is an example of a direct proof. Also, an L2 norm proof is not automatically a mesh-uniform maximum-norm proof: finite-dimensional norm-equivalence constants may grow with the number of grid points.

Finally, stability measures error propagation; consistency of a numerical method measures the defect of inserting the exact solution. For a well-posed linear initial-value problem in the chosen norm, the Lax equivalence theorem says that a consistent approximation is convergent exactly when it is stable. Boundary discretization and starting data must be included in that assertion. A Fourier calculation proves convergence only after consistency, well-posedness and the actual boundary treatment have also been checked.

Consider an initial value problem $y^{\prime }=f(t,y)$ on $[0,T]$, with a Lipschitz continuous vector field in the solution variable and enough smoothness for the stated error expansions. Convergence of a numerical method meanswith consistent starting values. An order of a numerical method $p$ quantifies this approximation: with sufficiently accurate starting data and the required stability, the global error is $O(h^p)$. Stability of a numerical method controls the amplification and accumulation of perturbations; it is the link between a small local defect and small global error.

For a one-step formula $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the unscaled local truncation errorThe method has local order $p$ if $\Vert d_{n+1}\Vert \le C{h}^{p+1}$ uniformly along sufficiently smooth solutions. If the defect is divided by $h$, its order is instead $O(h^p)$; the convention must be specified. Assume also the perturbation estimateThe global error $e_n=Y_n-y(t_n)$ then satisfiesIterating, or applying the discrete Gronwall inequality, givesThere are $O(1/h)$ steps on a fixed interval, so stable propagation turns $O(h^{p+1})$ local defects into $O(h^p)$ accumulated error. Starting error $O(h^p)$ is needed to retain order $p$. This explains why matching an exact Taylor expansion is insufficient unless errors are controlled during repeated steps.

For a fixed-coefficient linear multistep method, writeThe consistency of a numerical method conditions are $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. The higher order conditions for a linear multistep method arewith right side zero for $q=0$; exact order $p$ means first failure at $q=p+1$. They come from inserting the exact solution and comparing its Taylor expansion.

Errors at $h=0$ follow the homogeneous recurrence determined by $\rho$. Zero-stability is equivalent to the root condition for a multistep method: all roots of a polynomial have modulus at most one, and every unit-modulus root is simple. Roots outside the disk give exponential growth in the number of steps; repeated unit roots give polynomial growth in that number. Repeated interior roots are harmless because their polynomial factors are dominated by decay. Under the usual Lipschitz continuity and initialization assumptions, the Dahlquist equivalence theorem givesFor a zero-stable order-$p$ linear multistep method, local defects and starting errors of the stated orders yield global error $O(h^p)$ by a recurrence bound and the discrete Gronwall inequality. Every required starting value matters.

A concrete failure of “consistency implies convergence” isHere $\rho =(\zeta -1)(\zeta -2)$ and $\sigma =-1$, so $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=-1$: the method is consistent, of order one. But its root $2$ violates zero-stability. For $y^{\prime }=0$ with exact solution zero and starting values $Y_0=0$, $Y_1=h^2$, the computed solution is $Y_n=h^2(2^n-1)$. These starting errors vanish, yet at $n\approx T/h$ the error diverges. Local accuracy cannot compensate for an unstable recurrence.

A different meaning of stability concerns a fixed step on decaying equations. The Dahlquist test equation $y^{\prime }=\lambda y$ introduces $z=h\lambda$. A one-step method gives $Y_{n+1}=R(z)Y_n$ and is absolutely stable when $|R(z)|\le 1$. For a linear multistep method, the corresponding amplification polynomial of a multistep method is $\rho (\zeta )-z\sigma (\zeta )$, and every root must satisfy the bounded root condition, not just the root approximating $e^z$. Strict asymptotic decay requires the roots to have modulus less than one for $\mathrm{Re}z<0$. This distinction matters at reducible methods with an undamped parasitic mode, as in question 1(c).

The Forward Euler method has order one and $R(z)=1+z$, so its linear stability domain is $|1+z|\le 1$. For $y^{\prime }=-\kappa y$, with $\kappa >0$, it requires $0\le h\kappa \le 2$. It converges as $h\to 0$, yet an excessively large step can produce growing oscillations on an exactly decaying problem. The Backward Euler method has the same order, but $R(z)=(1-z{)}^{-1}$ makes it A-stable and L-stable. For a stiff differential equation, fast decay can impose a severe explicit time-step restriction even when accuracy of the slow dynamics would permit a much larger step.

The trapezoidal rule has order two and is A-stable, but is not L-stable, since its stability function tends to $-1$ along the negative real axis. The Second Dahlquist barrier limits an irreducible A-stable linear multistep method to order at most two. This does not limit implicit Runge-Kutta methods: question 4's Lobatto IIIA method is fourth-order and A-stable. Conversely, question 1's fourth-order member is convergent but not A-stable. These examples show that order of a numerical method, zero-stability and A-stability answer different questions.

For nonlinear dissipative problems, B-stability controls distances between numerical solutions; algebraic stability of a Runge-Kutta method is a sufficient coefficient criterion. Question 4 shows that A-stability alone does not imply that coefficient criterion. Variable time steps require additional analysis of step ratios and error propagation, and stiff problems may require uniform-in-stiffness error estimates beyond the fixed-problem convergence result. The useful chain is local accuracy plus the appropriate perturbation bound, giving global accuracy; absolute and nonlinear stability then determine which practical time steps preserve the intended dynamics.

Convergence of a numerical method on a fixed interval $[0,T]$ means that, as the step size tends to zero and the starting approximations tend to the exact initial values,For a vector-valued ordinary differential equation, replace the absolute value by a norm. The interval remains fixed as $h\to 0$; a small error at only one step is insufficient.

The order of a numerical method is measured by its defect on a smooth exact solution. With the unscaled residual convention for a linear multistep method, order at least $p$ means local residual $O(h^{p+1})$; equivalently, the residual divided by $h$ is $O(h^p)$. Order exactly $p$ means the next coefficient does not vanish in general. For a zero-stable method, under the usual smoothness and Lipschitz continuity assumptions and with starting errors $O(h^p)$, this yields global error $O(h^p)$ on the fixed interval. Local order by itself does not imply convergence without stability.

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.