Numerical analysis

On 25 February 1991 a Patriot missile battery at Dhahran failed to intercept an incoming missile, and 28 soldiers died. The system counted time in tenths of a second and converted to seconds by multiplying by $0.1$, stored in a 24-bit register. But $0.1$ has no finite binary expansion, and the stored value was too small by about $9.5 \times 10^{-8}$. After 100 hours of continuous operation the clock had drifted by $0.34$ seconds, during which the target travelled more than half a kilometre (US General Accounting Office report, 1992).

Rounding errors are usually harmless. Every weather forecast, every bridge designed by finite elements and every GPS fix solves enormous systems of equations in floating point, and the answers are excellent. Numerical analysis is the mathematics that explains when and why, and it rests on three questions.

• Is the problem sensitive? Some problems amplify any small change in their data. That is a property of the problem — its condition number — and no algorithm can repair it.
• Is the algorithm stable? A good algorithm adds no more error than the inevitable rounding of the data. Gaussian elimination with pivoting is good; without pivoting it can fail on an innocent $2 \times 2$ system.
• How fast does the method converge? Discretising a differential equation with step $h$ gives an error like $h^p$, where $p$ is the order, provided the method is stable — and for "stiff" problems, stability can matter far more than order.

This lesson develops these ideas for the two workhorses of scientific computing: linear algebra (linear systems, least squares and eigenvalues) and ordinary differential equations.

Think of a problem as a function $f$ from data to answers, and of an algorithm as a computed approximation $\tilde f$. The forward error $|\tilde f(x) - f(x)|$ is what we care about but usually cannot measure. The backward error asks a cleverer question: for which nearby data $x + \delta x$ is the computed answer exact? If $\tilde f(x) = f(x + \delta x)$ with $\delta x$ as small as the rounding of the data, the algorithm is backward stable: it has answered a question we could not distinguish from ours.

Then the error splits into two factors:

The algorithm controls the second factor and the problem controls the first. Geometrically, a matrix $A$ maps the unit sphere to an ellipsoid. If the ellipsoid is very thin, two very different inputs have nearly the same image, and inverting $A$ magnifies errors by the ratio of its longest to its shortest axis. That ratio is the condition number.

For differential equations the analogue of backward stability is stability of the scheme: errors made in early steps must not be amplified by later ones. A method that is consistent (each step is accurate to some order) and stable converges. Stiff problems contain a rapidly decaying component that is invisible in the solution, but an explicit method must still take tiny steps to keep that component from exploding.

Since rounding the data alone causes relative perturbations of size $u$, no algorithm can promise relative accuracy better than about $\kappa(A)\,u$. The rule of thumb: with $\kappa \approx 10^k$, expect to lose about $k$ of the $16$ significant decimal digits.

Gaussian elimination subtracts multiples $\ell_{ik} = a^{(k)}_{ik}/a^{(k)}_{kk}$ of the pivot row $k$ from the rows below it. Each step is multiplication by a unit lower triangular matrix $L_k$, and the product of their inverses is the unit lower triangular matrix $L$ of multipliers. So elimination computes $A = LU$. Solving $Ax = b$ then takes a forward substitution $Ly = b$ and a back substitution $Ux = y$.

The matrix $\begin{pmatrix}0 & 1\\1 & 1\end{pmatrix}$ is perfectly well conditioned but has no LU factorisation, since $A_1 = (0)$. A nearly zero pivot is just as bad in floating point. For $\begin{pmatrix}\varepsilon & 1\\1 & 1\end{pmatrix}$ the multiplier is $1/\varepsilon$ and $u_{22} = 1 - 1/\varepsilon$. With $\varepsilon = 10^{-20}$, $\mathrm{fl}(1 - 10^{20}) = -10^{20}$, and the computed solution of $Ax = (1, 2)^T$ is $(0, 1)$ instead of $(1, 1)$. Swapping the rows first gives multiplier $\varepsilon$ and an accurate answer. Partial pivoting always selects the largest available pivot. Then every $|\ell_{ik}| \le 1$, and every invertible $A$ has a factorisation $PA = LU$: at each step the pivot column below the diagonal contains a non-zero entry, since otherwise the matrix would be singular.

(Wilkinson, 1961; see Higham, Accuracy and Stability of Numerical Algorithms (2002), chapter 9.) The bound $2^{n-1}$ is attained, so partial pivoting is not backward stable in the worst case. In practice growth factors are small — of order $n^{2/3}$ on average for random matrices, according to experiments of Trefethen and Schreiber (1990) — and matrices with large growth seem almost never to arise. Why this is so is still not fully understood. It is an empirical fact supported by probabilistic analysis, not a theorem about all matrices. LU costs about $\frac23n^3$ flops; each additional right-hand side costs only $2n^2$.

Do not compute $A^{-1}$ to solve $Ax = b$. Forming $A^{-1}$ costs about $2n^3$ flops, three times as much as LU. The computed $\mathrm{fl}(A^{-1})b$ is also not backward stable: its residual $\|b - A\hat x\|$ can be about $\kappa(A)$ times larger than the residual of the LU solution, although both forward errors are of size $\kappa(A)u$ (Higham 2002, §14.1). And for sparse $A$, the inverse is typically dense while the LU factors often stay sparse.

Given $A \in \R^{m\times n}$ with $m > n$ and full column rank, and $b \in \R^m$, the least-squares problem asks for $x$ minimising $\|Ax - b\|_2$ — for example, fitting a line to $m$ noisy measurements.

So solving the normal equations in floating point loses twice as many digits as the problem deserves. In Läuchli's example $A = \begin{pmatrix}1 & 1\\ \varepsilon & 0\\ 0 & \varepsilon\end{pmatrix}$ with $\varepsilon = 10^{-8}$, $A^TA = \begin{pmatrix}1 + \varepsilon^2 & 1\\ 1 & 1 + \varepsilon^2\end{pmatrix}$, and $\mathrm{fl}(1 + 10^{-16}) = 1$ makes it exactly singular. QR computed by Householder reflections works on $A$ itself and is backward stable (Higham 2002, chapter 19). Orthogonal transformations do not amplify errors, since $\kappa_2(Q) = 1$.

Eigenvalues by the QR algorithm. Starting from $A_0 = A$, factor $A_k = Q_kR_k$ and set $A_{k+1} = R_kQ_k$. Since $A_{k+1} = Q_k^TA_kQ_k$, every $A_k$ is orthogonally similar to $A$ and has the same eigenvalues. Theorem (stated): if the eigenvalues of $A$ have distinct absolute values $|\lambda_1| > |\lambda_2| > \cdots > |\lambda_n|$ (plus a mild condition on the eigenvectors), then $A_k$ converges to an upper triangular matrix with the $\lambda_i$ on the diagonal. The subdiagonal entry $(i+1, i)$ decays like $|\lambda_{i+1}/\lambda_i|^k$. The reason is that $Q_0Q_1\cdots Q_{k-1}$ is the orthonormalisation of $A^k$ — simultaneous power iteration on all columns at once. Practical versions first reduce $A$ to Hessenberg form (cost $O(n^3)$ once, then $O(n^2)$ per step) and use shifts $A_k - \mu_kI$, which make convergence typically quadratic, or cubic for symmetric matrices (Trefethen and Bau 1997, lectures 28–29).

Consider $y' = f(t, y)$, $y(0) = y_0$ on $[0, T]$, with $f$ continuous and $L$-Lipschitz in $y$: $|f(t, u) - f(t, v)| \le L|u - v|$. Take steps $h = T/N$, $t_n = nh$.

The proof has the shape of every convergence proof in numerical analysis: consistency (small local errors) plus stability (errors are not amplified too much) gives convergence. The whole argument is formalised in Lean below.

Runge–Kutta methods evaluate $f$ several times per step to cancel more terms of the Taylor expansion. Heun's method $y_{n+1} = y_n + \frac h2\big(k_1 + k_2\big)$, with $k_1 = f(t_n, y_n)$ and $k_2 = f(t_n + h, y_n + hk_1)$, has order 2. The classical fourth-order method of Kutta (1901),

has order 4. Theorem (stated): a Runge–Kutta method whose local truncation error is $O(h^{p+1})$ for smooth $f$ (equivalently, whose coefficients satisfy Butcher's order conditions up to order $p$) has global error $O(h^p)$, by the same consistency-plus-stability argument as Theorem 6. Checking that the classical method satisfies the eight order conditions for $p = 4$ is a finite computation (Hairer, Nørsett and Wanner, Solving Ordinary Differential Equations I, §II.2). The experiment below measures the orders $1, 2, 4$.

Absolute stability. On the test equation $y' = \lambda y$ with $\operatorname{Re}\lambda < 0$, whose solutions decay, we want the numerical solution $y_n = R(h\lambda)^ny_0$ to decay too, i.e. $|R(h\lambda)| < 1$. The stability functions are:

A problem is stiff when it contains decay rates far faster than the time scale of interest. The equation $y' = -50(y - \cos t)$ has the smooth solution $y \approx \cos t$ after a transient of length $\approx 0.1$, yet explicit Euler needs $50h < 2$ for ever after. In chemical kinetics or circuit simulation, rates can differ by factors of $10^{10}$, and only implicit methods are practical. Germund Dahlquist (1963) proved the second Dahlquist barrier: an A-stable linear multistep method has order at most $2$. Among the A-stable multistep methods of order 2, the trapezoidal rule has the smallest error constant.

Smaller is not always better. Each step makes a rounding error of order $u$, and $N = T/h$ steps accumulate roughly $uT/h$ of it. So the total error behaves like $Ch^p + uT/h$, which is smallest at a moderate $h$ and grows as $h \to 0$. The same trade-off appears in a single difference quotient: $(f(x + h) - f(x))/h$ has truncation error $\approx h|f''|/2$ and rounding error $\approx 2u|f|/h$, optimal near $h \approx \sqrt u \approx 10^{-8}$.

Mathlib has operator norms, matrices and determinants, but no model of floating-point arithmetic and no library of numerical algorithms. So what can be formalised is the exact mathematics behind the error analysis: perturbation bounds, factorisations, and the convergence and stability estimates. The rounding-error statements (Theorem 4, backward stability of Householder QR) are not formalised.

The condition number of an invertible bounded linear map between normed spaces, and Theorem 1 in multiplied-out form (which avoids dividing by $\|x\|$ and $\|b\|$):

/-- The **condition number** `κ(A) = ‖A‖ ‖A⁻¹‖` of an invertible bounded linear map between normed
spaces (for matrices, with the operator norm induced by the chosen vector norm). -/
noncomputable def condNum {E F : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E]
    [NormedAddCommGroup F] [NormedSpace ℝ F] (A : E ≃L[ℝ] F) : ℝ :=
  ‖(A : E →L[ℝ] F)‖ * ‖(A.symm : F →L[ℝ] E)‖

/-- The defining property of the operator norm: `‖A x‖ ≤ ‖A‖ ‖x‖`. This applies Mathlib's
`ContinuousLinearMap.le_opNorm`. -/
theorem op_norm_bound {E F : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E]
    [NormedAddCommGroup F] [NormedSpace ℝ F] (A : E →L[ℝ] F) (x : E) : ‖A x‖ ≤ ‖A‖ * ‖x‖ :=
  A.le_opNorm x

/-- **`κ(A) ≥ 1`.** For `x ≠ 0`, `‖x‖ = ‖A⁻¹(A x)‖ ≤ ‖A⁻¹‖ ‖A‖ ‖x‖`; divide by `‖x‖`. A condition
number near `1` (an orthogonal matrix in the 2-norm) is as good as it gets. -/
theorem one_le_condNum {E F : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E]
    [NormedAddCommGroup F] [NormedSpace ℝ F] (A : E ≃L[ℝ] F) {x : E} (hx : x ≠ 0) :
    1 ≤ condNum A := by
  have hxpos : 0 < ‖x‖ := norm_pos_iff.mpr hx
  have h1 : ‖x‖ ≤ ‖(A.symm : F →L[ℝ] E)‖ * ‖A x‖ := by
    have := (A.symm : F →L[ℝ] E).le_opNorm (A x)
    simpa using this
  have h2 : ‖A x‖ ≤ ‖(A : E →L[ℝ] F)‖ * ‖x‖ := (A : E →L[ℝ] F).le_opNorm x
  have h3 : ‖x‖ ≤ condNum A * ‖x‖ := by
    calc ‖x‖ ≤ ‖(A.symm : F →L[ℝ] E)‖ * ‖A x‖ := h1
      _ ≤ ‖(A.symm : F →L[ℝ] E)‖ * (‖(A : E →L[ℝ] F)‖ * ‖x‖) :=
          mul_le_mul_of_nonneg_left h2 (norm_nonneg _)
      _ = condNum A * ‖x‖ := by unfold condNum; ring
  nlinarith

/-- **The perturbation bound.** If `A x = b` and `A (x + δx) = b + δb`, then
`‖δx‖ · ‖b‖ ≤ κ(A) · ‖δb‖ · ‖x‖`, i.e. `‖δx‖/‖x‖ ≤ κ(A) ‖δb‖/‖b‖`: the relative error in the
solution is at most the condition number times the relative error in the data. Proof: by
linearity `A δx = δb`, so `δx = A⁻¹ δb` and `‖δx‖ ≤ ‖A⁻¹‖ ‖δb‖`; and `‖b‖ = ‖A x‖ ≤ ‖A‖ ‖x‖`.
Multiply the two inequalities. -/
theorem relative_error_bound {E F : Type*} [NormedAddCommGroup E] [NormedSpace ℝ E]
    [NormedAddCommGroup F] [NormedSpace ℝ F] (A : E ≃L[ℝ] F) {x δx : E} {b δb : F}
    (hb : A x = b) (hpert : A (x + δx) = b + δb) :
    ‖δx‖ * ‖b‖ ≤ condNum A * ‖δb‖ * ‖x‖ := by
  have hAδ : A δx = δb := by
    rw [map_add, hb] at hpert
    exact add_left_cancel hpert
  have h1 : ‖δx‖ ≤ ‖(A.symm : F →L[ℝ] E)‖ * ‖δb‖ := by
    have := (A.symm : F →L[ℝ] E).le_opNorm δb
    rw [← hAδ] at this ⊢
    simpa using this
  have h2 : ‖b‖ ≤ ‖(A : E →L[ℝ] F)‖ * ‖x‖ := by
    rw [← hb]; exact (A : E →L[ℝ] F).le_opNorm x
  calc ‖δx‖ * ‖b‖ ≤ (‖(A.symm : F →L[ℝ] E)‖ * ‖δb‖) * (‖(A : E →L[ℝ] F)‖ * ‖x‖) :=
        mul_le_mul h1 h2 (norm_nonneg _) (by positivity)
    _ = condNum A * ‖δb‖ * ‖x‖ := by unfold condNum; ring

An LU factorisation produced by elimination, and the determinant as the product of the pivots:

/-- **An LU factorisation**, as produced by Gaussian elimination: subtract `2 ×` row 1 from row 2
and `4 ×` row 1 from row 3, then `3 ×` the new row 2 from row 3. The multipliers `2, 4, 3` fill
the unit lower triangular factor `L`; the result of elimination is the upper triangular `U`.
Checked by exact rational arithmetic. -/
theorem lu_example :
    !![(2 : ℚ), 1, 1; 4, 3, 3; 8, 7, 9] = !![(1 : ℚ), 0, 0; 2, 1, 0; 4, 3, 1] * !![(2 : ℚ), 1, 1; 0, 1, 1; 0, 0, 2] := by
  ext i j
  fin_cases i <;> fin_cases j <;> simp [Matrix.mul_apply, Fin.sum_univ_succ] <;> norm_num

/-- `det (L U) = det L · det U`. This applies Mathlib's `Matrix.det_mul`. With `L` unit lower
triangular (`det L = 1`) and `U` upper triangular, `det A` is the product of the pivots — the
cheapest way to compute a determinant (about `2n³/3` operations instead of `n!` terms). -/
theorem det_mul_lu {n : ℕ} (L U : Matrix (Fin n) (Fin n) ℚ) : (L * U).det = L.det * U.det :=
  Matrix.det_mul L U

/-- The determinant of the example matrix is the product of the pivots `2 · 1 · 2 = 4`, computed
through the factorisation `lu_example` (`det L = 1` and `det U = 4`, each by the formula for a
3 × 3 determinant). -/
theorem det_via_lu : (!![(2 : ℚ), 1, 1; 4, 3, 3; 8, 7, 9]).det = 4 := by
  rw [lu_example, Matrix.det_mul, Matrix.det_fin_three, Matrix.det_fin_three]
  simp
  norm_num

The pivoting example, exactly, with and without the row swap:

/-- **Why pivot: elimination without a row swap.** For `ε ≠ 0`,
`[[ε, 1], [1, 1]] = [[1, 0], [1/ε, 1]] · [[ε, 1], [0, 1 - 1/ε]]`. When `ε` is tiny the multiplier
`1/ε` and the entry `1 - 1/ε` are huge, and in floating point `1 - 1/ε` rounds to `-1/ε`,
destroying the information in the `(2, 2)` entry: the computed solution can be completely wrong
although the matrix is perfectly well conditioned. -/
theorem no_pivot_lu {ε : ℝ} (hε : ε ≠ 0) :
    !![ε, 1; 1, 1] = !![1, 0; 1 / ε, 1] * !![ε, 1; 0, 1 - 1 / ε] := by
  ext i j
  fin_cases i <;> fin_cases j <;> simp [Matrix.mul_apply, Fin.sum_univ_succ]
  field_simp

/-- **With partial pivoting.** Swapping the rows first (`P A` with `P` the swap) gives
`[[1, 1], [ε, 1]] = [[1, 0], [ε, 1]] · [[1, 1], [0, 1 - ε]]`: the multiplier is `ε`, at most `1`
in absolute value, and no entry grows. Partial pivoting always chooses the largest available
pivot, so every multiplier satisfies `|ℓᵢⱼ| ≤ 1`. -/
theorem pivot_lu (ε : ℝ) :
    !![1, 1; ε, 1] = !![1, 0; ε, 1] * !![(1 : ℝ), 1; 0, 1 - ε] := by
  ext i j
  fin_cases i <;> fin_cases j <;> simp [Matrix.mul_apply, Fin.sum_univ_succ]

Theorem 6, Euler's method has order 1, in three steps that follow the human proof — the error recursion, the discrete Gronwall inequality, and the global bound (with $1 + x \le e^x$ from Mathlib):

/-- **One step of the error recursion.** Let `yₙ₊₁ = yₙ + h f(tₙ, yₙ)` (Euler) and let the exact
solution satisfy `y(tₙ₊₁) = y(tₙ) + h f(tₙ, y(tₙ)) + τₙ` (`τₙ` is the local truncation error).
If `f` is `L`-Lipschitz in `y`, the error `eₙ = y(tₙ) - yₙ` satisfies
`|eₙ₊₁| ≤ (1 + h L) |eₙ| + |τₙ|`. -/
theorem euler_error_step {f : ℝ → ℝ → ℝ} {L h t yn Y τ : ℝ} (hh : 0 ≤ h)
    (hlip : ∀ s u v, |f s u - f s v| ≤ L * |u - v|) :
    |(Y + h * f t Y + τ) - (yn + h * f t yn)| ≤ (1 + h * L) * |Y - yn| + |τ| := by
  have hf := hlip t Y yn
  calc |(Y + h * f t Y + τ) - (yn + h * f t yn)|
      = |(Y - yn) + h * (f t Y - f t yn) + τ| := by ring_nf
    _ ≤ |Y - yn| + |h * (f t Y - f t yn)| + |τ| := abs_add_three _ _ _
    _ = |Y - yn| + h * |f t Y - f t yn| + |τ| := by rw [abs_mul, abs_of_nonneg hh]
    _ ≤ |Y - yn| + h * (L * |Y - yn|) + |τ| := by gcongr
    _ = (1 + h * L) * |Y - yn| + |τ| := by ring

/-- **Discrete Gronwall inequality.** If `eₙ₊₁ ≤ (1 + h L) eₙ + τ` with `h L ≥ 0` and `τ ≥ 0`,
then `eₙ ≤ (1 + h L)ⁿ e₀ + n τ (1 + h L)ⁿ`. Errors committed at each step are amplified by at most
`(1 + h L)` per subsequent step. Induction on `n`. -/
theorem discrete_gronwall {e : ℕ → ℝ} {h L τ : ℝ} (hhL : 0 ≤ h * L) (hτ : 0 ≤ τ)
    (hstep : ∀ n, e (n + 1) ≤ (1 + h * L) * e n + τ) (n : ℕ) :
    e n ≤ (1 + h * L) ^ n * e 0 + n * τ * (1 + h * L) ^ n := by
  have hq : 1 ≤ 1 + h * L := by linarith
  induction n with
  | zero => simp
  | succ n ih =>
    have hpow : 1 ≤ (1 + h * L) ^ (n + 1) := one_le_pow₀ hq
    calc e (n + 1) ≤ (1 + h * L) * e n + τ := hstep n
      _ ≤ (1 + h * L) * ((1 + h * L) ^ n * e 0 + n * τ * (1 + h * L) ^ n) + τ := by
          gcongr
      _ = (1 + h * L) ^ (n + 1) * e 0 + n * τ * (1 + h * L) ^ (n + 1) + τ := by ring
      _ ≤ (1 + h * L) ^ (n + 1) * e 0 + n * τ * (1 + h * L) ^ (n + 1)
            + τ * (1 + h * L) ^ (n + 1) := by nlinarith
      _ = (1 + h * L) ^ (n + 1) * e 0 + ((n + 1 : ℕ) : ℝ) * τ * (1 + h * L) ^ (n + 1) := by
          push_cast; ring

/-- `1 + x ≤ eˣ` for every real `x`, so `(1 + h L)ⁿ ≤ e^{n h L}`. This applies Mathlib's
`Real.add_one_le_exp`. -/
theorem one_add_le_exp (x : ℝ) : 1 + x ≤ Real.exp x := by
  rw [add_comm]; exact Real.add_one_le_exp x

/-- **Euler's method has order 1.** Suppose the errors satisfy `e₀ = 0` and
`eₙ₊₁ ≤ (1 + h L) eₙ + C h²` (the recursion of `euler_error_step`, with local truncation error
`|τₙ| ≤ C h²`, e.g. `C = max |y''| / 2`). Then at time `T = n h`,
`eₙ ≤ C h · T · e^{L T}`. For fixed `T`, the error is `O(h)`: halving the step halves the error. -/
theorem euler_global_error {e : ℕ → ℝ} {h L C : ℝ} (hh : 0 ≤ h) (hL : 0 ≤ L) (hC : 0 ≤ C)
    (h0 : e 0 = 0) (hstep : ∀ n, e (n + 1) ≤ (1 + h * L) * e n + C * h ^ 2) (n : ℕ) :
    e n ≤ C * h * (n * h) * Real.exp (L * (n * h)) := by
  have hhL : 0 ≤ h * L := mul_nonneg hh hL
  have hg := discrete_gronwall hhL (by positivity) hstep n
  rw [h0, mul_zero, zero_add] at hg
  have hpow : (1 + h * L) ^ n ≤ Real.exp (L * (n * h)) := by
    calc (1 + h * L) ^ n ≤ Real.exp (h * L) ^ n :=
          pow_le_pow_left₀ (by linarith) (one_add_le_exp _) n
      _ = Real.exp (L * (n * h)) := by rw [← Real.exp_nat_mul]; ring_nf
  calc e n ≤ n * (C * h ^ 2) * (1 + h * L) ^ n := hg
    _ ≤ n * (C * h ^ 2) * Real.exp (L * (n * h)) := by
        apply mul_le_mul_of_nonneg_left hpow; positivity
    _ = C * h * (n * h) * Real.exp (L * (n * h)) := by ring

Stability on the test equation: explicit Euler multiplies by $1 + h\lambda$ and is stable exactly when $h|\lambda| < 2$, while backward Euler contracts for every step size:

/-- **Explicit Euler on `y' = λ y`** is `yₙ₊₁ = yₙ + h λ yₙ = (1 + hλ) yₙ`, so after `n` steps
`yₙ = (1 + hλ)ⁿ y₀`. The exact solution is `e^{λ n h} y₀`. -/
theorem explicit_euler_iterate (h lam y₀ : ℝ) (n : ℕ) :
    (fun y => y + h * (lam * y))^[n] y₀ = (1 + h * lam) ^ n * y₀ := by
  induction n with
  | zero => simp
  | succ n ih => rw [Function.iterate_succ_apply', ih]; ring

/-- **Stability of explicit Euler.** For a decaying equation (`λ < 0`) and step `h > 0`, the
numerical solution decays (`|1 + hλ| < 1`) exactly when `h |λ| < 2`. For a stiff equation with
`λ = -1000` this forces `h < 0.002`, whatever accuracy is wanted. -/
theorem explicit_euler_stable_iff {h lam : ℝ} (hh : 0 < h) (hlam : lam < 0) :
    |1 + h * lam| < 1 ↔ h * (-lam) < 2 := by
  have : 0 < h * (-lam) := mul_pos hh (by linarith)
  rw [abs_lt]
  constructor
  · rintro ⟨h1, -⟩; linarith
  · intro h2; constructor <;> nlinarith

/-- **Backward Euler is unconditionally stable on decaying modes.** The implicit step
`yₙ₊₁ = yₙ + h λ yₙ₊₁` gives `yₙ₊₁ = yₙ / (1 - hλ)`, and for `λ < 0` and every `h > 0` the factor
`1/(1 - hλ)` lies strictly between `0` and `1`: the numerical solution decays, like the true
one, with no restriction on the step size (A-stability, on the negative real axis). -/
theorem backward_euler_contracts {h lam : ℝ} (hh : 0 < h) (hlam : lam < 0) :
    0 < 1 / (1 - h * lam) ∧ 1 / (1 - h * lam) < 1 := by
  have hpos : 1 < 1 - h * lam := by nlinarith
  constructor
  · positivity
  · rw [div_lt_one (by linarith)]; exact hpos

What is and is not formalised. Proved here: $\kappa \ge 1$, Theorem 1, the example factorisations, the Euler error recursion, discrete Gronwall and Theorem 6 (given the one-step recursion), and Theorem 7 on the negative real axis for the two Euler methods. Cited from Mathlib: the operator-norm inequality, $\det(LU) = \det L\det U$ and $1 + x \le e^x$. Not formalised: Theorem 2, Theorems 3–5 in general, the QR algorithm, the order of RK4, and anything about floating-point rounding.

The first program puts the linear-algebra theorems to work in IEEE double precision. Hilbert systems are solved with relative error close to $\kappa_2(H)u$, as Theorems 1–2 predict for a backward-stable solver. (The right-hand side $b = H\mathbf{1}$ is itself rounded, which already causes errors of this size.) Elimination without pivoting fails on the $\varepsilon$-matrix. Wilkinson's matrix achieves the worst-case growth $2^{n-1}$ and destroys all accuracy at $n = 55$, although its condition number is modest. And the normal equations for Läuchli's problem are exactly singular in floating point, while QR solves it.

"""Conditioning, pivoting and least squares in floating point.

Four numerical experiments in IEEE double precision (unit roundoff
u = 2^-53, about 1.1e-16). (1) Hilbert matrices H_ij = 1/(i + j - 1) are
notoriously ill conditioned: solving H x = b with a backward-stable solver
gives a relative error of roughly kappa(H) u, as the perturbation bound
predicts. (2) Gaussian elimination without pivoting on [[1e-20, 1], [1, 1]]
returns a wrong answer although the matrix is well conditioned; with
partial pivoting the answer is right. (3) Wilkinson's matrix makes partial
pivoting's growth factor equal 2^(n-1), the worst case. (4) For Lauchli's
least-squares problem the normal equations A^T A x = A^T b are exactly
singular in floating point, while a QR factorisation solves it accurately.

Evidence: numerical (floating point), compared with exact solutions.
"""

import numpy as np

UNIT_ROUNDOFF = 2.0 ** -53

def hilbert(n):
    """The n x n Hilbert matrix."""
    i = np.arange(1, n + 1)
    return 1.0 / (i[:, None] + i[None, :] - 1)

def lu_solve(a, b, pivot=True):
    """Solve a x = b by Gaussian elimination (optionally partial pivoting).

    Returns (x, growth) where growth = max |U_ij| / max |A_ij|.
    """
    a = np.array(a, dtype=float)
    b = np.array(b, dtype=float)
    n = len(b)
    scale = np.max(np.abs(a))
    for k in range(n - 1):
        if pivot:
            p = k + int(np.argmax(np.abs(a[k:, k])))
            a[[k, p]] = a[[p, k]]
            b[[k, p]] = b[[p, k]]
        for i in range(k + 1, n):
            m = a[i, k] / a[k, k]
            a[i, k:] -= m * a[k, k:]
            b[i] -= m * b[k]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - a[i, i + 1:] @ x[i + 1:]) / a[i, i]
    growth = np.max(np.abs(np.triu(a))) / scale
    return x, growth

def wilkinson(n):
    """1 on the diagonal and last column, -1 below the diagonal."""
    a = np.eye(n) - np.tril(np.ones((n, n)), -1)
    a[:, -1] = 1.0
    return a

def lauchli(eps):
    """Lauchli's 3 x 2 matrix [[1, 1], [eps, 0], [0, eps]]."""
    return np.array([[1.0, 1.0], [eps, 0.0], [0.0, eps]])

def rel_err(x, exact):
    """Relative error in the max norm."""
    return np.max(np.abs(x - exact)) / np.max(np.abs(exact))

def main():
    print("1. Hilbert systems H x = b with exact x = (1, ..., 1)")
    print("    n     kappa(H)   rel. error   kappa * u")
    for n in [4, 6, 8, 10, 12]:
        h = hilbert(n)
        x = np.linalg.solve(h, h @ np.ones(n))
        k = np.linalg.cond(h)
        print(f"   {n:2d}   {k:10.2e}   {rel_err(x, np.ones(n)):10.2e}"
              f"   {k * UNIT_ROUNDOFF:9.2e}")
    print()
    a, b = [[1e-20, 1.0], [1.0, 1.0]], [1.0, 2.0]
    x0, _ = lu_solve(a, b, pivot=False)
    x1, _ = lu_solve(a, b, pivot=True)
    print("2. [[1e-20, 1], [1, 1]] x = (1, 2); exact x = (1, 1) to 20 digits")
    print(f"   no pivoting:      x = ({x0[0]:.6f}, {x0[1]:.6f})")
    print(f"   partial pivoting: x = ({x1[0]:.6f}, {x1[1]:.6f})")
    print()
    print("3. Growth factor of partial pivoting on Wilkinson's matrix")
    for n in [10, 30, 55]:
        w = wilkinson(n)
        xs = np.linspace(1.0, 2.0, n)
        x, g = lu_solve(w, w @ xs)
        print(f"   n = {n:2d}: growth {g:9.3e} = 2^{np.log2(g):.0f},"
              f" rel. error {rel_err(x, xs):.1e}")
    print()
    eps = 1e-8
    m = lauchli(eps)
    rhs = m @ np.array([1.0, 1.0])
    gram = m.T @ m
    q, r = np.linalg.qr(m)
    xq = np.linalg.solve(r, q.T @ rhs)
    print("4. Lauchli least squares, eps = 1e-8, exact x = (1, 1)")
    print(f"   computed A^T A = {gram.tolist()}: det = {np.linalg.det(gram):.1e}")
    print(f"   QR solution: ({xq[0]:.8f}, {xq[1]:.8f})")
    print()
    print("Epistemic status: numerical experiments in IEEE double precision.")
    print("They illustrate the perturbation bound, the need for pivoting, the")
    print("worst-case growth factor and the squaring of the condition number")
    print("by the normal equations; the theorems themselves are in the lesson.")

if __name__ == "__main__":
    main()

The second program measures the order of three one-step methods by halving $h$ — the error ratios approach $2$, $4$ and $16$ — then shows stiffness: explicit Euler on $y' = -50(y - \cos t)$ explodes for $50h = 2.5$, is neutral at $50h = 2$, and is accurate for $50h \le 1.25$, while backward Euler is accurate for every step. Finally it shows the U-shaped total error of a difference quotient as $h$ shrinks past $\sqrt u$.

"""Order, stiffness and rounding: one-step methods for ODEs.

(1) Observed order: solve y' = y, y(0) = 1 on [0, 1] with Euler, Heun
(a second-order Runge-Kutta method) and classical RK4, halving the step
each time. If the global error is C h^p, halving h divides it by 2^p, so
log2 of successive error ratios estimates the order p: we see 1, 2 and 4.
(2) Stiffness: y' = -50 (y - cos t), y(0) = 0 has a smooth solution, but
explicit Euler is stable only for h * 50 < 2. With h = 0.05 it explodes;
backward Euler is stable for every h. (3) Rounding: the difference
quotient (e^(1+h) - e)/h approximates e with truncation error ~ h e/2 and
rounding error ~ u e/h, so the total error is smallest near h = sqrt(u)
and grows again for smaller h.

Evidence: numerical, compared with exact solutions.
"""

import math

def euler_step(f, t, y, h):
    """Explicit Euler."""
    return y + h * f(t, y)

def heun_step(f, t, y, h):
    """Heun's method (explicit trapezoid rule), order 2."""
    k1 = f(t, y)
    k2 = f(t + h, y + h * k1)
    return y + h * (k1 + k2) / 2

def rk4_step(f, t, y, h):
    """Classical fourth-order Runge-Kutta method."""
    k1 = f(t, y)
    k2 = f(t + h / 2, y + h * k1 / 2)
    k3 = f(t + h / 2, y + h * k2 / 2)
    k4 = f(t + h, y + h * k3)
    return y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6

def integrate(step, f, y0, t_end, n):
    """Take n equal steps of a one-step method from t = 0."""
    h = t_end / n
    t, y = 0.0, y0
    for _ in range(n):
        y = step(f, t, y, h)
        t += h
    return y

def backward_euler_stiff(lam, y0, t_end, n):
    """Backward Euler for y' = lam (y - cos t): solve the linear step exactly."""
    h = t_end / n
    t, y = 0.0, y0
    for _ in range(n):
        t += h
        y = (y - h * lam * math.cos(t)) / (1 - h * lam)
    return y

def stiff_exact(t):
    """Exact solution of y' = -50 (y - cos t), y(0) = 0."""
    a, b = 2500 / 2501, 50 / 2501
    return a * math.cos(t) + b * math.sin(t) - a * math.exp(-50 * t)

def main():
    f = lambda t, y: y
    print("1. y' = y on [0, 1]: error at t = 1 and observed order")
    print("     n      Euler   order       Heun   order        RK4   order")
    prev = None
    for n in [10, 20, 40, 80]:
        errs = [abs(integrate(s, f, 1.0, 1.0, n) - math.e)
                for s in (euler_step, heun_step, rk4_step)]
        cells = []
        for k, e in enumerate(errs):
            order = "" if prev is None else f"{math.log2(prev[k] / e):.2f}"
            cells.append(f"{e:9.2e}  {order:>5}")
        print((f"  {n:4d}   " + "   ".join(cells)).rstrip())
        prev = errs
    print()
    g = lambda t, y: -50 * (y - math.cos(t))
    print("2. Stiff: y' = -50 (y - cos t), y(0) = 0; error at t = 1")
    for n in [20, 25, 40, 100]:
        e1 = abs(integrate(euler_step, g, 0.0, 1.0, n) - stiff_exact(1.0))
        e2 = abs(backward_euler_stiff(-50, 0.0, 1.0, n) - stiff_exact(1.0))
        print(f"   h = {1 / n:5.3f} (50h = {50 / n:4.2f}): explicit {e1:9.2e},"
              f"  backward {e2:8.2e}")
    print()
    print("3. Forward difference for (e^x)' at x = 1: total error")
    for k in [2, 4, 6, 8, 10, 12, 14]:
        h = 10.0 ** -k
        err = abs((math.exp(1 + h) - math.e) / h - math.e)
        print(f"   h = 1e-{k:02d}: error {err:.1e}")
    print()
    print("Epistemic status: numerical experiments. The orders 1, 2, 4, the")
    print("stability limit 50h < 2 and the optimal h near sqrt(u) are")
    print("theorems or estimates proved in the lesson; these runs illustrate")
    print("them on single examples and prove nothing in general.")

if __name__ == "__main__":
    main()

These are numerical experiments on single examples. They illustrate the theorems, and the agreement with the predicted rates is good evidence that the code is right, but they prove nothing about other matrices or equations.

Ill-conditioning. This matrix has determinant $0.05$: it squashes the unit square into a sliver, and its condition number is about $80$. Move the $(2, 2)$ entry towards $1$ and watch the parallelogram collapse onto a line — the matrix becomes singular, and $\kappa \to \infty$.