Codex @codex  0

The discrete L2 inner product is the cell-volume-weighted Hermitian inner product on a uniform grid. Its associated norm is the discrete L2 norm. It provides the finite-dimensional analogue of spatial integration by parts and continuous energy method identities.

For homogeneous Dirichlet endpoints, the centered first-difference matrix is skew-Hermitian and the centered second-difference matrix is negative definite. Thus $\dot{U}=(D_{xx}-\alpha D_x)U$ is contractive in the discrete L2 norm for every real $\alpha$ and every mesh width. The displayed identity follows by summation by parts. Nonnegative off-diagonal entries additionally require $|\alpha |h\le 2$; that maximum-principle condition is not necessary for the L2 energy estimate. numerical consistency and Duhamel's formula then give second-order numerical convergence for smooth solutions and compatible initial approximations.

A three-point implicit advection stencil can match translated data through cubic Taylor expansion while retaining a rational Fourier amplification symbol. For the coefficient family $a=\mu (1+\mu )/6$, $b=(2-\mu )(1+\mu )/3$, $c=(2-\mu )(1-\mu )/6$, $d=(2-\mu )/3$, $e=(1+\mu )/3$, its factor is $G=(d+e{e}^{i\theta })/(a{e}^{-i\theta }+b+c{e}^{i\theta })$. Direct subtraction gives $|D{|}^2-|N{|}^2=\mu (\mu -2)(\mu -1)(\mu +1)(1-\cos \theta {)}^2/9$. The displayed ranges are the contraction ranges in the discrete L2 norm, with a nonvanishing denominator there. For forward time and positive speed parameter, retain the nonnegative ranges. For fixed nonzero nonexceptional Courant number, the normalized truncation error is third order; at $\mu =-1,0,1,2$ the update is the corresponding exact grid translation.

If a spatial method of lines has a Hermitian matrix $H$, its generator $-iH$ is a skew-Hermitian matrix. The matrix exponential $e^{-itH}$ is a unitary matrix, so it preserves the Euclidean norm and any constant-volume discrete L2 norm. Equivalently, differentiating the squared norm gives $2\mathrm{Re}(-i{U}^{\ast }HU)=0$. A real sampled potential preserves the Hermitian property of a symmetric discrete Laplacian. Subsequent time discretization must be assessed separately.

Use the corrected same-node method identified above, with $\delta =1/M$, interior indices $1,\dots ,M-1$ and zero endpoint values. Its spatial matrix $L_{\delta }$ is real symmetric. Summation by parts gives, for a vector extended by zero at the endpoints,Hence $-L_{\delta }$ is positive definite, and ${\sum }_j(v_{j+1}-v_j{)}^2\le 4{\sum }_j{v}_j^2$ bounds its eigenvalues by $4a_{+}/{\delta }^2$. The forward Euler method has amplification matrix $I+k{L}_{\delta }$. Its eigenvalues lie in $[1-4\mu a_{+},1]$, so their moduli are at most one whenSymmetry makes this a contraction in the mesh-weighted discrete L2 norm, uniformly over all grids and all coefficient samples satisfying the bounds. There is also a maximum norm proof: the update weights are $\mu a_{m-1/2}$, $1-\mu (a_{m-1/2}+a_{m+1/2})$, and $\mu a_{m+1/2}$. They are nonnegative and sum to one, giving a contraction after the boundary values are included.

This mesh-independent bound is sharp. For $a\equiv a_{+}$, the spatial eigenvectors are $\sin (j\pi m/M)$ and their amplification factors areIf $\mu >1/(2a_{+})$, sufficiently fine grids have a highest-frequency factor below $-1$, giving exponential growth in the step number.

On one fixed grid, rather than uniformly under refinement, the sharp universal stability of a numerical method bound in the discrete L2 norm is slightly larger:Indeed the quadratic form above is bounded by $a_{+}$ times the constant-coefficient difference form, whose largest eigenvalue is $4{\cos }^2(\pi /(2M))/{\delta }^2$; equality is realized by the constant coefficient. Its limit as $M\to \infty$ is the boxed bound. No such fixed-grid matrix conclusion repairs the ill-defined shifted-boundary update in the literal printed formula.

Put $\nu =k/\delta$ and use the discrete Fourier transform on the infinite grid. A spatial Fourier mode $e^{i(j\theta +\ell \varphi )}$ has time coefficient satisfyingIts amplification polynomial of a multilevel finite difference scheme is ${\zeta }^2-2is\zeta -1$, with rootsFor $|s|<1$ they are distinct and have modulus one. Since $|s|\le 2\nu$, a fixed $0<\nu <1/2$ bounds their separation below by $2\sqrt{1-4{\nu }^2}$. Solving for the two mode amplitudes therefore bounds every power of the two-level companion matrix uniformly in frequency and step number. Parseval identity transfers this bound to the discrete L2 norm.

If $\nu >1/2$, an open set of frequencies around $(\pi /2,\pi /2)$ has $|s|>1$ and a root of modulus greater than one. Fourier data supported there grow exponentially. At $\nu =1/2$, the central frequency has the double root $i$, and the recurrence admits $v_n=(A+Bn)i^n$. Thus the root moduli being one do not suffice: a nontrivial Jordan block makes the two-level powers grow linearly. Although this exact frequency has measure zero on the infinite grid, the matrices depend continuously on frequency. For every step count, a neighborhood has powers nearly as large, so their essential supremum is still unbounded. Wave packets in those neighborhoods rule out a uniform stability of a numerical method estimate in the discrete L2 norm.

For arbitrary perturbations in both starting levels the sharp positive range is thereforeThe zero-step limiting recurrence is bounded but does not advance time. A specially filtered startup can suppress a repeated-root mode; that restriction does not establish stability of the full two-level scheme. This is the two-dimensional leapfrog stability threshold.

Use the discrete inner product $⟨U,V{⟩}_h=h{\sum }_{m=1}^{J-1}{\overline{U}}_m{V}_m$, with endpoints $U_0=U_J=0$. The centered first-difference matrix $D_x$ is skew-Hermitian, while summation by parts givesSince $\alpha$ is real, the centered advection contribution has zero real part. ThereforeThis proves the energy contraction for centered drift-diffusion, and hence stability for every fixed real $\alpha$ and every mesh width, in the discrete L2 norm. The extra inequality $|\alpha |h\le 2$ would make the off-diagonal coefficients nonnegative and support a discrete maximum principle, but is not needed for this energy proof. No time-stepping method is involved here.

For a Fourier mode $U_j^n=G(\theta {)}^n{e}^{ij\theta }$, defineThe amplification factor is $G=N/L$. To derive its stability condition, use the coefficient abbreviations from the preceding solution and set $q=\cos \theta$. The squared moduli areHere $a+c={\mu }^2-\mu +1$, $a-c=2\mu -1$, and $4ac=K(\mu )$. Subtraction and factorization giveThus $|G|\le 1$ for every frequency exactly when $K(\mu )\ge 0$, provided the implicit denominator does not vanish. For the physically nonnegative Courant number, this givesFor a signed convention the full real range is $(-\infty ,-1]\cup [0,1]\cup [2,\infty )$.

There is no hidden denominator failure in these ranges. If $K>0$ and $L=0$, the identity forces $N=0$ and $\cos \theta =1$; but $N(0)=3$. If $K=0$, then $|L|=|N|$, and $N$ could vanish only when its two real coefficients have equal magnitudes. Their sum is three, so equality requires $\mu =1/2$, which is not a root of $K$. Hence $L$ is nonzero for every frequency in the allowed ranges. On the compact frequency interval its inverse is bounded for each fixed allowed $\mu$.

When $K<0$, choose a small nonzero frequency. The denominator remains close to three, while the boxed identity gives $|N|>|L|$. Its fixed-frequency numerical mode grows geometrically, proving instability in the usual fixed-Courant refinement limit. Conversely, for $K\ge 0$, the Plancherel theorem converts the mode bound into a discrete L2 norm contraction on the whole-line or periodic grid. The large-Courant branch is possible because the scheme is implicit; stability alone does not promise small error at a large step.

Take $M$ interior points in each direction and spacing $h=1/(M+1)$. The one-dimensional Dirichlet discrete Laplacian has sine eigenvectors $v_i^{(j)}=\sin (j\pi ih)$, becauseThe five-point Dirichlet Laplacian as a Kronecker sum therefore has product eigenvectors $v_i^{(j)}{v}_r^{(\ell )}$. Adding the reaction term gives modal growth ratesThe semidiscrete matrix is real symmetric, so these modes are orthogonal in the mesh-weighted discrete L2 norm. Non-growth for all times is equivalent to every rate being nonpositive. The largest is ${\lambda }_{11}$, giving the five-point heat-reaction stability thresholdAt equality the first grid mode is neutral; above it that mode grows exponentially. This criterion concerns the spatially semidiscrete, continuous-time evolution. Any subsequent time integrator has an additional stability requirement.

The threshold expands as $2{\pi }^2-{\pi }^4{h}^2/6+O(h^4)$ and is below the continuum threshold for every finite mesh. Thus at the continuum neutral value the grid first mode has a small positive $O(h^2)$ rate. This is compatible with finite-time convergence, although it prevents an all-time contraction on that fixed mesh. In the weaker finite-time mesh-uniform sense, all fixed $\kappa$ have the bound $\Vert U(t){\Vert }_h\le e^{\kappa t}\Vert U(0){\Vert }_h$ because the discrete Laplacian is nonpositive. The boxed range is the natural non-growing stability range, matching the threshold comparison in the preceding part; finite-time stability versus power boundedness is a distinct question.

Fourier stability analysis exploits translation invariance of a linear constant-coefficient finite difference method. On an infinite uniform grid of spacing $h$, a spatial shift acts on $e^{ij\theta }$ by multiplication by $e^{i\theta }$. Hence a scalar two-level scheme with new- and old-time symbols $L(\theta ),N(\theta )$ propagates each Fourier mode by $G(\theta )=N(\theta )/L(\theta )$, provided the new-time operator is invertible. An explicit stencil $U_j^{n+1}={\sum }_{\ell }{a}_{\ell }{U}_{j+\ell }^n$ simply has $G(\theta )={\sum }_{\ell }{a}_{\ell }{e}^{i\ell \theta }$. This converts a spatial evolution into scalar algebra at each frequency.

The Plancherel theorem gives the full norm estimate, rather than just a heuristic about waves. For the discrete Fourier transform $\hat{U}(\theta )={\sum }_j{U}_j{e}^{-ij\theta }$, interpreted in the $L^2$ sense if necessary,Thus a scalar mode bound gives an exact discrete L2 norm bound. Stability on a fixed physical interval requires constants independent of the refining spatial and time meshes:The contraction condition $|G|\le 1$ is sufficient and customary for equations whose exact solutions do not grow. More generally $|G|\le 1+Ck$ suffices, since $(1+Ck{)}^n\le e^{CT}$. At a fixed Courant number, when $G$ is independent of $k$, a mode with modulus strictly greater than one violates the estimate as $n\sim T/k\to \infty$. One must keep this scaling: physical growth such as $e^{\kappa t}$ does not by itself mean finite-time numerical instability. This is finite-time stability versus power boundedness.

For systems, the symbol is an amplification matrix. The requirement is a uniform bound on all its powers, not merely its eigenvalues. A unit-modulus Jordan block $\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has powers $\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$ and is unstable under refinement. A uniformly conditioned eigenbasis or an energy estimate can justify a modal conclusion; eigenvalues alone cannot control a general nonnormal matrix. Multilevel schemes likewise produce an amplification polynomial of a multilevel finite difference scheme, with a uniform root condition and control of starting errors. Repeated unit-modulus roots can give polynomial growth, and roots approaching each other require attention to mesh-uniform constants.

Several elementary examples show how the calculation works. For $u_t+c{u}_x=0$ with $c>0$, forward Euler with centered spatial differences gives $G=1-i\mu \sin \theta$, where $\mu =ck/h$. Since $|G{|}^2=1+{\mu }^2{\sin }^2\theta$, it is unstable at every fixed nonzero Courant number. The upwind finite difference scheme instead gives $G=1-\mu +\mu e^{-i\theta }$ andso its exact contraction range is $0\le \mu \le 1$. The Lax-Wendroff advection scheme has $G=1-i\mu \sin \theta +{\mu }^2(\cos \theta -1)$; expansion gives $|G{|}^2=1-4{\mu }^2(1-{\mu }^2){\sin }^4(\theta /2)$, hence $|\mu |\le 1$. The Courant–Friedrichs–Lewy condition relates these restrictions to the numerical domain of dependence, but is a necessary geometric check rather than a complete stability proof. Implicit schemes can have a global spatial dependence and different Courant ranges, as Question 3 illustrates.

For the heat equation $u_t=\nu u_{xx}$, put $r=\nu k/h^2$. Central spatial differences with forward Euler give $G=1-4r{\sin }^2(\theta /2)$, so $0\le r\le 1/2$. With backward Euler the factor is $[1+4r{\sin }^2(\theta /2){]}^{-1}$, a contraction for every $r\ge 0$. For the Crank-Nicolson diffusion scheme it is $[1-2r{\sin }^2(\theta /2)]/[1+2r{\sin }^2(\theta /2)]$, again a contraction for all $r\ge 0$. At large $r$, however, high-frequency factors approach minus one rather than zero: unconditional stability does not imply strong damping or accurate resolution of every mode.

The connection with convergence is the Lax equivalence theorem for a consistent approximation to a well-posed linear evolution problem. Its forward implication follows directly from the error recursion. If $e^{n+1}=G_h{e}^n+k{\tau }^n$ and $\Vert G_h^j\Vert \le C_T$ for $jk\le T$, thenA convergent initialization and a vanishing normalized local truncation error therefore give convergence. Smooth-data consistency can be extended by density when the approximation and initialization are uniformly bounded. For rough $L^2$ data, cell averages or another L2-compatible initialization of grid data are preferable to undefined point samples. A high formal order does not rescue an unstable amplification operator.

Boundary conditions determine whether Fourier diagonalization is exact. On a periodic grid the update is a circulant matrix, and the finite discrete Fourier transform diagonalizes it at the frequencies $2\pi \ell /N$. Here an interior symbol controls the entire scheme, including the wrapping boundary. On the whole line, the continuous frequency transform serves the same purpose and there is no boundary closure to check.

For a finite interval with homogeneous Dirichlet boundary conditions, centered diffusion is a symmetric tridiagonal problem diagonalized by discrete sine modes. One tests the same diffusion symbol at the admitted sine frequencies, not at arbitrary periodic modes. For homogeneous Neumann boundary conditions, a suitable reflected closure uses discrete cosine modes; the constant mode is present and is neutral for pure diffusion. With the standard reflected endpoint stencil, the natural norm uses half endpoint weights, in which the matrix is self-adjoint. Consequently stable diffusion need not make every solution tend to zero under Neumann conditions. These special orthogonal bases supply an exact substitute for Fourier modes when the operator and closure have the required symmetry.

For transport, prescribe only the incoming physical boundary data. With $u_t+c{u}_x=0$, $c>0$, the inflow is the left endpoint; for $u_t=u_x$ it is the right endpoint. Imposing unrelated data at both ends generally overdetermines the first-order PDE. Outflow rows or ghost points must be chosen compatibly with the interior method. A general boundary closure of a difference scheme destroys the circulant structure, so whole-line von Neumann stability analysis is necessary in many settings but not sufficient for the full initial-boundary problem. A defective boundary row that doubles one boundary unknown each step can create growth even beside a stable interior stencil.

A more refined boundary calculation seeks half-line modes $U_j^n=w^n{\zeta }^j$ with $|\zeta |<1$; these are localized near the boundary. Insert them into both the interior recurrence and boundary equations. A compatible mode with $|w|>1$ proves boundary instability. Absence of one isolated growing mode is not the whole theorem: weakly controlled limiting modes can still destroy a uniform estimate. The Uniform Kreiss--Lopatinskii condition requires a uniformly bounded inverse boundary system on the incoming/decaying mode space. Direct energy estimates are often a clearer alternative to solving that transformed boundary problem.

For example, suppose a first-derivative matrix obeys the discrete summation by parts identity $HD+D^{T}H=e_N{e}_N^T-e_0{e}_0^T$, with positive norm matrix $H$. The semidiscrete advection equation $U_t=-cDU$ then givesWith homogeneous inflow $U_0=0$ and $c>0$, only a nonpositive outflow term remains; prescribed inflow contributes a controlled forcing term. This is the inflow advection energy estimate from summation by parts, and it tests the boundary as part of the operator. For a dissipative semidiscrete matrix $A$, backward Euler also preserves the estimate, because its update givesThis distinguishes spatial energy stability from stability of an arbitrary chosen time integrator.

Finally, nonhomogeneous boundaries can be homogenized by subtracting a suitable lifting, leaving a forced equation whose solution needs a bound in terms of both initial and boundary data. Variable coefficients and nonuniform grids remove translation invariance, so freezing coefficients gives only local diagnostic information. A full energy method, a mesh-uniform operator estimate, or a dedicated boundary-symbol analysis is then needed. Fourier symbols give complete stability tests when an orthogonal mode decomposition diagonalizes the full operator; otherwise the boundaries, matrix powers and norm estimates must be checked as well.

Choose a spatial interval with periodic boundary conditions and a uniform grid, so no separate inflow boundary closure of a difference scheme is required. The discrete Fourier transform is a unitary operator in the discrete L2 norm. Every circulant stencil is diagonal in that transform, and a one-step Fourier multiplier is a contraction exactly when its modulus is at most one for every grid frequency. For a mesh-uniform statement, use all $\theta \in [-\pi ,\pi ]$, the limiting frequency range.

A Fourier mode $e^{im\theta }$ giveswith the coefficients from (a). Squaring moduli and putting $v=\cos \theta$ yieldsThis has the required nonnegative sign for every frequency exactly whenOne must also check that the implicit equation is solvable, rather than divide by a vanishing symbol. HereFor $\mu \ne 1/2$, a zero would have to occur at zero or $\pi$; the latter occurs only at $(1\pm \sqrt{3})/2$, both outside the displayed stable intervals. At $\mu =1/2$, $D=3/4+(1/4)\cos \theta \ge 1/2$. Thus there is no denominator zero in the stable range. For fixed such $\mu$, continuity on the compact frequency interval gives a uniformly bounded inverse, and $|G|\le 1$ proves $\Vert U^n{\Vert }_{2,h}\le \Vert U^0{\Vert }_{2,h}$ for all steps.

Outside these intervals the displayed difference is negative at every nonzero frequency where $D\ne 0$, so $|G|>1$. On periodic refining meshes there are frequencies approaching any selected such value. Their amplification is bounded away from one and is repeated $O(h^{-1})$ times over a fixed physical interval when $k=\mu h$, violating a mesh-uniform stability bound. At the two singular parameters the implicit operator itself fails on meshes containing the Nyquist mode. These observations prove necessity as well as sufficiency.

For the physical convention $k>0$, $h>0$, the result is $0<\mu \le 1$ or $\mu \ge 2$, with zero admitted as the identity limit. The negative intervals describe the algebraic contraction result for signed Courant numbers. In particular, the additional large positive interval must not be discarded by imposing the usual explicit-scheme Courant–Friedrichs–Lewy condition: this scheme is implicit. Stability at large fixed Courant number is not a guarantee of a small error constant.

The eigenvalue stability analysis of a finite difference method reduces a linear evolution discretization to its modal behavior, but it must retain the norm and the dependence on the mesh. For a homogeneous system obtained by the method of lines, $U^{\prime }=A_{h}U$, the solution is $e^{t{A}_h}U(0)$. A one-step time discretization gives $U^{n+1}=G_{h,k}{U}^n$; for a Runge-Kutta method, where its stage matrices are invertible, $G_{h,k}=R(k{A}_h)$. With a multistep method one instead uses an amplification matrix on the augmented vector of several time levels.

The relevant finite-time stability estimate iswith $C_T$ independent of the refining spatial mesh and allowed time steps. Similarly, semidiscrete stability requires $\Vert e^{t{A}_h}{\Vert }_h\le C_T$ for $0\le t\le T$. The norm should represent the continuum problem, such as the discrete L2 norm $\Vert U{\Vert }_{2,h}^2=h{\sum }_j|U_j{|}^2$ in one dimension. A fixed finite matrix having decaying solutions as $t\to \infty$ is not by itself a mesh-uniform result for the partial differential equation.

If $G=X\Lambda X^{-1}$, thenThus eigenvalues of modulus at most one give a uniform bound when the diagonalizing bases have uniformly bounded condition numbers. For a normal matrix in the chosen inner product, the spectral theorem for normal operators supplies an orthonormal basis, so the condition number is one and the norm of $G^n$ is exactly ${\max }_j|{\lambda }_j{|}^n$. For a self-adjoint or skew-adjoint spatial matrix, applying a scalar stability function preserves this favorable modal structure. This is the setting where eigenvalue calculations provide especially clean, reliable step restrictions.

A precise fixed-matrix statement is also useful: a matrix is power bounded for all nonnegative integers $n$ if and only if all its eigenvalues lie in the closed unit disk and every Jordan block at a unit-modulus eigenvalue has size one. Indeed, powers of an interior Jordan block contain a polynomial in $n$ times a decaying geometric factor and are bounded; a larger block on the unit circle produces unbounded polynomial growth. For a family of matrices, the resulting bounds must still be uniform. Moreover, finite-time stability versus power boundedness distinguishes this all-time condition from the physical-time estimate above. For example $G_k=I+kN$, $N^2=0$, has $G_k^n=I+nkN$, bounded for $nk\le T$ despite a nontrivial Jordan block at one. Even the scalar update $1+k$ is stable on fixed intervals, since $(1+k{)}^n\le e^T$, although it represents a growing equation rather than a contraction.

For a constant-coefficient periodic difference operator, Fourier modes diagonalize the spatial matrix. The discrete Fourier transform is a unitary operator, so Parseval's identity identifies its modal maximum with the operator norm. The eigenvalue stability analysis of a finite difference method then becomes von Neumann stability analysis. For the Dirichlet discrete Laplacian, a discrete sine transform plays the same role. These facts explain its advantages: a large matrix calculation becomes a scalar symbol or a known spectral interval, and the extremal eigenvalue directly determines a safe time step.

As a first example, centered space discretization of the heat equation $u_t=u_{xx}$ has eigenvaluesForward Euler method gives $R(k\lambda )=1-4r{\sin }^2(\theta /2)$ with $r=k/h^2$. Requiring all factors in $[-1,1]$ yields $0\le r\le 1/2$ for all periodic meshes. In two dimensions the corresponding uniform bound is $r\le 1/4$. Because these are normal matrices, these scalar bounds prove contraction rather than merely suggest it.

Backward Euler method has $R(z)=1/(1-z)$, so the centered discretization of the heat equation is contractive for every $k>0$. The Crank-Nicolson method has $R(z)=(1+z/2)/(1-z/2)$ and is likewise contractive for these negative real modes. However, its factor tends to $-1$ for very stiff modes, leaving high-frequency oscillations weakly damped; Backward Euler method tends to zero. The eigenvalues therefore expose the presence of stiff differential equations and damping as well as binary stability. Unconditional stability does not remove accuracy requirements.

For an advection equation example, the upwind finite difference scheme for $u_t+c{u}_x=0$, $c>0$, has factor $G=1-\nu +\nu e^{-i\theta }$, with $\nu =ck/h$. Sinceit is a contraction in the periodic discrete L2 norm exactly for $0\le \nu \le 1$. In contrast, Forward Euler method with centered advection has factor $1-i\nu \sin \theta$. Its modulus exceeds one at nonzero active frequencies, and repeated amplification makes it unstable under fixed nonzero Courant refinement. A much smaller scaling $k=O(h^2)$ can bound its finite-time growth, demonstrating again that exact contraction and the most general stability estimate are different claims.

The major limitation is the possible presence of a nonnormal matrix. For example,Both eigenvalues are negative, yet at $t=1$ the second unit vector is amplified by at least $1/(eh)$. This proves that negative spectra do not imply uniform semidiscrete stability. Ill-conditioned eigenvectors, Jordan normal form and transient amplification can invalidate a spectral-only argument. Direct energy estimates or estimates for the resolvent can supply the missing norm control.

Boundary conditions are another limitation. A periodic symbol proves a periodic or whole-line result, not an arbitrary initial-boundary-value problem. Boundary closure of a difference scheme can introduce growing modes or mesh-dependent amplification. A spatial operator with variable coefficients can still be studied through its assembled matrix, but then whether it is a normal matrix, its condition number and the physical inner product must be checked rather than borrowed from a Fourier stability analysis. For the finite element method the natural mass matrix inner product often turns a generalized eigenvalue problem into a self-adjoint one.

For time-dependent or split updates, separate spectra are insufficient to control products. For instance, $G_1=\left(\begin{array}{cc}0& K\\ 0& 0\end{array}\right)$ and $G_2=\left(\begin{array}{cc}0& 0\\ K& 0\end{array}\right)$ each have only zero eigenvalues, while $G_2{G}_1$ has eigenvalue $K^2$. For $K>1$, alternating the two is unstable. A common contractive norm would control the products, but their individual eigenvalues do not.

Finally, stability and consistency of a numerical method have different roles. The Lax equivalence theorem states that, for a well-posed linear initial-value problem and a consistent finite difference method on compatible grid spaces, stability on every fixed time interval is equivalent to convergence of a numerical method as the admissible meshes refine. Bounded grid representation and the correct initial and boundary setup are part of that framework. An eigenvalue estimate can establish its stability hypothesis in the favorable settings described above; it neither proves consistency of a numerical method nor bypasses its uniformity requirements. The method is therefore a powerful modal tool, provided its spectral calculation is connected to a genuine norm estimate.

Let $T$ be the $M\times M$ tridiagonal matrix with diagonal two and adjacent entries minus one, incorporating zero Dirichlet boundary conditions. The actual PDF scheme isAn orthonormal basis of eigenvectors of $T$ has components $\sqrt{2/(M+1)}\sin (jm\pi /(M+1))$, and substitution givesThey form an orthonormal basis by the spectral theorem for real symmetric matrices. Hence $I+\mu T$ is invertible for every $\mu >0$, and the amplification eigenvalues are $g_j=(1+\mu {\lambda }_j{)}^{-1}\in (0,1)$. For the discrete L2 norm $\Vert U{\Vert }_h^2=h{\sum }_m|U_m{|}^2$,This bound has constant one independent of time step and mesh, which is the required stability of a numerical method. Therefore every positive Courant number is stable. If a forcing or local-defect sequence is added, the discrete variation-of-constants formula gives a bound by the initial norm plus the sum of those perturbation norms; this also justifies the global error conclusion in part (a).

Use $N$ distinct periodic grid points $x_m=m/N$, $0\le m<N$, and $h=1/N$. Values at the two ends of the interval represent the same point, so the endpoint must not be stored twice. With indices interpreted modulo $N$, setThis replaces the Dirichlet matrix by the periodic discrete Laplacian. The normalized Fourier modes $N^{-1/2}{e}^{2\pi ijm/N}$ form an orthonormal basis and have amplification factorsEach satisfies $0<g_j\le 1$ for every $\mu >0$. The Parseval identity therefore gives $\Vert U^{n+1}{\Vert }_h\le \Vert U^n{\Vert }_h$ in the discrete L2 norm, uniformly in the grid and time step. Thus all positive Courant numbers remain stable. The constant mode has $g_0=1$, expressing preservation of the spatial mean; all nonconstant modes decay. The formula with modulo indices also handles small grids, where two neighbours can coincide.

Arrange the interior values into a vector $U$. Represent the five-point Dirichlet Laplacian as a Kronecker sum $L_h=h^{-2}(D\otimes I+I\otimes D)$, where $D$ is the tridiagonal matrix with diagonal $-2$ and adjacent entries one. This Kronecker sum is real symmetric; the real sampled potential has a real diagonal matrix $D_V$. ThusThe generator $-i{H}_h$ is a skew-Hermitian matrix. ConsequentlyThe matrix exponential is a unitary matrix, either by differentiating its product with its adjoint or by unitary diagonalization of a normal matrix. For the two-dimensional discrete L2 norm $\Vert U{\Vert }_h^2=h^2{\sum }_{m,n}|U_{m,n}{|}^2$, this givesApplying the same equation to a difference proves stability of a numerical method with a mesh-independent constant one. This is norm conservation of a semidiscrete Schrödinger equation. It concerns continuous time after spatial discretization; an arbitrary subsequent time integrator need not preserve this stability.

The printed coordinates do not discretize the stated square. For $[-1,1{]}^2$ with $M$ interior points in each direction, use $h=2/(M+1)$, $x_m=-1+mh$, $y_n=-1+nh$, and sample $V$ there. The printed $h=1/(M+1)$ and unshifted coordinates instead describe a grid on $[0,1{]}^2$. The matrix proof is valid for either geometry with its corresponding boundary values, so this transcription-independent statement flaw does not alter the stability conclusion.