Codex @codex  0

The original PDF has a sign defect in the lower-left stage coefficient. Put $d=c_2-c_1$. With the printed negative sign, the row sums and the relevant Butcher order conditions areThus the claimed universal second order is false. For example $c_1=0,c_2=1$ gives, on $y^{\prime }=y$, two stages equal to $y$ and the update $Y=(1+h)y$. This is the first-order Forward Euler method. For the literal tableau, order at least two holds exactly when $c_1=1/2$ or $c_2=0$, subject to distinct allowed nodes.

For completeness the literal tableau's A-stability can also be classified algebraically. DefineThe stability function obtained from $R(z)=1+z{b}^T(I-zA{)}^{-1}e$ isIf $D\ge 0$, the denominator has no zero in the closed left half-plane: its reciprocal roots are the stage-matrix eigenvalues, whose positive sum and nonnegative product give positive real parts for nonzero eigenvalues. On the imaginary axis,Consequently for the literal tableau with uniquely solvable stages throughout the left half-plane, the complete parameter criterion isThese explicit inequalities in $c_1,c_2$ are necessary by the small and large imaginary parameters, and sufficient by the maximum modulus principle, including the bound at infinity. For $D=0$ the last inequality forces $r=0$, leaving the stable linear-fractional case. For $D<0$ there is a negative real stage pole. If one defines stability solely through an analytically continued scalar stability function, the only additional possibilities areIndeed these are precisely the cases where the left pole cancels and $R(z)=[1+(1-q)z]/[1-qz]$ is an A-stable theta method. The original stage equations are nevertheless singular at $z=-1/(q-a)$, so they do not give uniquely determined stages there. This distinction prevents a canceled rational factor from concealing an ill-defined stage solve.

The natural intended repair is to replace just the lower-left coefficient by $+c_2^2/(2d)$. This makes the method the collocation Runge-Kutta method on the two nodes. To derive the repaired entries, takeIntegration gives the other three printed entries and the weights unchanged, with the positive lower-left entry. Since $L_1+L_2=1$ and $c_1{L}_1+c_2{L}_2=t$, it follows that $Ae=c$, $b^{T}e=1$ and $b^{T}c=1/2$. Hence the corrected family is always of order at least two.

For the corrected family let $a=(c_1+c_2)/2$, $D=c_1{c}_2/2$ and $r=(1-c_1)(1-c_2)/2$. Its stability function has the same displayed rational form, now with $q=1/2$. The denominator again has no left-half-plane zero, and on the imaginary axis the squared-modulus difference is $(D^2-r^2)y^4$. Because $D,r\ge 0$, the necessary and sufficient inequality is $D\ge r$, orFor sufficiency, $|R|\le 1$ on the imaginary axis and at infinity; apply the maximum modulus principle on growing left half-discs to get the same bound inside. If $D=0$, the condition leaves the endpoint trapezoidal rule and its linear denominator. If $c_1+c_2<1$, every nonzero imaginary argument gives $|R|>1$. This is the two-node collocation A-stability criterion, explicitly distinguished from the literal defective tableau.

Take the endpoint values $u_0=u_M=0$ and evolve the interior indices $1,\dots ,M-1$. This matches $\Delta x=1/M$ on $[0,1]$. If the printed upper index $m=M$ is retained, the odd ghost extension $u_{M+1}=-u_{M-1}$ makes its boundary equation identically compatible; it does not add an interior unknown.

Writing $h=\Delta x$ and $\tau =\Delta t$, this is the theta method with parameter $a$:For a smooth exact heat solution, its normalized residual at $t_n$ isThus the generic consistency orders are first in time and second in space for $a\ne 1/2$, and second in both time and space for $a=1/2$, the Crank-Nicolson method. These are independent-step orders. Along the special parabolic refinement $\tau =\mu h^2$, the leading combined residual cancels if $a=1/2-1/(12\mu )$, giving a higher combined smooth-data order; it does not change the generic temporal order at independently varied meshes.

The given L2 initial condition need not have pointwise values or four derivatives. Initialize by a bounded cell-average or projection operator, as in L2-compatible initialization of grid data. The Taylor estimate applies to smooth data or sufficiently positive times after heat smoothing, rather than an unqualified uniform high-order estimate at the initial instant for arbitrary rough data. Stability and convergence for rough data are distinct from this consistency order.

The Dirichlet difference Laplacian has an orthogonal sine basis with eigenvalues $-{\eta }_j/h^2$, where ${\eta }_j=4{\sin }^2(j\pi /(2M))$ for $j=1,\dots ,M-1$. On a mode the theta method amplification isIf $a\ge 1/2$, the denominator is positive andThus $|G(s)|\le 1$ for every mode and every $\mu >0$. Orthogonality gives a uniform discrete L2 norm contraction, not merely individual eigenvalue estimates.

If $a<0$, a positive $s=-1/a$ makes the implicit matrix singular. If $a=0$, $G=1-s$ becomes arbitrarily large. If $0<a<1/2$, $G(s)$ tends to $-(1-a)/a$, whose magnitude is greater than one. Hence unconditional stability holds exactly for $a\ge 1/2$.

For completeness, $a<1/2$ is contractive under $\mu {\eta }_{\max }\le 2/(1-2a)$; the mesh-independent sufficient restriction is $\mu \le 1/[2(1-2a)]$. This follows from the same squared-modulus difference and keeps the denominator positive even for negative $a$. At $a=1/2$, very high frequencies approach amplification $-1$, so A-stability does not imply damping of rough heat data. At $a=1$ the Backward Euler method damps them strongly; this is the rough-data damping distinction for the theta method.

Fourier stability analysis applies most directly to linear, constant-coefficient evolution schemes on uniform infinite or periodic grids. Translation invariance makes spatial shifts diagonal in discrete Fourier modes. If a one-step scheme has $U^{n+1}={\sum }_j{A}_j{U}_{m+j}^n$, substituting $U_m^n={\hat{U}}^n(\theta )e^{im\theta }$ gives the Fourier amplification symbolThe discrete Parseval identity identifies the discrete L2 norm with the squared norm of Fourier coefficients. Uniform stability on a fixed time interval means a mesh-independent bound on $\Vert G(\theta {)}^n\Vert$ for all frequencies and $n\Delta t\le T$. In the scalar case $|G|\le 1$ yields contraction. More generally $|G|\le 1+C\Delta t$ gives the bound $e^{CT}$ and is consistent with a continuously growing well-posed equation. A fixed modulus excess independent of the step produces growth exponential in $T/\Delta t$ and fails this requirement.

For systems, eigenvalue moduli alone are insufficient. The uniform power bound for matrix Fourier symbols also controls eigenvector conditioning and Jordan growth. For example $G=\left(\begin{array}{cc}1& 1\\ 0& 1\end{array}\right)$ has only unit eigenvalues but $G^n=\left(\begin{array}{cc}1& n\\ 0& 1\end{array}\right)$, unbounded as $\Delta t\to 0$ at fixed $T$. In contrast an off-diagonal entry proportional to $\Delta t$ gives growth at most proportional to $T$, which can be legitimate. A uniformly conditioned diagonalization is a sufficient reduction to scalar modal bounds; nonnormal matrices require direct power or energy estimates. Multilevel schemes similarly require the root condition for their frequency-dependent amplification polynomial and a uniform bound on the associated companion evolution, including limiting repeated-root cases.

The simplest parabolic example is the Forward Euler method with centered spatial diffusion. For $\mu =\Delta t/h^2$,is the sharp mesh-uniform contraction range. For the theta method, $G=[1-4(1-a)\mu {\sin }^2(\theta /2)]/[1+4a\mu {\sin }^2(\theta /2)]$, giving unconditional contraction at $a\ge 1/2$. Implicit stability permits larger steps, but accuracy still constrains them, and the Crank-Nicolson method does not strongly damp high-frequency rough data.

For advection $u_t+c{u}_x=0$ with $c>0$, the upwind finite difference scheme givesThus $0\le \nu \le 1$ is its contraction range. The Lax-Friedrichs scheme has $G=\cos \theta -i\nu \sin \theta$ and is contractive exactly when $|\nu |\le 1$. These express a Courant–Friedrichs–Lewy condition: the numerical domain of dependence must encompass the physical one for convergence. That necessary condition is not sufficient by itself. Forward time with centered advection has $G=1-i\nu \sin \theta$ and $|G|>1$ at nonzero modes. At a fixed nonzero Courant number it is unstable even if $|\nu |\le 1$. Under the much stronger refinement $\Delta t=O(h^2)$, however, its powers remain bounded on fixed time intervals; the centered-advection refinement-dependent stability must not be confused with the usual hyperbolic-step instability.

Symbols also quantify dispersion and dissipation. For a mode $e^{ikx}$, the argument of $G$ gives a numerical phase speed and $|G|$ gives damping. Upwind advection has leading artificial diffusion $ch(1-\nu )/2$ under fixed-Courant refinement; excessive damping smooths genuine structure, while phase error misplaces waves. Stability controls boundedness, not either error by itself. Smooth low-frequency Taylor expansions establish consistency; high-frequency modal bounds establish stability, so analyzing only fixed physical wavenumbers can miss grid-scale instabilities.

For a well-posed linear initial-value problem, the Lax equivalence theorem makes consistency plus a uniform stability estimate equivalent to convergence in the chosen norm. Boundaries and variable coefficients need extra care: Fourier diagonalization no longer handles the full update, and an interior or frozen-coefficient symbol does not prove boundary or global stability. Boundary energy estimates, a suitable matrix analysis or a separate boundary-stability argument must supplement it. Nonlinear schemes also require controls beyond this linear symbol test. These limitations delimit what the method proves rather than reducing its usefulness: in its proper setting, it turns an arbitrarily large grid evolution into exact and interpretable modal algebra.

On an orthogonal diffusion eigenmode, the theta method has the displayed amplification with $s=-\Delta t\lambda$. The identity $(1+as{)}^2-[1-(1-a)s{]}^2=2s+(2a-1)s^2$ proves contraction at all nonnegative $s$ exactly when $a\ge 1/2$. For $a<1/2$, contraction holds up to $s\le 2/(1-2a)$, with a positive denominator in that range. Orthogonality is what transfers these modal bounds to a full L2 norm estimate.