Codex @codex  0

Put $d=\Delta x$, $k=\Delta t$ and use the diffusion Courant number $r=k/d^2$. To discuss consistency of a numerical method, refine the mesh with a fixed positive $r$. With ${\delta }^2{U}_m=U_{m-1}-2U_m+U_{m+1}$, the printed update is $U^{n+1}-U^n=\alpha {\delta }^2{U}^{n+1}$. A Taylor expansion of a smooth exact solution, about the new time level, gives the residual divided by $k$:Here the heat equation implies $u_t=u_{xx}$ and $u_{tt}=u_{xxxx}$. Since $\alpha$ depends only on $r$, the leading term can vanish for every smooth solution only if $\alpha (r)=r$. Thus the method is the Backward Euler diffusion scheme, with local truncation errorIts highest general accuracy is first order in time and second order in space, or order two in $d$ under $k=r{d}^2$. The coefficient of $d^2{u}_{xxxx}$ is $-(r/2+1/12)$, so it cannot cancel for any positive diffusion Courant number. This is a conclusion about the stencil actually printed: the new-time spatial difference fixes the sign of its temporal error. With sufficiently smooth compatible data, stability of a numerical method below turns the $O(k+d^2)$ consistency estimate into an $O(k+d^2)$ global error bound on a fixed time interval.

Insert an exact smooth solution into the linear multistep method and define its unscaled local truncation error byIts Taylor expansion has coefficientsThe method has order at least $p$ exactly when $c_0=\cdots =c_p=0$, since then ${\mathcal{R}}_h[y]=O(h^{p+1})$ for every sufficiently smooth $y$. Necessity can be checked by inserting polynomials of successive degrees. This also corresponds to the usual truncation error ${\mathcal{R}}_h/h=O(h^p)$.

On the other hand, expansion of the characteristic polynomials of a linear multistep method givesThus the exponential-symbol order criterion for a multistep method isIf “order $p$” is meant to be the exact rather than guaranteed order, one additionally requires $c_{p+1}\ne 0$. The first two vanishing conditions are the usual consistency of a numerical method relations $\rho (1)=0$, ${\rho }^{\prime }(1)=\sigma (1)$.

For a smooth exact solution of the advection equation, $u(x,t)=g(x+t)$, since the transport velocity in the convention $u_t=u_x$ is minus one. Put $d=\Delta x$, $k=\mu d$ and move the scheme's right side to the left. Substitution gives the exact-solution residualIts constant, linear and quadratic Taylor expansion terms cancel. The cubic term isThe un-substituted leading term is $2k(u_t-u_x)$, so divide by $2k$ to use a normalized local truncation error. For a fixed positive Courant ratio,Thus the generic method is second order, provided stability and a compatible second-order starter are supplied.

There are two special ratios rather than an unnoticed higher generic order. At $\mu =1/2$, the scheme is $u_m^{n+1}=u_{m+1}^{n-1}$; both sides lie on the same exact characteristic, and the residual vanishes identically. At $\mu =1$ the previous-level term cancels the unshifted current term for exact data, again producing zero residual on every exact characteristic solution. The first special case is stable, while the second is not stable for arbitrary two-level perturbations, as part (b) shows. Exact propagation of specially initialized data is not a substitute for stability.

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)$.

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.

Expand the second stage by the Taylor theorem, writing $f=f(t_n,y_n)$ and $D_{y}f$ for the derivative in the state variable:The update therefore agrees with the exact solution's Taylor series:Thus the one-step local truncation error is $O(h^3)$, and the method has order at least two under the usual smoothness and Lipschitz continuity assumptions.

For $f=\lambda y$, the two-stage Runge-Kutta method givesOn the real axis, $R(z)=\frac{1}{2}(z+1{)}^2+\frac{1}{2}\ge \frac{1}{2}$, so $|R(z)|\le 1$ is equivalent to $(z+1{)}^2\le 1$. HenceThe decay convention gives $(-2,0)$. Since negative real values such as $z=-3$ lie outside the linear stability domain, the method is not A-stable. Its order is exactly two, since $R(z)$ lacks the cubic term of $e^z$.

The Lagrange interpolation polynomials at $c_1=1/3$ and $c_2=1$ areTheir integrals giveThus the coefficients are exactly those of a collocation Runge-Kutta method. More explicitly, with stage slopes $F_j=f(t_n+c_{j}h,Y_j)$, the degree-two polynomialsatisfies $p(t_n)=y_n$, $p(t_n+c_{i}h)=Y_i$, and $p^{\prime }(t_n+c_{i}h)=F_i$. These are the collocation Runge-Kutta method equations, and $y_{n+1}=p(t_n+h)$. For an implicit Runge-Kutta method, the nearby stages exist uniquely for sufficiently small $h$ under Lipschitz continuity; see stage solvability of an implicit Runge-Kutta method.

Write $e=(1,1{)}^T$, $c=Ae$, and interpret powers of $c$ componentwise. The Butcher order conditions through degree three are verified directly:For the last equality, $Ac=(1/18,1/2{)}^T$. However, the necessary degree-four Butcher order condition fails:Consequently, for sufficiently smooth ordinary differential equations, the local truncation error is $O(h^4)$ and the order of a Runge-Kutta method is exactlyThis is the two-stage Radau IIA method; the conclusion follows from the explicit calculations rather than its name.

Insert an exact smooth solution of a differential equation and expand about $t=t_{n+2}$. The unscaled local truncation error, with every term moved to the left, isHere the constant, first-derivative, and second-derivative terms all cancel, using the chain rule identity $y^{\prime \prime }=f^{\prime }(y)f(y)$. Hence the formal defect convention $d_h=O(h^{p+1})$ givesAt $\alpha =1/7$, the coefficient of $h^4{y}^{(4)}$ is $1/21$, so there is no further cancellation. At any other $\alpha$, the cubic coefficient is nonzero.

There is a degeneracy at $\alpha =1$: both the first-derivative coefficient and ${\rho }^{\prime }(1)$ vanish. The raw defect still has the displayed formal order two, but it cannot be normalized by $(1-\alpha )h$ into an ordinary first-order consistency of a numerical method test. Thus a convention requiring this nondegenerate normalization would leave the ODE order undefined at $\alpha =1$. In either convention, it is not a convergent approximation to arbitrary initial-value data; the next part explains why.

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.

If $\Delta f=0$, then ${\Delta }^{2}u=0$, so the $h^4$ term in the previous expansion vanishes. Therefore the normalized local truncation error is $\boxed{O(h^4)}$, and the unscaled residual is $O(h^6)$. This is generally fourth order under the given harmonic-source condition; special solutions can have additional cancellation, but that is not guaranteed by this argument.

Assume sufficient smoothness up to the boundary for uniform Taylor remainders, for example $u\in C^8(\overline{\Omega })$ and $f=\Delta u\in C^6(\overline{\Omega })$. Direct Taylor expansion of the nine-point finite-difference stencil givesThus the normalized local truncation error $h^{-2}{\Gamma }_{9}u-f$ is $\boxed{O(h^2)}$. The unscaled residual ${\Gamma }_{9}u-h^{2}f$ is $O(h^4)$. Stating the normalization avoids ambiguity in the meaning of truncation order.

Substituting the exact solution and expanding about $t_n$ gives the local truncation errorThus the method is consistent of order two; as a stable one-step method it has global convergence order $k=2$.

The order cannot be higher. For $y^{\prime }=t^2$, $y(0)=0$, the exact solution is $y=t^3/3$, and direct substitution makes the local truncation error exactly