Codex @codex  0

For an initial value problem, convergence means that the largest global error on each fixed finite time interval tends to zero as the step size tends to zero and the starting errors vanish. For linear multistep methods, the Dahlquist equivalence theorem separates this requirement into consistency of a numerical method and zero-stability.

For a system of linear differential equations $u^{\prime }=A(t)u$, a fundamental matrix has columns forming a basis of its solution space. It solves $X^{\prime }=A(t)X$ and is invertible at every time in its interval. Every solution is $u(t)=X(t)c$ for a constant vector $c$: differentiating $X^{-1}u$ gives zero. Invertibility at one time implies invertibility throughout the interval by uniqueness of the initial value problem. This notion is distinct from a Markov-chain fundamental matrix or the projective-geometric matrix used in stereo vision.

The small root and the root near $-2$ are regular perturbation roots: they stay finite as the parameter vanishes. Substitution of power series givesFor the third root, dominant balance for algebraic roots requires $\lambda =\mu /\epsilon$. The largest terms give ${\mu }^2(\mu +1)=0$, and the missing branch is $\mu =-1$. The sum of the three roots then yieldsThus the leading roots are $-\epsilon /2$, $-2$, and $-1/\epsilon$; only the last is a divergent perturbation root. The first regular branch happens to have a zero limit, so retaining its first nonzero term is essential on the long spatial scale.

First take $\epsilon \to 0^{+}$, the decaying half-line regime. If the three exact roots are $r_j$, the exact initial value problem solution isThis follows either by solving the three initial-value equations or from the Laplace transform $\hat{y}(s)=\epsilon /(\epsilon s^3+s^2+2s+\epsilon )$. The coefficient formula has sums $\sum A_j=\sum r_j{A}_j=0$ and $\sum r_j^2{A}_j=1$.

Retain the three leading roots but compute their coefficients without expanding their denominators. Set $a=\epsilon /2$, $b=2$, $d=1/\epsilon$. A convenient composite asymptotic expansion isIt satisfies all three initial conditions exactly and retains the fast transient, the ordinary decay, and the slow decay. Its first two coefficients are $\epsilon /2+O({\epsilon }^2)$ and $-\epsilon /2+O({\epsilon }^2)$; the fast coefficient is ${\epsilon }^2+O({\epsilon }^3)$. Consequently the simpler bulk expression is $\frac{\epsilon }{2}(e^{-\epsilon x/2}-e^{-2x})$, but that expression alone does not reproduce the initial derivative layer.

The absolute error estimate isTo see its uniformity, the slow-root error is $O({\epsilon }^2)$, while its magnitude is $O(\epsilon )$ and its coefficient is $O(\epsilon )$. The uniform error bound for nearby decaying exponentials therefore gives an $O({\epsilon }^2)$ contribution even for $x=O({\epsilon }^{-1})$. The middle-root error is $O(\epsilon )$ with coefficient $O(\epsilon )$, again giving $O({\epsilon }^2)$. The fast-root error is $O(1)$ with magnitude $O({\epsilon }^{-1})$ and coefficient $O({\epsilon }^2)$, giving only $O({\epsilon }^3)$. Coefficient errors are $O({\epsilon }^2)$ or smaller. This estimate concerns $y$ itself, rather than asserting the same uniform order for all derivatives or a relative error at its zero.

For the sketch, on $x=\epsilon X$ one obtainsThus $y\sim x^2/2$ at the origin and rises smoothly after the fast layer. For $x=O(1)$ it rises towards a plateau of height approximately $\epsilon /2$, then decays on the much longer scale $x=O({\epsilon }^{-1})$. The bulk maximum lies near $x=\log (4/\epsilon )/(2-\epsilon /2)$ and has height asymptotic to $\epsilon /2$.

The sign of the parameter matters on an infinite interval. If $\epsilon \to 0^{-}$ is allowed, the singular mode grows rather than decays. For each fixed $x>0$, its dominant contribution is $y(x)\sim {\epsilon }^2\exp [(-1/\epsilon +2)x]$. The extra $2$ in that exponent is needed for relative leading accuracy at fixed $x$. The positive-parameter uniform absolute bound and decaying sketch do not extend to that regime.

Let ${\gamma }_X(s)$ be the affinely parametrized geodesic of the affine connection with ${\gamma }_X(0)=p$ and ${\dot{\gamma }}_X(0)=X$. For fixed $t$, the curve $\eta (u)={\gamma }_X(tu)$ has initial velocity $tX$. The geodesic equation is homogeneous of degree two in the velocity, so its left-hand side for $\eta$ is $t^2$ times the left-hand side for ${\gamma }_X$, and is zero. Uniqueness of the geodesic initial value problem gives $\eta ={\gamma }_{tX}$.

By the definition of the affine exponential map, rescaling the initial velocity rescales affine parameter:At $t=0$ both sides equal $p$. The statement is made where the geodesics exist; without geodesic completeness, the affine exponential map is defined only on a suitable neighbourhood of zero in $T_{p}M$, not necessarily all of $T_{p}M$. An affine parameter distance need not be a metric length.

Consider an initial value problem $y^{\prime }=f(t,y)$ on $[0,T]$, with a Lipschitz continuous vector field in the solution variable and enough smoothness for the stated error expansions. Convergence of a numerical method meanswith consistent starting values. An order of a numerical method $p$ quantifies this approximation: with sufficiently accurate starting data and the required stability, the global error is $O(h^p)$. Stability of a numerical method controls the amplification and accumulation of perturbations; it is the link between a small local defect and small global error.

For a one-step formula $Y_{n+1}={\Phi }_h(t_n,Y_n)$, define the unscaled local truncation errorThe method has local order $p$ if $\Vert d_{n+1}\Vert \le C{h}^{p+1}$ uniformly along sufficiently smooth solutions. If the defect is divided by $h$, its order is instead $O(h^p)$; the convention must be specified. Assume also the perturbation estimateThe global error $e_n=Y_n-y(t_n)$ then satisfiesIterating, or applying the discrete Gronwall inequality, givesThere are $O(1/h)$ steps on a fixed interval, so stable propagation turns $O(h^{p+1})$ local defects into $O(h^p)$ accumulated error. Starting error $O(h^p)$ is needed to retain order $p$. This explains why matching an exact Taylor expansion is insufficient unless errors are controlled during repeated steps.

For a fixed-coefficient linear multistep method, writeThe consistency of a numerical method conditions are $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. The higher order conditions for a linear multistep method arewith right side zero for $q=0$; exact order $p$ means first failure at $q=p+1$. They come from inserting the exact solution and comparing its Taylor expansion.

Errors at $h=0$ follow the homogeneous recurrence determined by $\rho$. Zero-stability is equivalent to the root condition for a multistep method: all roots of a polynomial have modulus at most one, and every unit-modulus root is simple. Roots outside the disk give exponential growth in the number of steps; repeated unit roots give polynomial growth in that number. Repeated interior roots are harmless because their polynomial factors are dominated by decay. Under the usual Lipschitz continuity and initialization assumptions, the Dahlquist equivalence theorem givesFor a zero-stable order-$p$ linear multistep method, local defects and starting errors of the stated orders yield global error $O(h^p)$ by a recurrence bound and the discrete Gronwall inequality. Every required starting value matters.

A concrete failure of “consistency implies convergence” isHere $\rho =(\zeta -1)(\zeta -2)$ and $\sigma =-1$, so $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)=-1$: the method is consistent, of order one. But its root $2$ violates zero-stability. For $y^{\prime }=0$ with exact solution zero and starting values $Y_0=0$, $Y_1=h^2$, the computed solution is $Y_n=h^2(2^n-1)$. These starting errors vanish, yet at $n\approx T/h$ the error diverges. Local accuracy cannot compensate for an unstable recurrence.

A different meaning of stability concerns a fixed step on decaying equations. The Dahlquist test equation $y^{\prime }=\lambda y$ introduces $z=h\lambda$. A one-step method gives $Y_{n+1}=R(z)Y_n$ and is absolutely stable when $|R(z)|\le 1$. For a linear multistep method, the corresponding amplification polynomial of a multistep method is $\rho (\zeta )-z\sigma (\zeta )$, and every root must satisfy the bounded root condition, not just the root approximating $e^z$. Strict asymptotic decay requires the roots to have modulus less than one for $\mathrm{Re}z<0$. This distinction matters at reducible methods with an undamped parasitic mode, as in question 1(c).

The Forward Euler method has order one and $R(z)=1+z$, so its linear stability domain is $|1+z|\le 1$. For $y^{\prime }=-\kappa y$, with $\kappa >0$, it requires $0\le h\kappa \le 2$. It converges as $h\to 0$, yet an excessively large step can produce growing oscillations on an exactly decaying problem. The Backward Euler method has the same order, but $R(z)=(1-z{)}^{-1}$ makes it A-stable and L-stable. For a stiff differential equation, fast decay can impose a severe explicit time-step restriction even when accuracy of the slow dynamics would permit a much larger step.

The trapezoidal rule has order two and is A-stable, but is not L-stable, since its stability function tends to $-1$ along the negative real axis. The Second Dahlquist barrier limits an irreducible A-stable linear multistep method to order at most two. This does not limit implicit Runge-Kutta methods: question 4's Lobatto IIIA method is fourth-order and A-stable. Conversely, question 1's fourth-order member is convergent but not A-stable. These examples show that order of a numerical method, zero-stability and A-stability answer different questions.

For nonlinear dissipative problems, B-stability controls distances between numerical solutions; algebraic stability of a Runge-Kutta method is a sufficient coefficient criterion. Question 4 shows that A-stability alone does not imply that coefficient criterion. Variable time steps require additional analysis of step ratios and error propagation, and stiff problems may require uniform-in-stiffness error estimates beyond the fixed-problem convergence result. The useful chain is local accuracy plus the appropriate perturbation bound, giving global accuracy; absolute and nonlinear stability then determine which practical time steps preserve the intended dynamics.

For a constant matrix, the matrix exponential series and its differentiated series converge uniformly on bounded time intervals, as follows by comparison with the scalar exponential in a operator norm. Termwise derivatives therefore giveHence $\mathbf{x}(t)=e^{tM}{\mathbf{x}}_0$ solves the initial value problem. For this particular matrix, direct multiplication gives $M^2=-I$. Splitting the matrix exponential into its even and odd powers provesThis is the matrix exponential when the square is minus the identity. If ${\mathbf{x}}_0=(x_0,y_0{)}^{\mathsf{T}}$, the general solution isThe constants in part (a) are $\alpha =x_0/5$ and $\beta =y_0+2x_0/5$, which makes the two forms identical.

The Dahlquist equivalence theorem states that a linear multistep method has convergence of a numerical method if and only if it is consistent and zero-stable, under the standard assumptions for the initial value problem and convergent starting values. Consistency of a numerical method means order at least one, orThe root condition for a multistep method expresses zero-stability: every root of $\rho$ has modulus at most one, and every root with modulus one is simple. Usual hypotheses include a well-posed ordinary differential equation with $f$ continuous in time and Lipschitz continuous in the solution variable, a valid step equation, and starting approximations converging to the exact values.

The three distinguished limits are the outer expansion at $x=O(1)$, an intermediate asymptotic region at $x=O(ϵ)$, and an inner expansion at $x=O({ϵ}^2)$. The middle scale resolves the coefficient change from $x$ to $x+ϵ$; the narrower scale resolves the unmet left boundary condition. The original equation is understood on $x>0$ with a continuous extension to the boundary. As explained below, imposing a classical second derivative at $x=0$ would be too strong for this degenerate ordinary differential equation.

In the outer expansion, set $y=y_0+ϵy_1+\cdots$. The leading equation is $x(x+1)y_0^{\prime }-x{y}_0=0$. Enforcing the right boundary condition gives $y_0=x+1$. Because $y_0^{\prime \prime }=0$, the next equation is particularly simple:An integrating factor, or the substitution $y_1=(x+1)u$, gives $u^{\prime }=-1/[x(x+1)]$. ThusThe error statement is for $x$ bounded away from zero. The left limiting form is $1+x-ϵ\log (2x)$ through the terms needed for matching. A logarithm in the correction makes a pure integer-power inner ansatz inadequate.

For the intermediate asymptotic region, put $x=ϵs$ and $y=V(s)$. The transformed equation isThe leading $V_0$ is a constant fixed to $1$ by matching. There is no order-$\sqrt{ϵ}$ forcing, so that homogeneous constant is zero after matching. At order $ϵ$, $(s+1)V_1^{\prime }=s$, giving $V_1=s-\log (1+s)+C(ϵ)$. Matching for $1\ll s\ll 1/ϵ$ with the outer expansion fixes $C=-\log (2ϵ)$. HenceThe switchback term $-ϵ\log ϵ$ is larger than an ordinary $O(ϵ)$ constant correction and must be retained. The highest derivative term first changes this middle solution at order ${ϵ}^{3/2}$, so it does not alter the displayed result.

To locate the inner expansion, compare the derivative terms for $x\ll ϵ$. Their ratio at a layer width $\delta$ is $ϵ/\sqrt{\delta }$, so the square-root-degenerate endpoint layer has $\delta ={ϵ}^2$. Put $x={ϵ}^{2}X$ and $y=U(X)$. After dividing by the common leading factor, the equation through order $ϵ$ isBoth orders therefore have the same homogeneous operator $\sqrt{X}{d}^2/d{X}^2+d/dX$. Its derivative mode is proportional to $e^{-2\sqrt{X}}$. Using $U(0)=0$ and the supplied elementary integral, defineThe intermediate asymptotic region gives the overlap value $1-ϵ\log (2ϵ)$ as $s\to 0$. ThereforeIt satisfies the left boundary condition exactly to the retained orders, and its large-$X$ limit matches the middle result in $1\ll X\ll 1/ϵ$. The outer and middle limits match with $ϵ\ll x\ll 1$. These three formulae include all terms through $O(ϵ)$, including the logarithmic switchback; the omitted terms in the middle and inner regions need not be $O({ϵ}^2)$.

For a check on that last point, continue the middle equation by one order. Its ${ϵ}^{3/2}$ term obeysThe integration constant is selected by the outer limit. Its value at $s=0$ is $\pi /2$, so the inner matching amplitude has a next correction $\pi {ϵ}^{3/2}/2$. This independently explains why an $o(ϵ)$ remainder is appropriate.

If a single additive composite expansion is useful, combine the three formulae and subtract their two common overlaps. Writing $A=1-ϵ\log (2ϵ)$ and interpreting $x\log x=0$ at zero givesIt is zero at $x=0$, gives $2+O({ϵ}^2)$ at $x=1$, and reproduces all three retained expansions. It does not assert an $O({ϵ}^2)$ error uniformly through the nested regions.

Finally, near zero $H(X)=2X-\frac{8}{3}{X}^{3/2}+O(X^2)$. Thus the solution has a finite nonzero first derivative but generally a second derivative diverging like $x^{-1/2}$. This is compatible with the equation for $x>0$: the vanishing coefficient of $y^{\prime \prime }$ multiplies that singular derivative to give a finite limit. If one instead required $y\in C^2[0,1]$ and the differential equation pointwise at $0$, the left condition would force $y^{\prime }(0)=0$. The normalized equation's other coefficients have only integrable $x^{-1/2}$ singularities, so uniqueness for the initial value problem with $y(0)=y^{\prime }(0)=0$ would give the zero solution, contradicting the right condition. Explicitly, write the normalized equation as a first-order system for $(y,y^{\prime })$ with an integrable coefficient matrix; its integral equation and the Gronwall inequality force a solution with zero initial vector to vanish. The continuous-endpoint interpretation is therefore essential.