Codex @codex  0

Yes: this is an A-stable third-order method. Its second derivative is essential; the Second Dahlquist barrier applies to ordinary first-derivative linear multistep methods, not this multiderivative multistep method.

Apply the Dahlquist test equation $y^{\prime }=\lambda y$ and put $z=h\lambda$. The amplification polynomial isWe verify the root condition for a multistep method throughout $\mathrm{Re}z\le 0$. For a quadratic $a{\zeta }^2+b\zeta +c$ with $|a|>|c|$, the complex quadratic Schur criterion isfor both roots to be strictly inside the unit disk. Here is a brief justification of the criterion rather than just a root plot. Set $P^{\#}(\zeta )={\zeta }^2\overline{P(1/\overline{\zeta })}$. The polynomial $\overline{a}P-c{P}^{\#}$ is $\zeta (D\zeta +E)$, where $D=|a{|}^2-|c{|}^2$ and $E=\overline{a}b-c\overline{b}$. Its two roots are inside precisely when $|E|<D$. On the unit circle $|P^{\#}|=|P|$, and $|c|<|a|$, so the Rouche theorem equates its interior root count with that of $P$. A boundary zero of $P$ would also be a boundary zero of the transformed polynomial, which the strict inequality excludes. This proves the strict criterion; the non-strict form follows by continuity, with any unit root simple because $|c/a|<1$.

For the present polynomial it is enough to show $8|Q-1|<|Q{|}^2-1$, except at $z=0$. Write $z=-r+iy$, $r\ge 0$, and put $k=r^2+3r$, $s=y^2$. ThenDirect expansion gives the nonnegative factorizationEvery term is nonnegative, and the sum is strictly positive unless $r=y=0$. Also $|Q{|}^2=X+1>1$, so the criterion applies. For every nonzero $z$ in the closed left half-plane both amplification roots are strictly inside the unit disk. At $z=0$ they are $1$ and $1/7$, and the unit root is simple. ThereforeThis proof establishes stability of both the principal and parasitic modes, including on the imaginary axis.

For $y^{\prime }=f(t,y)$ and step $h$, a collocation Runge-Kutta method constructs a polynomial $P$ of degree at most $s$ on $[t_n,t_n+h]$ satisfyingThe step value is $y_{n+1}=P(t_n+h)$. Thus the differential equation holds exactly at the prescribed collocation nodes, while the initial value fixes the integration constant. For a smooth vector field and small enough $h$, the implicit stage equations have the local solution branch continuing the zero-step initial state; arbitrary large steps need not have a unique nonlinear solution.

Let ${\ell }_j(\tau )={\prod }_{k\ne j}(\tau -c_k)/(c_j-c_k)$ be the Lagrange interpolation basis. Set $Y_i=P(t_n+c_{i}h)$ and $F_j=f(t_n+c_{j}h,Y_j)$. Since $P^{\prime }$ is a polynomial of degree at most $s-1$, the collocation conditions giveIntegrate from zero to $c_i$ and to one to obtainThese are the stage and update equations of a Runge-Kutta method with nodes $c_i$. Conversely, stages satisfying these equations reconstruct $P$ by integrating the displayed interpolation polynomial, so the equivalence works in both directions. Because ${\sum }_j{\ell }_j=1$, they also satisfy the collocation tableau row-sum identity ${\sum }_j{a}_{ij}=c_i$.

The collocation order theorem identifies the classical order with the order of its interpolatory quadrature. If ${\sum }_i{b}_{i}g(c_i)$ integrates every polynomial of degree at most $p-1$ exactly, but fails for some polynomial of degree $p$, then the collocation method has order $p$ for sufficiently smooth differential equations and the local implicit stage branch. With $s$ nodes, $s\le p\le 2s$.

Equivalently, put ${\pi }_s(\tau )={\prod }_i(\tau -c_i)$. If its first $k$ moments vanish,and the next does not, then $p=s+k$. Polynomial division explains the quadrature part: any polynomial through degree $s+k-1$ is a remainder of degree below $s$ plus ${\pi }_s$ times a polynomial of degree below $k$. The remainder is integrated exactly by interpolation, and the remaining integral vanishes. Conversely a nonzero next moment is an explicit failed quadrature test. The differential-equation order assertion is the collocation theorem being stated here. In particular Gauss collocation methods achieve $2s$, Radau IIA methods achieve $2s-1$, and endpoint-including Lobatto IIIA methods achieve $2s-2$ for $s\ge 2$.

A-stability removes the scalar linear decay restriction on the step size; it does not replace accuracy, nonlinear solvability or a suitable system norm. Its importance is clearest for a stiff differential equation, where rapidly decaying modes coexist with slowly varying quantities of interest. On the Dahlquist test equation $y^{\prime }=\lambda y$, the exact solution decays for $\mathrm{Re}\lambda <0$. A stable time integrator should not turn that decay into numerical growth.

For a one-step method write $y_{n+1}=R(z)y_n$, where $z=h\lambda$ and $R$ is its stability function. Its linear stability domain consists of parameters with $|R(z)|\le 1$ and a well-defined update. The method is A-stable when this domain contains the closed left half-plane. A consistent order-$p$ one-step method has $R(z)=e^z+O(z^{p+1})$ near zero, but that local agreement alone says nothing about stability at large negative $z$.

For a Runge-Kutta method with stage matrix $A$, weights $b$ and $1=(1,\dots ,1{)}^T$, eliminating the stages on the test equation givesAn explicit finite-stage method has a polynomial $R$. A nonconstant polynomial is unbounded along the negative real axis, so no consistent explicit Runge-Kutta method is A-stable. The explicit Euler method, with $R(z)=1+z$, is stable in the disk $|1+z|\le 1$ and on the negative real axis only for $-2\le z\le 0$. A stiff mode with a large negative eigenvalue therefore forces an unnecessarily small step even after its exact contribution has almost disappeared.

The Backward Euler method has $R(z)=1/(1-z)$. For $\mathrm{Re}z\le 0$, $|1-z{|}^2=1-2\mathrm{Re}z+|z{|}^2\ge 1$, so it is A-stable. The trapezoidal rule hasand is also A-stable. The implicit midpoint method has the same scalar stability function, although its nonlinear update differs from the trapezoidal rule. The first method has order one and the latter two order two. Their implicit solves trade greater work per step for freedom from the scalar stiff stability restriction.

Damping of very stiff modes separates A-stability from L-stability. For a one-step method, the latter additionally requires $R(z)\to 0$ as $z$ tends to infinity in the left half-plane. Backward Euler method satisfies this. The trapezoidal and midpoint rules instead have $R(z)\to -1$ along the negative real axis: a very stiff mode can alternate in sign with little damping even though its amplitude remains bounded. L-stability is often preferable for diffusion and rapid relaxation, while retaining oscillatory modes can be desirable in conservative dynamics.

For a linear multistep method, writeIts test-equation amplification roots satisfy $\rho (\zeta )-z\sigma (\zeta )=0$. Absolute stability requires every root to lie in the unit disk or on its boundary with all unit-modulus roots simple; all parasitic modes matter. At $z=0$ this is the zero-stability root condition. Consistency requires $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$, and the Dahlquist equivalence theorem combines consistency and zero-stability to characterize convergence. Thus a claim of stability obtained by canceling a common factor is not a proof about the original recurrence and its arbitrary consistent starting data.

Explicit linear multistep methods cannot be A-stable either. If ${\beta }_k=0$ and ${\alpha }_k\ne 0$, some lower coefficient of $\rho -z\sigma$ becomes unbounded as $z\to -\infty$, while the leading coefficient stays fixed. If all roots remained in the unit disk, every coefficient-to-leading-coefficient ratio, an elementary symmetric function of those roots, would remain bounded. This contradiction forces an unstable root.

An important implicit example is BDF2 method, whose characteristic equation is $(3-2z){\zeta }^2-4\zeta +1=0$. A unit root $\zeta =e^{i\theta }$ givesThe leading coefficient vanishes only at $z=3/2$, outside the left half-plane. Near a small negative $z$, both roots are inside the disk; no unit-circle crossing can occur in the connected open left half-plane. On the imaginary axis the boundary equality can occur only at $z=0$, where the roots are $1,1/3$. Hence BDF2 is A-stable and its two roots tend to zero for very large negative $z$, giving strong stiff-mode damping. This is a complete boundary-locus argument, not just a boundary plot.

The Second Dahlquist barrier says that a convergent A-stable ordinary linear multistep method has order at most two. This theorem does not limit implicit Runge-Kutta methods or multiderivative multistep methods. In particular the A-stable third-order two-step multiderivative method in question 1 uses $y^{\prime \prime }$ and satisfies different order conditions. Among collocation Runge-Kutta methods, the standard family results give Gauss collocation methods order $2s$ with diagonal $[s/s]$ Padé approximants to $e^z$ as stability functions, and Radau IIA methods order $2s-1$ with subdiagonal $[s-1/s]$ approximants. Both families are A-stable; the Gauss functions have a nonzero limit at infinity, whereas Radau IIA is L-stable. The corrected two-stage right-Radau method, for example, haswhich matches the exponential through cubic order and tends to zero at infinity. Its poles $2\pm i\sqrt{2}$ lie in the right half-plane. For $z=-r+iy$, $r\ge 0$, direct subtraction givesproving its A-stability directly.

Sectorial A-alpha stability relaxes the requirement to a sector around the negative real axis, useful when the spectrum of a stiff problem remains in that sector; higher-order backward differentiation formulas illustrate the usefulness of this weaker requirement. For a normal linear system $y^{\prime }=Ly$, unitary diagonalization reduces the numerical update to its scalar factors $R(h{\lambda }_j)$. For a nonnormal system, a scalar eigenvalue test alone does not give a uniform bound in the chosen norm. Nonlinear stability is different again: B-stability asks for contractivity when $\mathrm{Re}⟨f(u)-f(v),u-v⟩\le 0$, and algebraic stability of a Runge-Kutta method is a standard sufficient criterion. Thus A-stability addresses a precise and valuable linear test, while step selection still has to meet accuracy requirements and the actual implicit equations must be solved on the relevant branch.

Eigenvalue analysis proves mesh-uniform stability when it is accompanied by control of eigenvectors or an energy norm; the signs of eigenvalues alone are insufficient. For an evolution equation, the method of lines first produces $u_h^{\prime }=L_h{u}_h$. A time integrator then gives an amplification matrix $S_{h,\delta }$, for example $S_{h,\delta }=R(\delta L_h)$ for a Runge-Kutta method. The useful stability requirement on a fixed physical time interval iswhere $C_T$ is independent of the refining spatial and temporal meshes. The semidiscrete analogue is $\Vert e^{t{L}_h}{\Vert }_h\le C_T$ for $0\le t\le T$. These are bounds on propagation of initial errors and residuals, not merely statements that each finite matrix has a bounded solution over its own finite time interval.

If $L_h$ is a normal matrix in the chosen discrete inner product, it has a unitary eigenbasis. ThenFor dissipative problems it is sufficient that every $\delta {\lambda }_j$ lie in the integrator's linear stability domain. A diagonalizable family also suffices if its eigenvector matrices have uniformly bounded condition numbers: the same maximum is multiplied by $\Vert V_h\Vert \Vert V_h^{-1}\Vert$. Without that uniform bound the argument does not establish mesh stability. For finite-time estimates, an amplification modulus $1+O(\delta )$ can also be acceptable; strict contraction at every step is a stronger property.

For a constant-coefficient stencil on the full lattice or a periodic grid, von Neumann stability analysis provides the eigenvectors explicitly. A Fourier mode diagonalizes each translation-invariant spatial difference, giving a scalar symbol ${\lambda }_h(\theta )$ or a small matrix symbol for a system. The discrete Fourier transform and Parseval identity turn a uniform modal bound into an actual discrete L2 norm bound. The two-dimensional stencil in question 3 has $\mathrm{Re}{\lambda }_h=-(1-\cos \eta {)}^2/h\le 0$, so its continuous-time amplification is contractive; a later time integrator must still be checked. With boundaries one must analyze the boundary rows as well: a whole-line Fourier symbol cannot detect an unstable boundary closure.

For the heat equation $u_t=\kappa u_{xx}$ on a unit interval with zero Dirichlet data, the centred second difference has a sine eigenbasis and eigenvaluesExplicit Euler method multiplies mode $j$ by $1+\delta {\lambda }_j$. Requiring all factors to have modulus at most one yields $\kappa \delta /h^2\le 1/[2{\sin }^2((N-1)\pi /(2N))]$ on a given grid, and the sharp grid-independent sufficient limit isThese parabolic time-step restrictions express the growth of the largest-magnitude diffusion eigenvalue like $h^{-2}$. Backward Euler method instead has factor $(1-\delta {\lambda }_j{)}^{-1}$ and the Crank-Nicolson method factor $(1+\delta {\lambda }_j/2)/(1-\delta {\lambda }_j/2)$, both bounded by one for every $\delta \ge 0$. Their unconditional stability follows from A-stability, with much stronger damping of the shortest scales for backward Euler.

For $u_t=-a{u}_x$, $a>0$, explicit time stepping with a backward spatial difference givesThus the upwind finite difference scheme is contractive exactly for $0\le \nu \le 1$, a Courant–Friedrichs–Lewy condition matching the direction and speed of information propagation. A centred spatial difference has purely imaginary eigenvalues, and forward Euler gives $|S{|}^2=1+{\nu }^2{\sin }^2\theta >1$. Under a fixed positive hyperbolic Courant number, high-frequency data grow by a fixed factor in each of $O(1/h)$ steps, proving instability as the mesh is refined. This is not a claim that no exceptionally small coupled step can ever give a finite-time bound: if $\delta =O(h^2)$, the excess per step is only $O(\delta )$, but the usual hyperbolic scaling is unstable. Higher-order explicit integrators can possess nontrivial imaginary-axis stability intervals and are then suitable for skew-adjoint advection or the energy-form first-order wave equation.

Variable coefficients usually prevent Fourier diagonalization, but may preserve a self-adjoint spatial operator. The positive-face diffusion discretization in question 5 is an example. Its discrete energy identity makes $A+B$ symmetric negative definite, so a unitary eigenbasis still exists even without a closed eigenvalue formula. More generally, a dissipative operator $L_h$ in a discrete energy inner product obeyswhich directly proves a uniform semigroup estimate. The dissipative Cayley-transform contraction also follows without normality: for $\alpha \ge 0$,Taking $v=(I-\alpha L_h{)}^{-1}u$ proves contraction of the Crank–Nicolson factor. The contractivity of split Crank-Nicolson diffusion then follows from multiplying two contraction bounds; simultaneous eigenvectors of the two directional matrices are unnecessary.

The negative spectra do not imply uniform semidiscrete stability example shows what can go wrong. LetBoth eigenvalues are $-1$ for every mesh, but applying the propagator to the second unit vector at $t=1$ gives norm at least $1/(eh)$. The evolution is therefore not mesh-uniformly stable. For a discrete matrix, a unit-modulus eigenvalue with a nontrivial Jordan block gives powers growing polynomially in the step number. Even strictly interior eigenvalues can allow large transient growth when eigenvectors are badly conditioned. A pseudospectrum or numerical range of an operator, and especially a direct energy estimate, can reveal behavior missed by eigenvalues.

A resolvent estimate provides another way to see the needed uniformity. If $\Vert S_h^n\Vert \le C$ for every $n$, then for $|z|>1$,This shows directly that a uniform power estimate imposes more than the location of the spectrum. Conversely, resolvent-to-power estimates in varying dimension must themselves be uniform; a bound with constants depending on the number of grid points is not enough.

For a linear multistep method, each spatial eigenvalue gives an amplification polynomial, and all roots must satisfy the root condition for a multistep method, including simplicity of any unit-modulus root. Merely following the root approximating $e^{\delta \lambda }$ misses parasitic modes. The temporal and spatial analyses must therefore be combined with the appropriate zero-stability and uniform modal bounds.

Finally stability connects to convergence. If the exact grid samples have normalized residual ${\tau }^n$ and the error obeys $e^{n+1}=S_h{e}^n+\delta {\tau }^n$, iteration givesA residual tending uniformly to zero then gives convergence on $[0,T]$. For a well-posed linear initial-value problem and a consistent linear finite-difference approximation, this is the forward implication of the Lax equivalence theorem. The practical conclusion is to analyze the full operator including boundaries, select a time method whose stability region contains the relevant scaled spectrum, and verify that eigenvectors or an energy estimate make the resulting bounds uniform as the grid is refined.

One convenient version of Watson's lemma is the following. Suppose $G(t)$ is integrable away from zero and has, as $t\downarrow 0$, an asymptotic expansion $G(t)\sim {\sum }_{j\ge 0}{c}_j{t}^{\beta +j-1}$, where $\beta >0$. On a finite interval, or on an infinite interval with suitable exponential growth control,To prove a finite truncation, split the integral at a fixed small $\delta >0$. The part away from zero is exponentially small. On $(0,\delta )$ subtract the first $N$ local terms and bound the remainder by $C{t}^{\beta +N-1}$; its integral is $O({\lambda }^{-\beta -N})$. Extending each polynomial integral to infinity changes it only exponentially, and the substitution $u=\lambda t$ produces the Gamma integral. A little-o local remainder gives the corresponding little-o asymptotic remainder. This is a proof by localization and remainder bounds, rather than merely formal integration of a series.

For the monotone-phase integral, assume the derivatives needed below exist and are bounded near $a$, and the rest of the integral is controlled. Put $t=\varphi (a)-\varphi (x)$. Since ${\varphi }^{\prime }(a)<0$, this is a regular endpoint coordinate andDifferentiating with $dx/dt=-1/{\varphi }^{\prime }$ gives $G(0)=-g(a)/{\varphi }^{\prime }(a)$ and $G^{\prime }(0)=g^{\prime }(a)/{\varphi }^{\prime }(a{)}^2-g(a){\varphi }^{\prime \prime }(a)/{\varphi }^{\prime }(a{)}^3$. Watson's lemma therefore gives the two-term monotone-endpoint Laplace expansionAn $O(e^{\lambda \varphi (a)}{\lambda }^{-3})$ remainder follows, for example, from a sufficiently smooth transformed amplitude with bounded second derivative near zero. The displayed second coefficient can vanish for special amplitudes without changing the expansion.

The unheaded transition request concerns a cubic endpoint-to-saddle transition. Near $x=0$ its exponent has the expansionAt the cubic endpoint $x$ has size ${\lambda }^{-1/3}$. Requiring the linear term to contribute on the same scale gives $1-\alpha =O({\lambda }^{-2/3})$. Thus the distinguished limit and its uniform leading integral areFor fixed $\nu$, the substitution $x={\lambda }^{-1/3}t$ makes the omitted exponent terms $O({\lambda }^{-2/3})$ on bounded $t$ intervals. Cubic decay controls the remaining tail, while the original interval away from zero is exponentially negligible. The estimate is uniform for $\nu$ in compact sets. In the printed range $\alpha \le 1$, $\nu \ge 0$; the same transition integral also makes sense for negative fixed $\nu$.

At $\nu =0$, the Gamma integral gives $h(0)=6^{1/3}\Gamma (1/3)/3$, recovering the cubic endpoint regime. For $\nu \to +\infty$, scale $t=u/\nu$ and expand $e^{-u^3/(6{\nu }^3)}$ to obtain $h(\nu )={\nu }^{-1}-{\nu }^{-4}+O({\nu }^{-7})$. Its leading term gives ${\lambda }^{-1/3}/\nu =1/[\lambda (1-\alpha )]$, recovering the ordinary endpoint regime in their overlap. Both matching limits are therefore checked explicitly.

Use the real diffusive interpretation $\gamma >0$, $U,{\mu }_0\in \mathbb{R}$, with $\nu >0$. Without a positive real diffusion coefficient the printed ordering inequality and real oscillator coordinate would need different hypotheses. Remove the drift by $\eta =e^{Ux/(2\gamma )}v$. Direct differentiation reduces the linear complex Ginzburg-Landau equation toSet $v=b(\xi )e^{-i({\omega }_0+ϵ{\omega }_1)t}$ and $\xi =\sqrt{ϵ}fx$. The resulting equation isCanceling its order-one coefficient and normalizing the quadratic term givesThe supplied decay quantization has $\lambda =2j+1$, $j=0,1,2,\dots$. Equivalently the Hermite functions $b_j=e^{-{\xi }^2/2}{H}_j(\xi )$ solve this harmonic-oscillator problem, where $H_j$ is the physicists' Hermite polynomial. Thus ${\omega }_1=-i(2j+1)\sqrt{\nu \gamma }$ and the exact temporal growth rates areThe global modes of a quadratically confined Ginzburg-Landau equation areFor every fixed positive $ϵ$, the negative quadratic spatial exponent dominates the linear drift exponent and the polynomial at both infinities; all these normal modes are spatially bounded. The largest growth rate is $s_0$, so the least-node normal mode is temporally non-growing, and all normal modes are non-growing, whenEquality is a neutral ground-state normal mode; strict inequality gives decay. This proves the printed sufficient condition for a bounded stable global normal mode and identifies its stronger meaning: stability of the leading global normal mode. Existence of merely some stable higher-$j$ normal mode is not equivalent to the inequality, since sufficiently large $j$ gives $s_j<0$ even above the leading-mode threshold. In the quadratic-coefficient model the reduction and eigenvalues are exact, not just a first-order approximation to a more general profile.

Use a Fourier transform parallel to the boundary. Splitting its traces into right- and left-supported parts gives Half-range Fourier transforms $F_{+}$ and $F_{-}$, analytic respectively above and below a common inversion contour. A small positive imaginary part of the wavenumber implements the limiting absorption principle and separates the incident poles and outgoing branch points. Solve the transformed differential equation normal to the boundary using its outgoing or decaying solution. The boundary conditions then yield a Wiener-Hopf equation of the formFactor its Wiener-Hopf kernel as $K=K_{+}{K}_{-}$, with each factor analytic and nonzero in its designated half-plane. Divide by $K_{-}$ and split $H/K_{-}=H_{+}+H_{-}$, allocating poles according to the incident-wave prescription. ThenThe two sides continue to a common entire function. Their growth at infinity, inherited from the field's edge behavior, determines that function, usually a constant or a polynomial; in a decaying case it is zero. Solve for the unknown transforms and apply Fourier inversion, preserving the outgoing contour prescription.

The principal analytic results used are as follows. The Paley–Wiener theorem and its half-line Laplace-transform version relate support on a half-line to holomorphy in a half-plane, with growth or square-integrability bounds. The Cauchy integral formula and Sokhotski–Plemelj theorem provide additive splitting: a Cauchy integral of a sufficiently decaying jump has upper and lower boundary values with that prescribed difference. Applying the same construction to a continuous logarithm of a nonvanishing scalar Wiener-Hopf kernel gives Wiener-Hopf factorization when its winding number is zero; a nonzero index requires an explicit index factor. The identity theorem gives analytic continuation across the common strip. Finally Liouville's theorem makes a bounded entire function constant; its polynomial-growth extension bounds the possible degree. The physical edge and radiation conditions select the remaining constants. Zeros, branch cuts, nonzero index and insufficient decay must be accounted for rather than discarded during splitting.

For $\alpha >0$, minimize $J_{\alpha }(x)=\Vert Ax-y{\Vert }^2+\alpha \Vert x{\Vert }^2$. Its first variation in every direction givesThe self-adjoint operator on the left has lower bound $\alpha I$, so it is invertible and the minimizer is unique. The regularized solution has the singular expansionSince ${\sigma }^2+\alpha \ge 2\sigma \sqrt{\alpha }$, its gain is at most $1/(2\sqrt{\alpha })$. ConsequentlyThis proves existence, uniqueness and continuous data dependence at fixed positive parameter. For the particular perturbation from part (i), the error is $\delta {\sigma }_j/({\sigma }_j^2+\alpha )$ instead of $\delta /{\sigma }_j$. To recover a minimum-norm exact solution as noise tends to zero, choose $\alpha \to 0$ while $\delta /\sqrt{\alpha }\to 0$; dominated convergence of the exact-data singular coefficients removes the bias.

To connect this with the provided inverse expansion, put $B=A{A}^{\ast }$ and $C=A^{\ast }A$. The exact intertwining identity isIt follows by multiplying $(C+\alpha I)A^{\ast }=A^{\ast }(B+\alpha I)$ on either side by the appropriate resolvent. A useful exact inverse identity isExpanding $(I+{\alpha }^{-1}C{)}^{-1}$ to two terms gives the printed approximation, namelyThis is the truncated Tikhonov resolvent. Multiplication gives $(\alpha I+B)R_{\alpha ,2}=I+{\alpha }^{-3}{B}^3$, so the expression is not an exact inverse. For example $A=1$ and $\alpha =1$ give the approximate value $1$ instead of the exact $1/2$. Its exact remainder isIn the regime $\Vert C\Vert /\alpha <1$ it is a controlled Neumann expansion. Using it to reconstruct givesFor each fixed $\alpha >0$ this is a bounded polynomial in bounded operators, withThus the supplied approximation also has continuous data dependence in its domain of use. It must not be extrapolated to $\alpha \downarrow 0$ with $A$ fixed: its gain on any nonzero singular mode diverges, so it would not define a convergent regularization family there. The exact Tikhonov filter supplies the well-posed regularized problem; the truncated formula is only a controlled computational approximation.

The Vandermonde determinant isIndeed, the determinant is an alternating polynomial, hence divisible by every $x_j-x_i$. Both sides have degree $r(r-1)/2$, and comparison of the coefficient of $x_0^0{x}_1^1\cdots x_{r-1}^{r-1}$ gives constant one. This is an identity over the integers, so reducing its coefficients modulo a prime number preserves it. For $r>p$, two residues coincide and the determinant vanishes. For $r=p$, choose all elements of ${\mathbb{F}}_p$ once to obtain a nonzero determinant.

For the ordering $0,1,\dots ,p-1$, this determinant is $D={\prod }_{j=1}^{p-1}j!$. If $p$ is odd, put $h=(p-1)/2$. Wilson theorem gives $j!(p-1-j)!=(-1{)}^{j+1}$. Pairing the factors with indices $j=0,\dots ,h-1$ and $p-1-j$, and leaving $h!$ unpaired, givesEvery distinct ordering changes the determinant by its permutation sign; repeated residues give zero. Thus the possible values are $0,\pm h!$ for odd $p$, and $0,1$ for $p=2$.

Use Shamir's secret sharing: choose independent uniform $a_1,\dots ,a_{r-1}\in {\mathbb{F}}_p$, and engrave the shares $(x_i,P(x_i))$, where $P(X)=N+{\sum }_{j=1}^{r-1}{a}_j{X}^j$ and the $x_i$ are distinct nonzero field elements. We may take $p>n$; otherwise use a sufficiently large finite field extension and an agreed integer encoding. Any $r$ shares determine $P$ by polynomial interpolation and hence $N=P(0)$. Given any $r-1$ shares and any proposed $N$, their equations in the coefficients have a unique solution: the coefficient matrix has determinant ${\prod }_i{x}_i$ times a nonzero Vandermonde determinant. Since the coefficients were uniform, the observed shares have the same distribution for every $N$. Thus they give no information about the secret.

For the fake, reduce its integers to field elements and first check that its abscissa is distinct and nonzero. Random collisions or a zero abscissa are immediately detectable and have probability $O(n/p)$; the following conclusions concern the overwhelmingly likely distinct-abscissa case.

With $r$ shares, there is always a unique degree-at-most-$r-1$ interpolating polynomial. If all shares are genuine it gives the correct secret. If one is fake, it normally gives the wrong secret, with probability $1-1/p$ conditional on a uniform fake ordinate, but the inconsistency is undetectable from those $r$ shares. The data do not certify which case has occurred.

With $r+1$ shares, genuine shares are consistent with one such polynomial. A fake is inconsistent with probability $1-1/p$, so it is usually detected as an inconsistency. However, omitting each of the $r+1$ shares in turn gives $r+1$ possible polynomials. If the data are inconsistent, their constant terms are distinct: two candidate polynomials with the same constant term agree at zero and at their $r-1$ common abscissae, hence agree identically, which would make all the data consistent. The secret is narrowed to $r+1$ possibilities, but the fake is not identified.

With $r+2$ shares, at least $r+1$ are genuine. There is a unique degree-at-most-$r-1$ polynomial agreeing with at least $r+1$ shares: two candidates would agree on at least $r$ abscissae and hence coincide. Find it by deleting each share and checking consistency. The secret is recovered uniquely, and the fake is identified unless it happens to lie on the genuine polynomial, in which case it is harmless.