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

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 thatHere 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 $\rho (1)=0$ and ${\rho }^{\prime }(1)=\sigma (1)$. Zero-stability is equivalent to the root condition for a multistep method: every root of a polynomial $\xi$ of $\rho$ satisfies $|\xi |\le 1$, and roots with $|\xi |=1$ are simple. Order $p\ge 1$, together with zero-stability and starting errors $O(h^p)$, gives global error $O(h^p)$ under the usual smoothness assumptions.

To prove the necessity of the root condition for multistep convergence, use the test problem $y^{\prime }=0$, $y(0)=0$. Its discrete error satisfies $\rho (E)e_n=0$, where $E{e}_n=e_{n+1}$. Set $h=T/N$. If $|\xi |>1$, the exact recurrence solutionhas starting errors tending to zero at every fixed starting index, but $|e_N|=1$. Thus convergence fails.

If a unit-modulus root has multiplicity $m\ge 2$, the recurrence instead admitsAgain every fixed starting error tends to zero while $|e_N|=1$. These polynomial-times-exponential solutions follow from the repeated factor $(E-\xi {)}^m$. 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.

Factor the first characteristic polynomial asThe other two roots of a polynomial have product one. If $-1<\alpha <1$, they are $\alpha \pm i\sqrt{1-{\alpha }^2}$: distinct unit-modulus roots, both different from one. At $\alpha =-1$, the root $-1$ is double; at $\alpha =1$, the root one is triple. Outside $[-1,1]$, one of the reciprocal real roots has modulus greater than one. Thus the root condition for a multistep method holds exactly for $-1<\alpha <1$.

For the order conditions for a linear multistep method, expand the exponential defect:The constant, linear and quadratic coefficients vanish for every $\alpha$. The cubic coefficient vanishes only at $\alpha =11/5$, where the quartic coefficient is $-3/5\ne 0$. Consequently the formal order isFor all parameters satisfying the root condition for a multistep method, ${\rho }^{\prime }(1)=\sigma (1)=2(1-\alpha )\ne 0$, so this formal calculation gives ordinary second-order consistency of a numerical method. The Dahlquist equivalence theorem therefore yieldsThe exceptional third-order formula is not zero-stable. At $\alpha =1$, $\sigma \equiv 0$ and $\rho =(w-1{)}^3$: the recurrence contains no vector-field evaluations. It satisfies the displayed Taylor identities only formally and is a degenerate, nonconvergent formula, not a usable second-order ODE solver. This distinction avoids interpreting formal cancellation as convergence.

Only $-1<\alpha <1$ needs consideration. On the Dahlquist test equation, the amplification polynomial of a multistep method is $P_z(w)=\rho (w)-z\sigma (w)$, where $z=h\lambda$. A-stability requires all amplification roots to have modulus at most one for every $\mathrm{Re}z\le 0$, with simple unit-modulus roots.

Choose the negative real value $z=-2/(1-\alpha )$. The leading coefficient of $P_z$ iswhereas $P_z(1)=-z\sigma (1)=4>0$. Hence $P_z(w)\to -\infty$ as real $w\to +\infty$. The intermediate value theorem supplies a real amplification root $w>1$, so this negative test value is unstable. ThereforeThis is a direct absolute stability obstruction, independent of the Second Dahlquist barrier, which by itself would not rule out the convergent second-order members.

The nodes $0,1/2,1$ have Lagrange interpolation polynomialsThe coefficients obtained from $a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds$ and $b_j={\int }_0^1{\ell }_j(s)ds$ reproduce the given Butcher tableau. In particular, the second row's last entry is $-1/24$, which is missing from the TeX transcription. Thus this is the three-stage Lobatto IIIA method.

One applicable collocation theorem is that $s$-stage Lobatto IIIA collocation has order $2s-2$. Alternatively, the fourth-order conditions for a Runge-Kutta method give a direct verification. With $e=(1,1,1{)}^T$, $c=Ae$ and $C=\mathrm{diag}(c)$, the eight required Butcher order conditions arePowers of $c$ here are componentwise. Substitution satisfies all eight. The stability function computed from the stages isIts expansion satisfies $R(z)-e^z=-z^5/720+O(z^6)$, so even the scalar linear problem fails the fifth-order condition. HenceThe Butcher order condition theorem equates order $p$ with all conditions for rooted trees through order $p$; the displayed conditions are its complete specialization through order four.

Use the stability function $R(z)=(z^2+6z+12)/(z^2-6z+12)$. Its denominator vanishes only at $3\pm i\sqrt{3}$, both in the open right half-plane. For $z=x+iy$,This is nonnegative whenever $x\le 0$, so $|R(z)|\le 1$ throughout the closed left half-plane. The internal stage equations are also solvable there, since $\det (I-zA)=1-z/2+z^2/12$. Therefore the method is A-stable.

The inequality is strict in the open left half-plane, while $|R(iy)|=1$. Moreover $R(z)\to 1$ as $|z|\to \infty$, so the method is not L-stable. These conclusions follow directly from the definitions of A-stability and L-stability, rather than from the order of the method.

For algebraic stability of a Runge-Kutta method, the weights must be nonnegative andmust be a positive semidefinite matrix. Here all weights are positive, but direct calculation givesIn particular $e_1^{T}M{e}_1=-1/36<0$; equivalently, its eigenvalues are $0,\pm \sqrt{3}/36$. Therefore the method is not algebraically stable. This does not contradict its A-stability: algebraic stability is a stronger condition designed to guarantee nonlinear B-stability through the Runge-Kutta contractivity identity.

The original PDF puts both neighbouring values at time $n+1$. Thus the actual scheme is Backward Euler method in time with a centered three-point approximation to the second spatial derivative. The old-time neighbours in the TeX transcription describe a different method.

Let $h=\Delta x$ and $k=\Delta t$. Insert a sufficiently smooth exact solution and divide by $k$. Expanding at $(x_m,t_{n+1})$ givesThe PDE cancels the $u_t-u_{xx}$ term. Therefore the normalized local truncation error is $O(k+h^2)$ andFor compatible smooth data with exact nodal initialization, the stability of a numerical method established below gives global error $O(k+h^2)$ on fixed time intervals, by summing the one-step defects with the contraction bound. Thus no relation between $k$ and $h$ is needed for stability; choosing $k=O(h^2)$ makes both consistency contributions $O(h^2)$.

Let $T$ be the $M\times M$ tridiagonal matrix with diagonal two and adjacent entries minus one, incorporating zero Dirichlet boundary conditions. The actual PDF scheme isAn orthonormal basis of eigenvectors of $T$ has components $\sqrt{2/(M+1)}\sin (jm\pi /(M+1))$, and substitution givesThey form an orthonormal basis by the spectral theorem for real symmetric matrices. Hence $I+\mu T$ is invertible for every $\mu >0$, and the amplification eigenvalues are $g_j=(1+\mu {\lambda }_j{)}^{-1}\in (0,1)$. For the discrete L2 norm $\Vert U{\Vert }_h^2=h{\sum }_m|U_m{|}^2$,This bound has constant one independent of time step and mesh, which is the required stability of a numerical method. Therefore every positive Courant number is stable. If a forcing or local-defect sequence is added, the discrete variation-of-constants formula gives a bound by the initial norm plus the sum of those perturbation norms; this also justifies the global error conclusion in part (a).

Use $N$ distinct periodic grid points $x_m=m/N$, $0\le m<N$, and $h=1/N$. Values at the two ends of the interval represent the same point, so the endpoint must not be stored twice. With indices interpreted modulo $N$, setThis replaces the Dirichlet matrix by the periodic discrete Laplacian. The normalized Fourier modes $N^{-1/2}{e}^{2\pi ijm/N}$ form an orthonormal basis and have amplification factorsEach satisfies $0<g_j\le 1$ for every $\mu >0$. The Parseval identity therefore gives $\Vert U^{n+1}{\Vert }_h\le \Vert U^n{\Vert }_h$ in the discrete L2 norm, uniformly in the grid and time step. Thus all positive Courant numbers remain stable. The constant mode has $g_0=1$, expressing preservation of the spatial mean; all nonconstant modes decay. The formula with modulo indices also handles small grids, where two neighbours can coincide.

For a sufficiently regular solution with zero Dirichlet boundary conditions, differentiate its squared L2 norm. Since $u_t=-i(\Delta u-Vu)$,Integration by parts and the zero boundary trace give $\int \overline{u}\Delta u=-\int |\nabla u{|}^2$. The potential contribution is real because $V$ is real. The quantity in parentheses is consequently real, and its product with $-i$ has zero real part. HenceThis is the energy method for a skew-adjoint evolution. For general L2 space data, the same result follows from the strongly continuous unitary group generated by a self-adjoint realization of $\Delta -V$. A sufficient assumption is a bounded real potential with the Dirichlet Laplacian; arbitrary real potentials require appropriate domain and self-adjointness assumptions. The differentiated proof presumes the regularity needed for the displayed integrals, rather than asserting it for every pointwise real function $V$.

Arrange the interior values into a vector $U$. Represent the five-point Dirichlet Laplacian as a Kronecker sum $L_h=h^{-2}(D\otimes I+I\otimes D)$, where $D$ is the tridiagonal matrix with diagonal $-2$ and adjacent entries one. This Kronecker sum is real symmetric; the real sampled potential has a real diagonal matrix $D_V$. ThusThe generator $-i{H}_h$ is a skew-Hermitian matrix. ConsequentlyThe matrix exponential is a unitary matrix, either by differentiating its product with its adjoint or by unitary diagonalization of a normal matrix. For the two-dimensional discrete L2 norm $\Vert U{\Vert }_h^2=h^2{\sum }_{m,n}|U_{m,n}{|}^2$, this givesApplying the same equation to a difference proves stability of a numerical method with a mesh-independent constant one. This is norm conservation of a semidiscrete Schrödinger equation. It concerns continuous time after spatial discretization; an arbitrary subsequent time integrator need not preserve this stability.

The printed coordinates do not discretize the stated square. For $[-1,1{]}^2$ with $M$ interior points in each direction, use $h=2/(M+1)$, $x_m=-1+mh$, $y_n=-1+nh$, and sample $V$ there. The printed $h=1/(M+1)$ and unshifted coordinates instead describe a grid on $[0,1{]}^2$. The matrix proof is valid for either geometry with its corresponding boundary values, so this transcription-independent statement flaw does not alter the stability conclusion.

Work in the real zero-boundary Sobolev space $X=H_0^1(0,1)$, whose members have square-integrable weak derivatives and zero endpoint trace. Define the bounded bilinear form and linear functionalThe weak formulation is to find $u\in X$ with $a(u,v)=\ell (v)$ for every $v\in X$. Integration by parts recovers the differential equation for a smooth solution; conversely, this identity is the definition of its weak solution.

Use the full $H^1$ norm. The Cauchy-Schwarz inequality gives $|a(u,v)|\le \Vert u{\Vert }_{H^1}\Vert v{\Vert }_{H^1}$ and $|\ell (v)|\le \Vert v{\Vert }_{H^1}/\sqrt{3}$. Also $a(v,v)=\Vert v{\Vert }_{H^1}^2$, so $a$ is a coercive bilinear form with constant one. The Lax-Milgram theorem states that a bounded, coercive bilinear form on a real Hilbert space determines a unique weak solution for every bounded linear functional. It therefore applies here.

The equivalent minimization problem isIf $u$ is the weak solution, symmetry and $a(u,v-u)=\ell (v-u)$ giveThis is positive unless $v=u$, proving a unique global minimizer. Conversely, differentiating $\mathcal{E}(u+tv)$ at $t=0$ gives the weak identity. Thus the unique energy minimizer is exactly the unique weak solution. As a check, the classical representative isIt satisfies both endpoint conditions and $-u^{\prime \prime }+u=x$.

The trial space must be a subspace of $H_0^1(0,1)$, not merely of $L^2(0,1)$. The printed assumptions on the basis do not ensure square-integrable weak derivatives, and boundary traces are not defined for arbitrary L2 space functions. For example, $\phi (x)=x^{1/4}(1-x)$ is continuous, vanishes at both endpoints and belongs to $L^2$, but $|{\phi }^{\prime }(x){|}^2\sim x^{-3/2}/16$ is not integrable at zero. Assume the basis belongs to the zero-boundary Sobolev space, as it does for the hat functions in part (c).

Write $u_H={\sum }_{j=1}^M{c}_j{\phi }_j$. Differentiating the restricted energy with respect to $c_i$, or imposing the Galerkin method with test function ${\phi }_i$, givesThe derivative term is the usual diffusion stiffness matrix and the second is the mass matrix. The full matrix is real symmetric, and for every nonzero coefficient vector $c$,because a basis is linearly independent. Thus it is a positive-definite matrix, so the restricted Ritz method has a unique solution and minimum.

The piecewise-linear hat functions have disjoint supports unless their nodes coincide or are neighbours. Direct integration on the two adjacent intervals givesSince $x=x_i+(x-x_i)$ and the hat is symmetric about $x_i=ih$, its load is $F_i=x_i\int {\phi }_i=x_{i}h=i{h}^2$. Thus the reaction-diffusion finite element matrix gives the explicit equationsThe coefficient matrix is a tridiagonal matrix and a positive-definite matrix by the energy identity in part (b), hence nonsingular. A separate spectral verification uses the sine basis with ${\theta }_j=j\pi /(M+1)$ and givesBoth terms are positive for the Dirichlet modes. This also identifies precisely how the diffusion stiffness matrix and reaction mass matrix combine.

The starting point of a conforming finite element method is a weak formulation in a Hilbert space $X$. For example, a homogeneous Dirichlet boundary condition on a bounded Lipschitz domain gives $X=H_0^1(\Omega )$. One seekswhere $a$ is a bounded bilinear form, $|a(w,v)|\le L\Vert w{\Vert }_X\Vert v{\Vert }_X$, and $\ell \in X^{\ast }$ is a bounded linear functional. A conforming finite element space $X_h\subset X$ is finite dimensional and is typically built from functions that are polynomial on each element of a finite element mesh. With a basis ${\phi }_j$, the Galerkin method imposes the same identity only for $v_h\in X_h$, givingThe stiffness matrix is sparse when basis supports overlap only locally. Local element integrals assemble its entries; the mass matrix represents an $L^2$ term. This separates the choice of trial functions from the variational principle determining their coefficients.

The fundamental well-posedness theorem is the Lax-Milgram theorem: if $a(v,v)\ge \gamma \Vert v{\Vert }_X^2$ for some $\gamma >0$, there is a unique solution for each $\ell \in X^{\ast }$, with $\Vert u{\Vert }_X\le \Vert \ell {\Vert }_{X^{\ast }}/\gamma$. The same theorem on $X_h$ gives a unique discrete solution. The constants must be independent of the mesh for uniform stability of a numerical method. Boundary conditions are built into the trial space when they are essential constraints; flux conditions arise through boundary terms in integration by parts and are natural constraints in the weak formulation.

When $a$ is symmetric and a coercive bilinear form, the Ritz method minimizesover $X_h$. Its first variation is exactly the Galerkin method, so Ritz minimization and Galerkin testing coincide for symmetric coercive problems. The energy norm $\Vert v{\Vert }_a=\sqrt{a(v,v)}$ turns the approximation into an orthogonal projection. Indeed, subtracting the exact and discrete weak equations gives Galerkin orthogonality, $a(u-u_h,v_h)=0$. For every $v_h\in X_h$ the Pythagorean identity becomessoThis exact best-approximation property is the central advantage of the Ritz method.

Symmetry is unnecessary for the Galerkin method. For any bounded, coercive form, the Céa lemma givesTo prove it, put $e=u-u_h$. Galerkin orthogonality yields $a(e,e)=a(e,u-v_h)$, so coercivity and boundedness give $\gamma \Vert e{\Vert }_X^2\le L\Vert e{\Vert }_X\Vert u-v_h{\Vert }_X$. Divide when $e\ne 0$ and take the infimum. Thus approximation capability plus uniform coercivity gives convergence. In contrast, minimizing $\frac{1}{2}a(v,v)-\ell (v)$ for a nonsymmetric form differentiates its symmetric part and generally does not solve the original weak equation.

A concrete symmetric example is the one-dimensional reaction-diffusion problem from Question 5. Continuous piecewise-linear hat functions give the diffusion stiffness matrix $S$ with diagonal $2/h$ and neighbouring entries $-1/h$, and the mass matrix $M$ with diagonal $2h/3$ and neighbouring entries $h/6$. Therefore $K=S+M$ and $F_i=i{h}^2$. Its positive-definite matrix is the coordinate form of the continuous energy. The exact smooth solution $u(x)=x-\sinh (x)/\sinh 1$ makes this an explicit test of the method, not just an abstract existence result.

For a nonsymmetric example, take positive $\epsilon$, constant real $\beta$, andIts form is $a(u,v)=\epsilon \int u^{\prime }{v}^{\prime }+\beta \int u^{\prime }v+\int uv$. Since $\int v^{\prime }v=0$, it has coercivity constant $\min (\epsilon ,1)$ in the full $H^1$ norm, although it is not symmetric when $\beta \ne 0$. The Galerkin method remains well posed. On the uniform hat basis, $K=\epsilon S+\beta C+M$, where $C_{i,i+1}=1/2$, $C_{i,i-1}=-1/2$ and $C_{ii}=0$. The skew-symmetric matrix $C$ cancels from $c^{T}Kc$, proving nonsingularity. A Ritz energy would omit precisely this advection contribution. Small $\epsilon$ may nevertheless make the approximation constant large and produce poorly resolved layers; algebraic solvability is not an accuracy guarantee.

Quantitative convergence uses a finite element interpolation estimate. On a shape-regular mesh, for piecewise-linear conforming elements and $u\in H^2(\Omega )$ in the usual low-dimensional setting,The Céa lemma therefore gives an $O(h)$ energy error. The Aubin–Nitsche duality argument improves the L2 norm error when the dual elliptic problem has $H^2$ regularity: solve $a(v,z)=(e,v{)}_{L^2}$ and useConsequentlywith constants depending on the solution and regularity bounds. The second estimate is conditional on dual regularity; corners or insufficient data regularity can reduce the rate.

The same Galerkin method handles evolution equations. For the homogeneous heat equation, the semidiscrete finite element heat equation is $M\dot{c}+Sc=0$. Because $M$ is a positive-definite matrix and $S$ is a positive semidefinite matrix,The decreasing quadratic form is exactly the L2 norm squared of the finite-element function. This separates spatial stability of a numerical method from the subsequent time integrator, just as in the Schrödinger equation example. Conformity, coercivity and approximation estimates together explain convergence.

An $s$-stage Runge-Kutta method for $y^{\prime }=f(t,y)$ is defined byThe Butcher tableau records $A=(a_{ij})$, $b$ and $c$; internal consistency usually sets $c=Ae$, and first-order consistency requires $b^{T}e=1$. A strictly lower triangular $A$ gives an explicit method, while other stage dependencies generally require an implicit Runge-Kutta method. The stage solvability of an implicit Runge-Kutta method must be checked separately: if $f$ is globally Lipschitz continuous in $y$ with constant $L$, the stage fixed-point map is a contraction in the maximum stage norm whenever $hL{\max }_i{\sum }_j|a_{ij}|<1$. This is a sufficient small-step condition, not a restriction intrinsic to the definitions of A-stability or B-stability.

Linear stability begins with the Dahlquist test equation $y^{\prime }=\lambda y$. Solving the stages givesThe linear stability domain consists of test values for which the stage system is well defined and $|R(z)|\le 1$. A-stability means this includes the closed left half-plane. For $\mathrm{Re}\lambda <0$, this gives stability with no scalar decay-mode step restriction. For a normal matrix $L$ in $y^{\prime }=Ly$, unitary diagonalization of a normal matrix reduces the norm estimate to these scalar factors. For a non-normal matrix, eigenvalues alone do not establish a uniform norm bound: eigenvector conditioning and transient amplification matter. Mesh-uniform estimates for discretized PDEs must control operator powers, not only their spectra.

Every consistent explicit Runge-Kutta method has a nonconstant polynomial stability function and therefore cannot be A-stable, since that polynomial is unbounded on the negative real axis. For example, forward Euler method has $R(z)=1+z$, stability disk $|1+z|\le 1$, and restriction $0\le hq\le 2$ for $\lambda =-q<0$. For a pure imaginary mode $\lambda =i\omega \ne 0$, $|1+ih\omega |>1$, so it is unstable for every positive step. The classical fourth-order Runge-Kutta method instead hasIt is stable on the imaginary axis exactly for $|q|\le 2\sqrt{2}$, an illustrative conditional stability range for oscillatory evolution.

Backward Euler method has $R(z)=(1-z{)}^{-1}$ and is A-stable. L-stability adds $R(z)\to 0$ as $|z|\to \infty$ in the left half-plane; backward Euler satisfies this and strongly damps unresolved rapidly decaying modes in a stiff differential equation. The implicit midpoint rule and the trapezoidal rule both have $R(z)=(1+z/2)/(1-z/2)$, so they are A-stable but not L-stable. They preserve the modulus of a pure imaginary test mode, but their stiff-decay limit is $-1$, leaving oscillatory numerical remnants. The three-stage Lobatto IIIA method has the fourth-order rational function found in Question 2 and likewise lacks L-stability. Thus high order and A-stability do not themselves imply efficient stiff damping.

Nonlinear stability measures differences between solutions of the same equation. A one-sided Lipschitz condition isDifferentiating the squared difference and applying the Gronwall inequality gives $\Vert y(t)-\tilde{y}(t)\Vert \le e^{\nu (t-t_0)}\Vert y(t_0)-\tilde{y}(t_0)\Vert$. A dissipative vector field has $\nu \le 0$ and is contractive. A B-stable method preserves this property for every positive step size for which the stages are well defined:This includes time-dependent vector fields when dissipativity holds at each common time argument. Applying it to scalar linear dissipative fields shows that B-stability implies A-stability when the linear stages are well defined; the converse fails.

The standard sufficient theorem is algebraic stability implies B-stability, subject to stage solvability. Define $B=\mathrm{diag}(b)$ and $M=BA+A^{T}B-b{b}^T$. Algebraic stability of a Runge-Kutta method means $b_i\ge 0$ and $M$ is a positive semidefinite matrix. To prove the theorem, set $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$ and $F_i=f(t_n+c_{i}h,Y_i)-f(t_n+c_{i}h,{\tilde{Y}}_i)$. Expanding the output difference and using $D_i=d+h{\sum }_j{a}_{ij}{F}_j$ yields the Runge-Kutta contractivity identityDissipativity makes the first sum nonpositive, while positive semidefiniteness makes the last quadratic form nonnegative. The latter follows by factoring $M=Q^{T}Q$ and writing it as a sum of squared norms of linear combinations of the $F_i$. This proves the claimed contraction. It is a sufficient theorem; absence of algebraic stability is not, by itself, a proof of failure of B-stability for every representation.

Several examples make the distinctions concrete. For backward Euler, $b=1$, $A=(1)$, so $M=(1)$: it is both L-stable and B-stable. For implicit midpoint, $b=1$, $A=(1/2)$, so $M=(0)$: it is B-stable and preserves squared norms for a skew-Hermitian linear equation, but lacks stiff damping. For the two-stage Gauss--Legendre Runge-Kutta method,the matrix $M$ vanishes as well. The collocation theorem gives order four, and its stability function is exactly the same as the three-stage Lobatto IIIA function. Yet the Gauss method is algebraically stable, whereas that Lobatto method is not. Identical scalar linear stability functions need not imply identical nonlinear stability properties.

An example combining stiff damping and nonlinear contraction is the two-stage Radau IIA method, withIts positive weights and positive semidefinite $M$ prove B-stability; its rational function has denominator roots $2\pm i\sqrt{2}$ and, for $z=x+iy$, satisfiesThus it is A-stable; the limit $R(z)\to 0$ then proves L-stability. The standard Radau IIA collocation theorem gives order $2s-1$, here three. These properties explain its usefulness for stiff nonlinear equations.

Finally, A-stability alone does not imply B-stability. The trapezoidal rule fails B-stability even for the scalar dissipative equation $y^{\prime }=-y^3$. With $h=1$, its step map $T$ is uniquely defined byAt $y=\sqrt{2}$, $T(y)=0$. Implicit differentiation givesThus nearby starting values are separated by approximately twice their original distance after one step, although the exact flow is contractive. This example and the contractivity theorem show why scalar linear analysis, stiff damping and genuinely nonlinear norm estimates answer different stability questions.