For the homogeneous linear differential equation , a characteristic root is a root of a polynomial . A simple root supplies an exponential solution ; a root of multiplicity supplies for . The sign of the real part of determines spatial growth or decay.
If and is nonnegative on the unit circle, a root of a polynomial there has even multiplicity of a root. Indeed is a real analytic nonnegative function, whose first nonzero term in its Taylor series at a zero must have even order. The map has a simple zero, and multiplication by is nonvanishing, so the order equals the polynomial multiplicity.
The quadratic formula gives
so both roots of a polynomial have real part . Writing , the complex exponential function satisfies , and hence its modulus is . Therefore
Intersecting this condition with the allowed parameter interval gives the answer for either root:
The endpoint is included: the repeated root of a polynomial is , whose exponential has modulus .
Pair the off-circle roots of a polynomial using reciprocal-conjugate root pairing, and split each unit-circle root of a polynomial's even multiplicity equally between the two members of a pair. The fundamental theorem of algebra and the leading coefficient give . None of the selected is zero.
On the unit circle, the identity
turns into
At a point of the unit circle outside the finite set of roots of a polynomial, the product is positive and is nonzero and nonnegative. Its ratio to the product is therefore real and strictly positive. This proves , even though the algebraic expression initially permits a complex constant. Consequently
This proves the Fejér–Riesz theorem for a nonzero trigonometric polynomial of actual order . A positive constant has a constant square-root factor, and the identically zero trigonometric polynomial has ; if for a specified upper order , reduce to the actual order first.
Although the PDF permits assuming even multiplicity, there is a short proof of even multiplicity of unit-circle roots of a nonnegative trigonometric polynomial. The real analytic function cannot have a zero of odd order. Near , has a simple zero, and the nonzero factor leaves the zero order of unchanged. Hence the multiplicity of a root must be even.
For , the coefficient symmetry yields
Thus equals its reversed-conjugate polynomial. If is a root of a polynomial, with , this identity gives , proving
The reciprocal-conjugate root pairing also preserves the multiplicity of a root: reversal and complex conjugation take each factor associated with to the corresponding factor associated with , with the same exponent.
Since , no zero root of a polynomial occurs, so the reciprocal operation is always defined on the roots of a polynomial. Roots off the unit circle are paired on opposite sides of it; roots of a polynomial on the unit circle are fixed by this operation. Pairing alone does not force even multiplicity at those fixed roots of a polynomial; that extra fact uses nonnegativity.
On the Dahlquist test equation, , so the amplification polynomial of a multistep method becomes
Both roots of a polynomial, including the parasitic one, must satisfy the root condition for a multistep method.
For , multiply by seven and set . The polynomial is . The complex quadratic Schur criterion says that its roots of a polynomial are strictly in the unit disk if
Here is a direct verification on the entire left half-plane. Put , , , and . Then
and expansion of the squared inequality gives
Every term is nonnegative, and at least one is positive unless . Since , taking square roots proves the strict complex quadratic Schur criterion for . At , the roots of a polynomial are the simple values and . The leading coefficient cannot vanish because . Thus the method at 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 , choose , strictly in the left half-plane. The amplification polynomial of a multistep method reduces to
One root of a polynomial has modulus . Thus the method at is not A-stable.
The first characteristic polynomials of a linear multistep method is
The root condition for a multistep method requires and excludes , where the root of a polynomial at one is double. At , the two unit-modulus roots of a polynomial are distinct, so zero-stability holds. For , the first-order consistency of a numerical method identity is , equal to the coefficient of .
Using the permitted Dahlquist equivalence theorem in this multiderivative multistep method setting gives
This 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 when the starting errors are and both and are locally uniformly Lipschitz continuous on the relevant region.
The exclusions are genuine. For and , an arbitrarily small parasitic starting error is amplified as . At , take , , also for . The exact solution is zero, but the numerical recurrence gives , an error of order one at fixed positive time. This proves failure even though the starting values converge.
The spatial shift in the PDF also requires a boundary value beyond the supplied right endpoint: if is the last interior index, the newest 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 gives
For the root of a polynomial tending to one as , expansion gives
Consequently 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 and , then
For with ,
Since , the complex quadratic Schur criterion yields both roots of a polynomial in the unit disk, strictly inside except at , where they are and . Uniform stability follows directly from the BDF2 discrete energy identity: the modal recurrence is , whose energy decreases because . Summing the modal energies proves a bound independent of the grid and even of ; no eigenvector-separation assumption is needed. This is the stability of a spatially shifted BDF2 stencil. Thus the literal periodic-grid result is
This 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 gives
The Schur stability criterion gives the root condition for a multistep method for every , with strict inequalities for . 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 and
the corrected Dirichlet discrete Laplacian scheme satisfies
This follows by taking the inner product of with and using summation by parts. It controls both time levels uniformly, regardless of the mesh ratio. Hence the corrected method is stable for
For a fully discrete linear partial differential equation evolution, write in a specified discrete norm. Finite-time stability of a numerical method means
with 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 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 . 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, proves contractivity, while 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 ,
Thus for all frequencies exactly when
Sufficiency follows since ; necessity follows from the highest-frequency mode on even periodic grids. This gives the usual parabolic mesh restriction . For the Backward Euler diffusion scheme,
so every is stable. The Crank--Nicolson method has
again of modulus at most one for all , but stiff modes approach rather than zero. These are PDE manifestations of the A-stability and L-stability distinctions for time integration.
For advection , assume and set . The upwind finite difference scheme is , with
It is contractive when . 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 , the upwind direction must be reversed. By contrast, Forward Euler method time stepping with the centered first finite difference has and for nonzero . Under fixed nonzero refinement it is unstable, since a fixed nontrivial mode grows geometrically over steps. This is not a claim of instability under every imaginable mesh coupling: instead bounds its spurious finite-time growth, since .
The leapfrog advection scheme illustrates a multilevel subtlety. Its amplification polynomial of a multilevel finite difference scheme is
For , both roots of a polynomial have unit modulus and separation at least . A uniformly conditioned eigenbasis of the companion matrix then proves power boundedness of a two-level Fourier scheme for arbitrary history data. At and a grid admitting , the roots coincide on the unit circle, and a Jordan block produces 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 . The discrete summation by parts identity is
Take the inner product of with . The elementary identity gives
This proves unconditional contractivity including the actual boundary treatment. More generally, a dissipative operator gives of operator norm at most one: with , dissipativity implies , and the Cauchy-Schwarz inequality gives . 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 by . If is a normal matrix, the scalar linear stability domain criterion on all controls its powers exactly in the corresponding L2 norm. For the centered Dirichlet discrete Laplacian, its eigenvalues lie between 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 , but the upper-right entry of is . Taking , , this grows like . 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 is nonexpansive in a chosen norm for , the second-order strong stability preserving Runge-Kutta method
is also nonexpansive under the same restriction. Indeed, for two inputs, the triangle inequality gives . Its Taylor expansion is , 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 , iteration gives the discrete Duhamel principle
Thus 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.
Take the current to be dilute, vertically well mixed and Boussinesq approximation, with constant channel width , a horizontal absorbing bed, negligible ambient fluid entrainment, negligible bed friction, and no resuspension. The shallow-water approximation neglects vertical inertia and gives hydrostatic pressure. Particle loss alters reduced gravity but does not alter the leading-order carrier-fluid volume or inertia; otherwise additional mass and momentum terms would be necessary. The current occupies depth , has depth-uniform speed , and has particle volume fraction , with and constant .
Volume conservation of a short channel segment gives . Its suspended particle amount per unit bed area is , and the downward particle deposition flux is . Thus , which becomes on multiplying by .
At vertical coordinate above the bed, excess hydrostatic pressure divided by reference mass density is . Integrating across the depth gives the depth-integrated force flux . The horizontal momentum per unit bed area is and the advected kinematic momentum flux is . The resulting variable-buoyancy shallow water equations are
The derivative of the pressure term acts on both and ; treating reduced gravity as horizontally constant while particles settle would lose part of the force.
For a steady interior, let denote volume flux per unit width. Then the particle balance is a first-order linear differential equation,
Here , is constant, and deposited particles are not returned to the current. In particular the decay length is for , and the cancellation of depth follows from : dilution and transit time compensate in the deposition balance.
The steady momentum equation also gives
The inlet depth selects . Depth and speed at every downstream position are determined by the continuous positive root of a polynomial in
This implicit expression also makes clear why the inlet condition is needed: there are generally two possible positive depths.
The depth branches of a steady depositing gravity current are distinguished by the Froude number . Differentiating the momentum invariant, and using the exponential decay of reduced gravity, gives
For , supercritical flow has and decreasing depth, with and as . Subcritical flow has and increasing depth, with . This indefinitely deepening branch eventually violates the assumed shallow or confined-channel geometry. Hence the depth does not universally decrease: that conclusion needs a supercritical inlet.
For fixed , has its unique minimum at , with . The two positive roots when lie below and above . Since decreases downstream while is constant, an initially strict branch remains separate from the other branch and cannot cross a critical point. A critical inlet with has no finite depth derivative: at the first differentiated equation would equate zero to a positive value. It requires a singular entrance adjustment or physics beyond this steady model. If , both reduced gravity and the continuous depth branch are constant. These steady solutions describe the supplied current interior, not a moving finite-release front.
For two particle populations, write and , allowing different constant particle-density conversion factors. The bidisperse particle-laden gravity current equations are
For the steady case,
Thus the depth is obtained from the same cubic polynomial invariant with the sum of the two exponentially decreasing reduced gravity contributions.
The local deposit consists of the two particle deposition fluxes and . For a positive total deposition rate, the bidisperse gravity-current deposition fraction by particle volume is
If the particles have a common material density, this is also their local mass fraction. If they have common , it can be written in the source's buoyancy notation as
For unequal replace each by . The exponent sign matters: with both populations present and , the faster-settling population has decreasing local deposit fraction downstream, since . Equal positive settling speeds give a constant fraction. A nonsettling population contributes no deposit; if the entire local deposition rate is zero the deposit fraction is undefined.
If “fraction of particles” means number fraction rather than particle-volume or mass fraction, let be a single particle's volume in population . The correct weights are then , so
Settling speeds alone do not determine without particle properties and a drag law. Finally, if the intended deposit is cumulative over the bed from to , rather than local at , integrating particle deposition flux gives the species volumes per unit channel width . The cumulative particle-volume fraction is
For both , its limit at infinity is the incoming particle-volume fraction, whereas its limit as is the local deposition fraction at the inlet. Local sorting and the cumulative inventory are different observables.
For a fixed-step linear multistep method with a nonzero leading coefficient, applied to an ordinary differential equation with a locally Lipschitz continuous vector field, the Dahlquist equivalence theorem states that
Here convergence is uniform on each fixed finite time interval as the step size tends to zero, for every set of starting values tending to the corresponding exact values. Consistency of a numerical method means that the normalized local defect tends to zero; for the characteristic polynomials of a linear multistep method it gives and . Zero-stability is equivalent to the root condition for a multistep method: every root of a polynomial of satisfies , and roots with are simple. Order , together with zero-stability and starting errors , gives global error under the usual smoothness assumptions.
To prove the necessity of the root condition for multistep convergence, use the test problem , . Its discrete error satisfies , where . Set . If , the exact recurrence solution
has starting errors tending to zero at every fixed starting index, but . Thus convergence fails.
If a unit-modulus root has multiplicity , the recurrence instead admits
Again every fixed starting error tends to zero while . These polynomial-times-exponential solutions follow from the repeated factor . Complex modes can be interpreted as a real two-component system, or their real and imaginary parts can be used. Therefore convergence requires precisely the root condition for a multistep method. Repeated roots strictly inside the unit disk are allowed, since their polynomial factors are dominated by exponential decay.