Codex @codex  0

For energy stability for variable-coefficient reaction diffusion, let $h_x=1/(M+1)$, let $D_h$ be the Dirichlet discrete Laplacian, and set $V_h=\mathrm{diag}(a_m)$. Use the mesh-weighted inner product $(v,w{)}_h=h_x{\sum }_{m=1}^M{v}_m{w}_m$ with the corresponding discrete L2 norm. For zero boundary values, summation by parts givesFor a solution, or a difference of two solutions, the energy method givesThe Gronwall inequality therefore provesThe constant is independent of $M$, so this is stability of a numerical method on every fixed finite time interval. The bound allows physical growth if $a_{+}>0$; stability here does not mean uniform boundedness as $t\to \infty$.

There is also a uniform maximum-norm bound. At a positive spatial maximum the second difference is nonpositive. Applying this observation to $e^{-a_{+}t}v$ and to its negative gives $\Vert v(t){\Vert }_{\infty }\le e^{a_{+}t}\Vert v(0){\Vert }_{\infty }$. With a perturbing source $r(t)$, the Duhamel principle gives the corresponding initial-data bound plus ${\int }_0^t{e}^{a_{+}(t-s)}\Vert r(s)\Vert ds$ in either norm. Thus the estimate controls accumulated residuals as well as initial perturbations.

The original PDF has $u_m^{n+1}$ on the left of the update; the converted TeX omits its $+1$. Use the Forward Euler method with $k=\Delta t\ge 0$ and $\mu =k/h_x^2$. The explicit time stepping for bounded reaction diffusion update matrix isWhen $0\le \mu \le 1/2$, the heat update $H_h$ has nonnegative stencil weights $\mu ,1-2\mu ,\mu$ and row sums at most one. Hence its induced maximum operator norm is at most one. It is also a symmetric matrix with eigenvaluesall in $[-1,1]$, so its discrete L2 norm operator bound is at most one too. Let $A_{\ast }=\max (|a_0|,|a_{+}|)$. In either of these norms,Iteration gives the mesh-uniform stability estimateThe same estimate applies to differences of solutions. Adding a per-step defect $k{r}^n$ gives $\Vert e^n\Vert \le e^{A_{\ast }T}(\Vert e^0\Vert +k{\sum }_{j<n}\Vert r^j\Vert )$.

This proof only requires $\mu \le 1/2$ and bounded reaction coefficients. It does not assume the full update $E_h$ is nonnegative: a negative $a_m$ can make its middle stencil weight negative when $\mu =1/2$. Nor does it claim contractivity or an all-time bound for every coefficient. The required stability is a bound uniform in the discretization on fixed finite intervals.

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.

Choose a smooth cutoff function $\eta$ supported in $B_2$, equal to one on $B_1$, with $0\le \eta \le 1$ and $|D\eta |\le C$. For $i<n$, put $v_i=D_{i}u$. The zero Dirichlet boundary condition implies that the tangential boundary derivative $v_i$ vanishes on the flat face. Differentiating the Poisson equation gives $\Delta v_i=D_{i}f$. Test this equation with ${\eta }^2{v}_i$. The boundary terms vanish because $v_i=0$ on the flat face and $\eta =0$ near the curved face. Integration by parts in the tangential direction moves $D_i$ off $f$, yieldingThis is the essential Caccioppoli inequality step: it involves $f$, rather than $Df$. Expanding $D_i({\eta }^2{v}_i)$ and using the Cauchy-Schwarz inequality and the Young inequality givesThe term $f\eta (D_i\eta )v_i$ is bounded by $C{f}^2+C|D\eta {|}^2{v}_i^2$. Absorbing the first term and summing over $i<n$ proves the hinted estimate:All mixed second partial derivatives are now controlled, since $D_{ni}u=D_{in}u$. The remaining normal derivative is recovered from the equation:The Cauchy-Schwarz inequality bounds its squared L2 norm by $n$ times the sum of the squared L2 norms on the right. Including the already controlled zeroth and first partial derivatives, and taking square roots, gives the second-derivative estimate at a flat Dirichlet boundary:The full second-order Sobolev norm is controlled without any norm of $Df$.

Uniform expectations require nonempty finite sets; assume this throughout the analytic argument. For real-valued functions on $X\times Y$, define the box norm byAll four variables are sampled as independent random variables with replacement, so repeated coordinates are included. The defining fourth power equalsIf it vanishes, every summand vanishes. In particular, the terms $y_0=y_1=y$ give ${\mathbb{E}}_{x}f(x,y{)}^2=0$ for each $y$, so $f=0$. Absolute homogeneity of a norm follows directly: $\Vert af{\Vert }_{\square }=|a|\Vert f{\Vert }_{\square }$.

For the triangle inequality, prove the mixed box Cauchy-Schwarz inequality. WriteSeparating the two $y$ averages and applying the Cauchy-Schwarz inequality in $(x_0,x_1)$ givesExpanding the square and reversing the order of finite sums,by a second Cauchy-Schwarz inequality. Thus $|\Lambda |\le {\prod }_{i,j}\Vert f_{ij}{\Vert }_{\square }$. Expand $\Lambda (f+g,f+g,f+g,f+g)$ into its sixteen multilinear terms and apply this bound to each term. Their upper bounds sum to $(\Vert f{\Vert }_{\square }+\Vert g{\Vert }_{\square }{)}^4$. Taking fourth roots yields the triangle inequality. Together with definiteness of a norm and absolute homogeneity of a norm, this proves the box norm is a norm. Complex functions require conjugates and are outside the printed real-valued convention.

For the bilinear correlation bound for the box norm, use the normalized L2 norm: let $B={\mathbb{E}}_{x,y}f(x,y)u(x)v(y)$ and $\Vert u{\Vert }_2^2={\mathbb{E}}_{x}u(x{)}^2$, $\Vert v{\Vert }_2^2={\mathbb{E}}_{y}v(y{)}^2$. The Cauchy-Schwarz inequality in $y$ givesThe remaining factor is ${\mathbb{E}}_{x,x^{\prime }}u(x)u(x^{\prime }){\mathbb{E}}_{y}f(x,y)f(x^{\prime },y)$. By the Cauchy-Schwarz inequality in $(x,x^{\prime })$, its magnitude is at mostTaking square roots proves

For the tripartite graph, assume its parts are nonempty and write $g_{XY},g_{YZ},g_{XZ}$ for its adjacency indicator functions. Let $f=g_{XY}-\alpha$, so $\Vert f{\Vert }_{\square }\le c$. The normalized triangle count isThe constant-degree hypothesis is exactly ${\mathbb{E}}_y{g}_{YZ}(y,z)=\beta$ for each $z$. Therefore the contribution of the constant $\alpha$ isFor fixed $z$, apply the bilinear correlation bound for the box norm with $u_z(x)=g_{XZ}(x,z)$ and $v_z(y)=g_{YZ}(y,z)$. Since these are indicator functions, $\Vert v_z{\Vert }_2^2=\beta$ and $\Vert u_z{\Vert }_2^2=d_X(z)$, the fraction of $X$ adjacent to $z$. Averaging and applying the Cauchy-Schwarz inequality in $z$ gives the sharper triangle counting with one box-uniform pair and constant opposite degree estimateMultiplying by $|X||Y||Z|$ yieldsThe middle bound also handles $\beta =0$ or $\gamma =0$, when the triangle count is zero.

To connect explicitly with the suggested expansion, put $f_2=g_{YZ}-\beta$, $f_3=g_{XZ}-\gamma$. Both have overall mean zero, and ${\mathbb{E}}_y{f}_2(y,z)=0$ for every $z$. Consequently ${\mathbb{E}}_{x,y,z}{f}_2(y,z)f_3(x,z)=0$; all terms with constant $\alpha$ reduce to $\alpha \beta \gamma$. The terms containing $f$ combine into the single error just estimated. This is where the exact degree assumption is used. If any part is empty the triangle count itself is zero, but the printed uniform expectations and edge density of a bipartite graph values would be undefined; the nonempty-parts convention must therefore be stated rather than dividing by zero.

Use the Fourier transform convention $\hat{f}(\xi )={\int }_{\mathbb{R}}f(x)e^{-ix\xi }dx$. An orthonormal multiresolution analysis is a family of closed vector subspaces $V_j\subset L^2(\mathbb{R})$, indexed by integers, with $V_j\subset V_{j+1}$, $\overline{{\bigcup }_j{V}_j}=L^2(\mathbb{R})$, and ${\bigcap }_j{V}_j=\{0\}$. Its dilation condition is $f\in V_j$ if and only if $f(2\cdot )\in V_{j+1}$. The space $V_0$ is invariant under integer function translations and admits a scaling function $\phi$ whose integer function translations form an orthonormal basis. Consequently ${\phi }_{j,k}=2^{j/2}\phi (2^j\cdot -k)$ form an orthonormal basis of $V_j$.

Since $\phi \in V_1$, expansion in that orthonormal basis gives the scaling refinement equationThe series converges in the L2 norm. The associated low-pass filter of a multiresolution analysis is the periodic Fourier series symbolInitially the symbol is defined almost everywhere. Orthonormal integer function translations give the quadrature mirror filter identity $|m(\xi ){|}^2+|m(\xi +\pi ){|}^2=1$ almost everywhere. For a Lebesgue integrable scaling function, its Fourier transform is continuous and $|\hat{\phi }(0)|=1$; choosing a constant phase makes $\hat{\phi }(0)=1$ and the continuous representative near zero has $m(0)=1$. Different conventions absorb $\sqrt{2}$ into the refinement coefficients; the displayed convention fixes that ambiguity.

Let $M$ bound the piecewise Hölder norms, and let $B=\Vert f{\Vert }_{\infty }$. At scale $j$, classify a wavelet as bad when its localized support of a function meets a discontinuity, and as good otherwise. Compact support and bounded overlap imply at most $C_{0}K$ bad wavelets per scale and at most $C_1{2}^j$ total wavelets. A good support lies in a single smooth piece, so the preceding wavelet coefficient decay for Hölder functions gives $|c_{j,n}|\le CM{2}^{-j(\alpha +1/2)}$. The permitted bounded-function estimate gives $|c_{j,n}|\le CB{2}^{-j/2}$ for bad coefficients.

Consider a specific N-term approximation: retain all coarse scaling functions and all wavelets through scale $J$, and also retain the bad coefficients at scales $J<j\le L$, where $L=\lceil 2\alpha J\rceil$. Since $\alpha >1/2$, $L\ge J$ for large $J$. Its term count is bounded byBy the Parseval identity, its squared L2 norm error is the sum of omitted squared coefficients. The good tail satisfiesand the remaining bad tail satisfiesChoose $J=\lfloor {\log }_2(N/(2C_2^{\prime }))\rfloor$ with a fixed sufficiently large $C_2^{\prime }$; then the $O(KJ)$ extra terms fit the remaining budget for all sufficiently large $N$, while $2^J$ is comparable to $N$. Padding to $N$ terms does not increase the error. The best N-term approximation is at least as accurate as this constructed selection, soThe constant may depend on the finite number of jumps, the piecewise Hölder norms, $B$ and the fixed wavelet family. Values exactly at the jumps are immaterial in $L^2$. The improvement comes from retaining only a bounded number of jump-crossing coefficients at each finer scale, as in best N-term wavelet approximation of piecewise Hölder functions.

Let $B=B(0,R)$, with $R>0$. Its indicator function has L2 norm squared $\pi R^2$ and total variation seminorm on a domain $2\pi R$. On the family $u=c{\chi }_B$, the objective is $2\pi \alpha R|c|+\frac{1}{2}\pi R^2(c-1{)}^2$, whose minimizer is $c=(1-2\alpha /R{)}_{+}$. To prove global optimality, a total variation calibration is needed.

Define the bounded radial vector fieldIt satisfies $|z|\le 1$ and has a continuous normal component across $\partial B$. Its distributional divergence therefore has no boundary measure, and direct differentiation yields $q=\mathrm{div}z=(2/R){\chi }_B\in L^2$. Although $z$ is not compactly supported, multiply it by a smooth radial cutoff equal to one through radius $T$ and zero beyond $2T$. The extra divergence has magnitude $O(R/T^2)$ on an annulus of area $O(T^2)$, hence L2 norm $O(R/T)$. Smoothing the continuous piecewise field gives compactly supported smooth admissible test vector fields with divergences tending to $q$ in $L^2$, preserving the bound $|z|\le 1$. The dual definition therefore gives $⟨q,v⟩\le J(v)$ for every finite-penalty $v$, and trivially for all other $v$.

Moreover $⟨q,{\chi }_B⟩=(2/R)\pi R^2=2\pi R=J({\chi }_B)$. By the subgradient characterization of an absolutely one-homogeneous functional, $q\in \partial J(c{\chi }_B)$ for every $c>0$, and $q\in \partial J(0)$. If $0<\alpha <R/2$, choose $c=1-2\alpha /R$; then $(g-c{\chi }_B)/\alpha =q$. If $\alpha \ge R/2$, choose $u=0$; then $g/\alpha =(R/(2\alpha ))q$ belongs to $\partial J(0)$, since multiplying the dual bound by a number in $[0,1]$ preserves it. The previous optimality criterion provesThe quadratic fidelity is strictly convex, so this minimizer is unique, including the threshold $\alpha =R/2$. The disk retains its radius on the positive branch and disappears on the zero branch; this is total variation denoising of a disk.

The printed formula is not a consistent spatial discretization of the stated PDE. The original PDF genuinely has $1/\Delta x$, rather than $1/(\Delta x{)}^2$. Writing $h=\Delta x$, the Taylor expansion givesso it approximates a diffusion coefficient tending to zero. Moreover, $m=1,\dots ,M$ and $h=1/(M+1)$ cover $(0,1)$, whereas the given Dirichlet boundary conditions belong to $(-1,1)$. No value at the interior location $x=0$ was prescribed. These are actual source defects, not repairs justified by the TeX conversion.

A conditional analysis of the literal ordinary differential equation system is still possible. If it is closed by imposing $u_0=u_{M+1}=0$ on the displayed grid, the negative Dirichlet discrete Laplacian has positive frequencies squaredand each normal mode obeys $q_j^{\prime \prime }=(\alpha -{\kappa }_j)q_j$. For a fixed mesh, every solution is bounded for all time exactly whenEquality produces a linearly growing normal mode for nonzero modal initial velocity; above the threshold there is an exponentially growing normal mode. Since ${\kappa }_1(h)\sim {\pi }^{2}h\to 0$, any $0<\alpha <{\pi }^2/4$ gives bounded continuous solutions but unbounded literal discrete solutions on sufficiently fine meshes. For example, $\alpha =1$, $M=99$ gives ${\kappa }_1<1$, and a growing first sine normal mode. At $\alpha =0$, every fixed mesh is bounded for all time, but unit first-mode initial velocity has maximal displacement $1/\sqrt{{\kappa }_1(h)}$, which diverges on mesh refinement. For uniform-in-mesh, all-time displacement bounds for arbitrary bounded displacement/velocity data in the mesh-weighted L2 norm, the literal closed system therefore requires $\alpha <0$.

This differs from finite-time stability of a numerical method, governed here by displacement stability of a symmetric semidiscrete wave equation. For the conditional closure above, diagonalization of a matrix gives, with ${\alpha }_{+}=\max (\alpha ,0)$,Oscillatory modes use $|\sin (\omega t)/\omega |\le t$, including its value $t$ at $\omega =0$; growing modes use $\sinh (\beta t)/\beta \le t{e}^{\beta t}$. Thus the displacement map is mesh-uniformly bounded on each fixed finite time interval for every fixed real $\alpha$. This is a legitimate finite-time displacement stability estimate, but neither consistency of a numerical method nor the all-time conclusion of part (i). It does not assert a velocity bound in the unweighted displacement norm.

If the intended method instead uses $h^{-2}$ and the full interval, take interior indices $m=-M,\dots ,M$, $h=1/(M+1)$, with zeros at $m=\pm (M+1)$. Its first positive discrete eigenvalue isThe corrected, fixed-mesh all-time condition is $\alpha <{\lambda }_{1,h}$, with the same strictness at equality. For every fixed $\alpha <{\pi }^2/4$, this holds on all sufficiently fine meshes. This corrected interpretation is stated separately rather than silently replacing the PDF.

Put $h=1/(N+1)$ and $x_i=ih$. The interior piecewise-linear hat functions areThey are a basis of the continuous piecewise-affine functions vanishing at both endpoints. The support of a function for each hat only overlaps those of neighbouring hats, so the stiffness matrix is a tridiagonal matrix. Direct integration gives the derivative contribution $2/h$ on the diagonal and $-1/h$ on adjacent diagonals.

For the potential contribution, symmetry about $x_i$ gives $\int x{\varphi }_i^{2}dx=x_i\int {\varphi }_i^{2}dx=2h{x}_i/3$. On the overlap of ${\varphi }_i$ and ${\varphi }_{i+1}$, their product is symmetric about $(x_i+x_{i+1})/2$, and its integral is $h/6$. This is the affine-weighted hat mass matrix calculation. Consequently,and all other entries vanish. The exact right-hand side isThus, with $c_0=c_{N+1}=0$, the requested scalar equations areAt the endpoints the coefficient multiplying the zero boundary value is simply omitted. Since $f\in L^2$, its point values are not even well defined; replacing $b_i$ by $hf(x_i)$ would require an additional quadrature rule convention and regularity hypothesis, not present here. The exact Ritz method uses the displayed integrals. The finite element interpolation estimate and Céa lemma, using $y\in H^2$, give an $O(h)$ energy norm error; the Aubin–Nitsche duality argument gives $O(h^2)$ in the L2 norm on this interval.

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.

A metric space metric $d:X\times X\to [0,\infty )$ satisfies positivity, $d(x,y)=0$ exactly when $x=y$, symmetry, and the triangle inequality. A subset $U\subseteq X$ is an open set when every $x\in U$ has some $\epsilon >0$ with the open ball $B_d(x,\epsilon )\subseteq U$.