Past exam of the mathematics course of the University of Cambridge / 2017 / iii / Paper 341 (Numerical Solution of Differential Equations)

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.

For the Dahlquist test equation, put $z=h\lambda$. Eliminating the stages of the Runge-Kutta method gives the stability functionIts poles are $2\pm i\sqrt{2}$, strictly in the right half-plane. To prove A-stability directly, write $z=-r+iy$ with $r\ge 0$. With $D(z)=1-2z/3+z^2/6$ and $N(z)=1+z/3$, elementary expansion givesTherefore $|R(z)|\le 1$ throughout the closed left half-plane; the inequality is strict except at $z=0$. The method is A-stable. Moreover $R(z)\to 0$ as $|z|\to \infty$, so the Radau IIA method is also L-stable.

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 first characteristic polynomials of a linear multistep method isThe root condition for a multistep method requires $|\alpha |\le 1$ and excludes $\alpha =1$, where the root of a polynomial at one is double. At $\alpha =-1$, the two unit-modulus roots of a polynomial are distinct, so zero-stability holds. For $\alpha \ne 1$, the first-order consistency of a numerical method identity is ${\rho }^{\prime }(1)=1-\alpha$, equal to the coefficient of $hf$.

Using the permitted Dahlquist equivalence theorem in this multiderivative multistep method setting givesThis means convergence on every fixed finite time interval where the smooth exact solution of a differential equation exists, with consistent starting values and the nearby branch of each implicit solve. More specifically, the convergence of a zero-stable multiderivative method gives order $p$ when the starting errors are $O(h^p)$ and both $f$ and $f^{\prime }f$ are locally uniformly Lipschitz continuous on the relevant region.

The exclusions are genuine. For $y^{\prime }=0$ and $|\alpha |>1$, an arbitrarily small parasitic starting error is amplified as ${\alpha }^n$. At $\alpha =1$, take $y_0=0$, $y_1=h$, also for $y^{\prime }=0$. The exact solution is zero, but the numerical recurrence gives $y_n=nh$, an error of order one at fixed positive time. This proves failure even though the starting values converge.

On the Dahlquist test equation, $g(y)={\lambda }^{2}y$, so the amplification polynomial of a multistep method becomesBoth roots of a polynomial, including the parasitic one, must satisfy the root condition for a multistep method.

For $\alpha =1/7$, multiply by seven and set $a=7-6z+2z^2$. The polynomial is $a{\xi }^2-8\xi +1$. The complex quadratic Schur criterion says that its roots of a polynomial are strictly in the unit disk ifHere is a direct verification on the entire left half-plane. Put $z=-r+iy$, $r\ge 0$, $t=y^2$, and $a_0=7+6r+2r^2$. Thenand expansion of the squared inequality givesEvery term is nonnegative, and at least one is positive unless $z=0$. Since $q-1>0$, taking square roots proves the strict complex quadratic Schur criterion for $z\ne 0$. At $z=0$, the roots of a polynomial are the simple values $1$ and $1/7$. The leading coefficient cannot vanish because $|a{|}^2\ge 49$. Thus the method at $\alpha =1/7$ is A-stable. Its order three is compatible with the Second Dahlquist barrier: that barrier concerns first-derivative linear multistep methods, whereas this method uses a second derivative.

For $\alpha =1/2$, choose $z=-4$, strictly in the left half-plane. The amplification polynomial of a multistep method reduces toOne root of a polynomial has modulus $(3+\sqrt{17})/4>1$. Thus the method at $\alpha =1/2$ is not A-stable.

The Dirichlet Laplacian eigenfunctions and positive Dirichlet Laplacian eigenvalues on the interval of length two areA Sturm-Liouville eigenfunction expansion gives $u={\sum }_j{q}_j(t){\varphi }_j$, with $q_j^{\prime \prime }+({\lambda }_j-\alpha )q_j=0$. Thus all normal modes oscillate precisely when $\alpha <{\lambda }_1$.

To justify a spatially uniform bound for arbitrary appropriate data, assume the usual wave equation energy class $u(0)\in H_0^1(-1,1)$ and $u_t(0)\in L^2(-1,1)$, or smoother compatible data. The energy method conservesThe sharp Poincaré inequality is $\Vert u{\Vert }_{L^2}^2\le {\lambda }_1^{-1}\Vert u_x{\Vert }_{L^2}^2$. For $0\le \alpha <{\lambda }_1$, this makes $E$ coercive with constant $1-\alpha /{\lambda }_1$; for $\alpha <0$, coercivity is immediate. The one-dimensional Sobolev representative satisfies $|u(x,t)|\le \sqrt{2}\Vert u_x(t){\Vert }_{L^2}$ by the fundamental theorem of calculus and the Cauchy-Schwarz inequality. Hence the conserved energy bounds $u$ uniformly in both variables.

At $\alpha ={\lambda }_1$, the smooth solution $u=t{\varphi }_1$ is unbounded. At $\alpha >{\lambda }_1$, $u=e^{\sqrt{\alpha -{\lambda }_1}t}{\varphi }_1$ is unbounded. ThereforeThe bound depends on the initial data; it is not a single bound for all arbitrarily rescaled solutions. The energy-class hypothesis supplies pointwise meaning and excludes undefined rough-data interpretations.

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.

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.

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

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

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

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

Use the real Hilbert space $V=H_0^1(0,1)$, the zero-boundary Sobolev space, with norm $\Vert v{\Vert }_V=\Vert v^{\prime }{\Vert }_{L^2}$. The Poincaré inequality makes this equivalent to its usual Sobolev space norm. Define the symmetric bilinear form and continuous linear functionalThe Cauchy-Schwarz inequality and the sharp interval Poincaré inequality implyThus $a$ is a bounded bilinear form and a coercive bilinear form, and $\ell$ is bounded because $f$ is square-integrable. The Lax-Milgram theorem supplies a unique weak solution:This is obtained from the differential equation by integration by parts; the test functions have zero boundary trace. Conversely, the weak solution satisfies $y^{\prime \prime }=xy-f$ in the sense of distributions. Since $xy-f\in L^2$, it belongs to $H^2(0,1)\cap H_0^1(0,1)$ and solves the equation almost everywhere.

Equivalently, the Ritz method minimizesover $V$. Indeed $J(y+w)-J(y)=a(w,w)/2$, since the mixed term is $a(y,w)-\ell (w)=0$. This proves existence and uniqueness of the minimizer directly from the weak solution, as well as its equivalence to the variational problem.

Choose a basis ${\varphi }_1,\dots ,{\varphi }_N$ of a conforming finite element space $V_N\subset H_0^1(0,1)$ and write $y_N={\sum }_j{c}_j{\varphi }_j$. Testing the variational problem with each basis function givesThese are also the equations $\partial J(y_N)/\partial c_i=0$ for the Ritz method. The stiffness matrix is symmetric and a positive-definite matrix: for a nonzero coefficient vector $c$, linear independence gives $v_N=\sum c_j{\varphi }_j\ne 0$, andThus there is a unique coefficient vector.

Subtracting the exact and discrete variational problems proves Galerkin orthogonality, $a(y-y_N,v_N)=0$ for every $v_N\in V_N$. In the energy norm $\Vert v{\Vert }_a=\sqrt{a(v,v)}$, the Pythagorean identity yieldsTherefore the Ritz method is the best approximation in the energy norm. The Céa lemma also gives $\Vert y-y_N{\Vert }_V\le (1+{\pi }^{-2}){\inf }_{v_N\in V_N}\Vert y-v_N{\Vert }_V$, so any dense sequence of conforming trial spaces converges.

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.

A Runge-Kutta method advances $y^{\prime }=f(t,y)$ through stagesThe Butcher tableau records $A,b,c$, normally with $c=Ae$. Linear absolute stability tests the amplification of a linear decay mode; nonlinear B-stability compares distances between two solutions of a dissipative vector field. Both describe propagation of perturbations, but they impose different conditions.

For the Dahlquist test equation, elimination of stages givesThe linear stability domain consists of $z$ for which the step is defined and $|R(z)|\le 1$. Poles or singular stage systems must be excluded, even if an output formula has a formal cancellation. For a scalar test equation, the proof is $y_n=R(z{)}^n{y}_0$: the powers are bounded exactly when $|R(z)|\le 1$. For $y^{\prime }=Ly$, the update is $R(hL)$. If $L$ is a normal matrix, unitary diagonalization of a normal matrix makes the operator norm equal to the largest scalar amplification modulus. For a diagonalizable matrix, powers are bounded by the condition number of the eigenvector matrix times the maximal scalar power. These constants must remain uniform for mesh-dependent systems. A defective unit-modulus eigenvalue can instead produce a growing Jordan block; the scalar spectral test alone does not control a general matrix.

An A-stable method accepts the entire closed left half-plane. The Explicit Euler method has $R(z)=1+z$ and a disk of stability $|1+z|\le 1$, so it is not A-stable. No nonconstant polynomial stability function can be A-stable, since it is unbounded on the negative real axis; this excludes all consistent explicit Runge-Kutta methods. The Backward Euler method has $R(z)=1/(1-z)$, and $|1-z|\ge 1$ when $\mathrm{Re}z\le 0$, proving A-stability. It is also L-stable, because $R(z)\to 0$ for large left-half-plane $z$.

The implicit midpoint rule and the trapezoidal rule have the same scalar stability function,The identity $|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z$ proves A-stability. On the imaginary axis their amplification has unit modulus, and along the negative real axis it tends to $-1$. Thus neither is L-stable; very stiff decay can persist as alternating numerical values. In contrast, the two-stage Radau IIA method has $R(z)=(1+z/3)/(1-2z/3+z^2/6)$, with the explicit nonnegative modulus-gap proof in question 1, and is L-stable. This illustrates why L-stability matters beyond A-stability for a stiff differential equation.

For nonlinear problems, assume a real or complex Euclidean inner product and the dissipativity hypothesisfor all relevant states at each common time. Exact solutions satisfyA B-stable Runge-Kutta method reproduces this nonexpansion of distances whenever its stages exist. Put $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$, and $G_i=f(t_n+c_{i}h,Y_i)-f(t_n+c_{i}h,{\tilde{Y}}_i)$. Then $D_i=d+h{\sum }_j{a}_{ij}{G}_j$. Expansion of the output squared norm, followed by substituting $d=D_i-h{\sum }_j{a}_{ij}{G}_j$ in its linear terms, proves the Runge-Kutta contractivity identityIf $b_i\ge 0$ and $M=(m_{ij})$ is a positive semidefinite matrix, the first sum is nonpositive by dissipativity, and the second is nonnegative: expand each coordinate of the stage vectors to express it as a sum of nonnegative quadratic forms. These are exactly the conditions for algebraic stability of a Runge-Kutta method. Hence algebraic stability implies B-stability. Applying B-stability to the linear dissipative problem $y^{\prime }=\lambda y$, or its two-dimensional real form, also proves A-stability whenever stages are well defined.

The Backward Euler method has $b=A=1$, hence $M=1$; the implicit midpoint rule has $b=1$, $A=1/2$, hence $M=0$. Both are algebraically stable and B-stable. For the two-stage Radau IIA method,a positive semidefinite matrix, while both weights are positive. It therefore combines third-order accuracy, B-stability, and L-stability. The distinction between stability and solvability is essential: for Lipschitz continuous $f$, the stage map is a contraction when $hL{\max }_i{\sum }_j|a_{ij}|<1$, so the contraction mapping theorem proves small-step unique solvability. The preceding algebraic stability proof is an estimate for existing stages; it is not a theorem asserting stage existence for arbitrary step sizes and domains.

Linear A-stability is insufficient for nonlinear B-stability. The trapezoidal rule fails B-stability, even though its scalar stability function is identical to that of the B-stable implicit midpoint rule. Its Runge-Kutta method coefficients areFailure of this sufficient algebraic stability test alone would not prove failure of B-stability. An actual counterexample does: take $f(y)=-y^3$ and $h=2$. Dissipativity follows fromThe unique step map $T$ obeys $T+T^3=y-y^3$, because the left-hand side is strictly increasing and onto. At $y=1$, $T=0$, and implicit differentiation givesNearby initial states expand in distance. This proves a true nonlinear failure rather than merely failure of a coefficient criterion.

Finally, stability of a numerical method on finite intervals connects these contractivity properties to accuracy. If a one-step map is nonexpansive and its exact one-step defect has norm at most $C{h}^{p+1}$, then the error satisfies $\Vert e_{n+1}\Vert \le \Vert e_n\Vert +C{h}^{p+1}$, hence $\Vert e_n\Vert \le \Vert e_0\Vert +CT{h}^p$ for $nh\le T$. A Lipschitz step bound $1+Kh$ gives the same order using the discrete Gronwall inequality. For a stiff differential equation, the value of B-stability is that the propagation estimate need not grow with a large negative dissipative rate. The accuracy constant still needs appropriate smoothness and uniform derivative bounds; stability alone does not establish a mesh-uniform error order.

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.