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 three vectors $a,b,c$ in an inner product space, expansion of each squared norm provesTherefore the second-order backward differentiation formula $3U^{n+2}-4U^{n+1}+U^n=2kL{U}^{n+2}$ with a dissipative operator $L$ decreases the discrete energy ${\mathcal{E}}_n=\Vert U^{n+1}{\Vert }^2+\Vert 2U^{n+1}-U^n{\Vert }^2$. This quadratic form is equivalent to the product norm of the two time levels with constants independent of $k$ and $L$. It proves uniform stability even when amplification roots of a polynomial coalesce inside the unit disk, where an eigenvector-separation proof may lose its bound.

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.

Interpolate through the $s+1$ values at $t_n,\dots ,t_{n+s}$, and differentiate the interpolation polynomial at its newest node. In the backward coordinate $\tau =(t-t_{n+s})/h$, the Newton form iswith the $j=0$ factor interpreted as one and $\nabla y_m=y_m-y_{m-1}$. For $j\ge 1$, the derivative of the factor at zero is $(j-1)!/j!=1/j$. Thus the backward differentiation formula isIf $w$ denotes the forward-shift variable, ${\nabla }^j{y}_{n+s}$ corresponds to $w^{s-j}(w-1{)}^j$. Consequently the requested polynomial isThe first form is a polynomial valid also at $w=0$. The second form is convenient for relating it to the truncated series of $\log w$. Differentiation of the degree-$s$ interpolant has local derivative error $O(h^s)$, giving order $s$ for smooth solutions.

A stiff differential equation can contain rapidly decaying modes as well as a slowly changing solution of interest. An explicit scheme may need a very small step solely to keep those already small fast modes from numerical growth. A-stability addresses this restriction through the Dahlquist test equation $y^{\prime }=\lambda y$, whose exact solution decays when $\mathrm{Re}\lambda <0$.

For a one-step method, write $y_{n+1}=R(z)y_n$, $z=h\lambda$. Its linear stability domain consists of the points where the update is defined and $|R(z)|\le 1$. A-stability means that this domain contains the closed left half-plane. This is a scalar linear stability property, not an order condition, an existence theorem for an implicit solve, or a general nonlinear contractivity assertion.

For a rational approximation to the exponential, consistency requires $R(z)=1+z+O(z^2)$, and order $p$ requires $R(z)-e^z=O(z^{p+1})$. Poles must be absent from the relevant half-plane. The Forward Euler method has $R=1+z$ and the disk $|1+z|\le 1$, so negative real modes require $h|\lambda |\le 2$; it is not A-stable. The Backward Euler method has $R=(1-z{)}^{-1}$, and $|1-z|\ge 1$ on the left half-plane, so it is A-stable. The trapezoidal rule and the implicit midpoint rule both have, on the scalar linear problem,They are second-order and A-stable despite being different nonlinear methods. No nonconstant polynomial approximation is A-stable, because its modulus grows without bound on the negative real axis. In particular every nontrivial explicit Runge-Kutta method has a bounded-step stability restriction.

A-stability prevents growth but need not eliminate extremely stiff modes: the trapezoidal factor tends to minus one as $z\to -\infty$. L-stability additionally requires $R(z)\to 0$ within the left half-plane. Backward Euler is L-stable; trapezoidal and implicit midpoint are not. The distinction explains persistent alternating transients in otherwise stable Crank-Nicolson calculations. Accuracy still constrains useful step sizes even for an A-stable method.

A linear multistep method is characterized by polynomials $\rho ,\sigma$, and its scalar modes satisfy $\rho (\zeta )-z\sigma (\zeta )=0$. There is generally no single scalar amplification factor: every root must lie in the closed unit disk and unit roots must be simple. At $z=0$, this is the zero-stability root condition. Consistency plus zero-stability gives convergence by the Dahlquist equivalence theorem, with suitable starting values. A-stability demands the amplification root condition throughout the left half-plane, including the zero-step condition. The Second Dahlquist barrier says that an irreducible A-stable linear multistep method has order at most two, and an explicit linear multistep method cannot be A-stable.

Examples include backward Euler, trapezoidal, and the second-order backward differentiation formulaIts zero-step roots are $1,1/3$. Its unit-circle boundary quotient has real part $(1-\cos \theta {)}^2\ge 0$, and the leading coefficient is nonzero on the left half-plane, giving A-stability by the same continuation argument as in Question 2. Conversely the third-order member of that question fails A-stability. A formal high order without zero-stability, or after silently cancelling a problematic unit factor, does not evade the barrier.

For an implicit Runge-Kutta method, the stages on the test equation satisfy $Y=e{y}_n+zAY$. Whenever these stages are solvable, its stability function isThis rational function is used to test scalar A-stability; the Butcher order conditions separately establish nonlinear order. The Lobatto IIIA method in Question 3 is an A-stable fourth-order example with $R$ the diagonal degree-two Padé approximant. More generally the $s$-stage Gauss--Legendre Runge-Kutta method has order $2s$ and is A-stable, with the diagonal Padé approximation to $e^z$; its stiff-limit factor is nonzero. The $s$-stage Radau IIA method has order $2s-1$ and is L-stable. Thus Runge-Kutta stages permit arbitrarily high A-stable order, unlike irreducible linear multistep formulas.

For a normal matrix, scalar amplification bounds transfer through unitary diagonalization. For a non-normal matrix, eigenvalue information alone does not assert Euclidean contraction: transient growth and eigenvector conditioning can matter. Direct energy method arguments, such as the dissipative Cayley-transform proof in Question 1, use the symmetric part and are more informative in that setting.

For nonlinear dissipative vector fields, the stronger B-stability property asks that distances between two numerical solutions not increase. A standard sufficient condition for Runge-Kutta methods is algebraic stability: $b_i\ge 0$ andis positive semidefinite. It follows from the Runge-Kutta contractivity identity when the implicit stages exist. The three-stage Lobatto IIIA method is A-stable but has the first diagonal entry of that matrix equal to $-1/36$, so scalar A-stability should not be confused with that sufficient nonlinear condition. Stage solvability of an implicit Runge-Kutta method remains a separate issue and computational cost of the implicit equations must be considered. A-stability removes a scalar decay-mode stability restriction; order, stiff damping, nonlinear stability and solvability remain distinct properties.

The spatial shift in the PDF also requires a boundary value beyond the supplied right endpoint: if $m$ is the last interior index, the newest $m+2$ value is outside the interval. Thus the printed Dirichlet method has no specified boundary closure, and no unique stability range can be assigned to it as stated.

For completeness, the literal interior formula has an exact von Neumann stability analysis on an infinite or periodic grid. Substitution of a Fourier mode $U_m^n={\xi }^n{e}^{im\theta }$ givesFor the root of a polynomial tending to one as $\theta \to 0$, expansion givesConsequently $\mu <3$ is unstable: sufficiently small nonzero frequencies have a root of a polynomial of modulus greater than one, and arbitrarily fine periodic grids contain such frequencies.

Conversely, if $\mu \ge 3$ and $s=\cos \theta$, thenFor $A=u+iv$ with $u\ge 3$,Since $|A|\ge 3$, the complex quadratic Schur criterion yields both roots of a polynomial in the unit disk, strictly inside except at $\theta =0$, where they are $1$ and $1/3$. Uniform stability follows directly from the BDF2 discrete energy identity: the modal recurrence is $3v^{n+2}-4v^{n+1}+v^n=-(A-3)v^{n+2}$, whose energy decreases because $\mathrm{Re}(A-3)\ge 0$. Summing the modal energies proves a bound independent of the grid and even of $\mu \ge 3$; no eigenvector-separation assumption is needed. This is the stability of a spatially shifted BDF2 stencil. Thus the literal periodic-grid result isThis cannot establish stability of an unspecified Dirichlet boundary condition closure. It also does not restore consistency of a numerical method.

For the intended unshifted backward differentiation formula, a spatial eigenvalue $-\beta \le 0$ givesThe Schur stability criterion gives the root condition for a multistep method for every $k\beta \ge 0$, with strict inequalities for $\beta >0$. A useful uniform energy method proof, including the repeated interior roots that a bare separation argument would miss, is the BDF2 discrete energy identity. With ${\delta }^{2}U=U^{n+2}-2U^{n+1}+U^n$ andthe corrected Dirichlet discrete Laplacian scheme satisfiesThis follows by taking the inner product of $3U^{n+2}-4U^{n+1}+U^n=2k{L}_h{U}^{n+2}$ with $U^{n+2}$ and using summation by parts. It controls both time levels uniformly, regardless of the mesh ratio. Hence the corrected method is stable for

The original PDF has a spatial shift $u_{m+2}^{n+2}$, so the printed scheme is not the usual second-order BDF discretization. The PDF does have $-2u_m^{n+2}$ on the right; the TeX aid loses that factor two, and that transcription must not be used.

Let $h=\Delta x$, $k=\Delta t$, and expand an exact smooth heat equation solution about $(x_m,t_{n+2})$. Relative to the standard unshifted backward differentiation formula, the printed newest value addsDividing the raw residual by $2k/3$, its expansion isIn particular, for the usual refinement $k=\mu h^2$ with fixed $\mu >0$, the term $3u_x/(\mu h)$ generally diverges. A concrete smooth counterexample is $u=e^{-{\pi }^{2}t/4}\cos (\pi x/2)$, which satisfies the Dirichlet boundary conditions and has $u_x\ne 0$ at $x=1/2$. Thus the printed scheme has no consistent order under fixed-$\mu$ refinement. Its normalized shift error is $O(h/k)$, so one could only investigate a special coupled limit such as $h=o(k)$ after specifying a valid boundary closure; it has no standard joint second-order accuracy.

For the likely intended correction $u_m^{n+2}$, the extra shift terms disappear. The normalized local truncation error isWith a stable, suitably accurate starter, the corrected backward differentiation formula combined with the second-order central difference therefore hasThis is a correction to the printed method, not a claim about it.

For a fully discrete linear partial differential equation evolution, write $U^{n+1}=C_{h,k}{U}^n$ in a specified discrete norm. Finite-time stability of a numerical method meanswith $C_T$ independent of the allowed meshes. For a multistep method, augment the state with its time-history values and include a stable starter. Bounds may grow as $e^{\omega T}$ when the continuous problem has growth; demanding decay for all time would be a stronger assertion. The treatment of initial data, forcing, Dirichlet boundary conditions, and the norm is part of the hypothesis, not a detail supplied by an interior calculation.

For constant coefficients on an infinite or periodic uniform grid, von Neumann stability analysis substitutes the Fourier mode $U_j^n=G(\theta {)}^n{e}^{ij\theta }$. The discrete Fourier transform and the Parseval identity turn a uniform multiplier estimate into a discrete L2 norm estimate. In a one-step scalar scheme, $|G(\theta )|\le 1$ proves contractivity, while $|G(\theta )|\le 1+\omega k$ gives finite-time stability. In a multilevel scheme, every root of a polynomial in its amplification polynomial of a multilevel finite difference scheme matters, as does uniform control of the associated companion matrix.

For the Forward Euler diffusion scheme with $\mu =k/h^2$,Thus $|G|\le 1$ for all frequencies exactly whenSufficiency follows since $G\in [-1,1]$; necessity follows from the highest-frequency mode $\theta =\pi$ on even periodic grids. This gives the usual parabolic mesh restriction $k\le h^2/2$. For the Backward Euler diffusion scheme,so every $\mu \ge 0$ is stable. The Crank-Nicolson method hasagain of modulus at most one for all $\mu \ge 0$, but stiff modes approach $-1$ rather than zero. These are PDE manifestations of the A-stability and L-stability distinctions for time integration.

For advection $u_t+a{u}_x=0$, assume $a>0$ and set $\nu =ak/h$. The upwind finite difference scheme is $U_j^{n+1}=(1-\nu )U_j^n+\nu U_{j-1}^n$, withIt is contractive when $0\le \nu \le 1$. This also has a direct maximum-norm proof: each newest value is a convex combination of two old values. The resulting Courant–Friedrichs–Lewy condition expresses that the numerical domain of dependence covers the physical one. For negative $a$, the upwind direction must be reversed. By contrast, Forward Euler method time stepping with the centered first finite difference has $G=1-i\nu \sin \theta$ and $|G|>1$ for nonzero $\nu \sin \theta$. Under fixed nonzero $\nu$ refinement it is unstable, since a fixed nontrivial mode grows geometrically over $T/k$ steps. This is not a claim of instability under every imaginable mesh coupling: $k=O(h^2)$ instead bounds its spurious finite-time growth, since $\log |G|\le a^2{k}^2/(2h^2)$.

The leapfrog advection scheme illustrates a multilevel subtlety. Its amplification polynomial of a multilevel finite difference scheme isFor $|\nu |<1$, both roots of a polynomial have unit modulus and separation at least $2\sqrt{1-{\nu }^2}$. A uniformly conditioned eigenbasis of the companion matrix then proves power boundedness of a two-level Fourier scheme for arbitrary history data. At $|\nu |=1$ and a grid admitting $\theta =\pi /2$, the roots coincide on the unit circle, and a Jordan block produces $n{\xi }^n$ growth. Thus the endpoint fails the ordinary arbitrary-history root condition for a multistep method; checking only that both moduli equal one misses the instability. Restricting the starter or filtering the parasitic mode is an additional hypothesis.

The energy method handles variable coefficients and finite boundaries for which Fourier stability analysis may be unavailable. For the Backward Euler diffusion scheme with homogeneous Dirichlet boundary conditions, define $\Vert U{\Vert }_h^2=h{\sum }_j|U_j{|}^2$. The discrete summation by parts identity isTake the inner product of $U^{n+1}-U^n=k{L}_h{U}^{n+1}$ with $2U^{n+1}$. The elementary identity $2\mathrm{Re}⟨x-y,x⟩=\Vert x{\Vert }^2-\Vert y{\Vert }^2+\Vert x-y{\Vert }^2$ givesThis proves unconditional contractivity including the actual boundary treatment. More generally, a dissipative operator $L_h$ gives $(I-k{L}_h{)}^{-1}$ of operator norm at most one: with $V=(I-k{L}_h)U$, dissipativity implies $\mathrm{Re}⟨V,U⟩\ge \Vert U{\Vert }^2$, and the Cauchy-Schwarz inequality gives $\Vert U\Vert \le \Vert V\Vert$. In finite dimensions this also proves invertibility. For the second-order backward differentiation formula, the BDF2 discrete energy identity used in question 4 controls both time levels and proves unconditional diffusion stability. A repeated root strictly inside the disk is harmless here; the energy argument supplies the uniform bound without a singular eigenvector formula.

A second technique is eigenvalue stability analysis of a finite difference method. Under the method of lines, a time integrator advances $U^{\prime }=L_{h}U$ by $R(k{L}_h)$. If $L_h$ is a normal matrix, the scalar linear stability domain criterion on all $k{\lambda }_j(L_h)$ controls its powers exactly in the corresponding L2 norm. For the centered Dirichlet discrete Laplacian, its eigenvalues lie between $-4/h^2$ and zero; intersecting this interval with a time integrator's stability interval determines its mesh restriction. A uniform bound on the diagonalization of a matrix is needed for nonnormal diagonalizable systems. Eigenvalues alone can be misleading. For example,has both eigenvalues inside the disk for $0<k<1$, but the upper-right entry of $C_{h,k}^n$ is $nk(1-k{)}^{n-1}/h$. Taking $k=h$, $nk\to T>0$, this grows like $T{e}^{-T}/h$. Hence stable scalar eigenvalues do not give mesh-uniform stability. Boundary closures can introduce precisely the extra growth that an interior Fourier symbol overlooks, so boundary stability of a finite-difference method must be checked separately.

For nonlinear spatial discretizations, a useful replacement for Fourier stability analysis is a convexity argument. If the Forward Euler method map $E_k(U)=U+kF(U)$ is nonexpansive in a chosen norm for $k\le k_{\mathrm{FE}}$, the second-order strong stability preserving Runge-Kutta methodis also nonexpansive under the same restriction. Indeed, for two inputs, the triangle inequality gives $\Vert U_{\mathrm{new}}-{\tilde{U}}_{\mathrm{new}}\Vert \le \frac{1}{2}\Vert U-\tilde{U}\Vert +\frac{1}{2}\Vert E_k(Y)-E_k(\tilde{Y})\Vert \le \Vert U-\tilde{U}\Vert$. Its Taylor expansion is $U+kF(U)+k^2{F}^{\prime }(U)F(U)/2+O(k^3)$, proving order two. Such strong stability preserving Runge-Kutta methods transfer suitable forward-step bounds without relying on a linear spectral calculation.

Finally, consistency of a numerical method explains what stability buys. For a linear Hadamard well-posed problem, the Lax equivalence theorem equates convergence of a consistent discretization with its stability, under the stated approximation-space and norm hypotheses. Directly, if the error satisfies $e^{n+1}=C_{h,k}{e}^n+k{\tau }^n$, iteration gives the discrete Duhamel principleThus vanishing normalized local truncation error and initial error yield convergence. With a nonlinear Lipschitz step estimate, the discrete Gronwall inequality plays the same role. A stable but inconsistent stencil, such as the literal defective formulas in questions 3 and 4 under their usual refinement, is not rescued by any amplification bound. The decisive checks are a mesh-uniform evolution bound, a valid boundary closure, and consistency with the actual PDE.