Formula & theorem reference

Continuity on the whole closed interval is what does the work — it rules out the graph ever jumping over a height without passing through it. The polynomial-only version from College Algebra (the sign-change test, N = 0) is the special case where f(a) and f(b) have opposite signs; this is that same argument with the hypothesis weakened from "polynomial" to "continuous," and the conclusion widened from "crosses zero" to "hits any value in between." A rigorous construction (typically the least-upper-bound / completeness argument, showing c is the supremum of \{x \in [a,b] : f(x) < N\}) belongs to real analysis; here it's enough to know continuity is the ONLY property the proof actually needs.

Existence of p is the Intermediate Value Theorem, applied once to the starting bracket. For the bound: each bisection step replaces the current bracket with a new one exactly half as wide, and the replacement is always the half that still contains a sign change — so it still contains p (since the endpoint that gets discarded shares c_n's sign, the discarded half cannot be the one holding a root by the same IVT argument that guaranteed one to begin with). By induction, the bracket containing p after n steps has width (b-a)/2^n. The midpoint c_n of that bracket is, at worst, half the bracket's own width away from either endpoint, and p lies inside the bracket — so |c_n - p| is at most half the bracket's width, which is already (b-a)/2^n for the bracket produced BEFORE computing c_n; tracking the indices carefully against the standard convention (where c_n is the midpoint used to produce the n-th bracket) gives exactly |c_n-p|\le (b-a)/2^n. Since (b-a)/2^n \to 0, the squeeze forces c_n \to p.

Existence. h(x) = g(x) - x is continuous on [a,b]. Since g maps [a,b] into itself, g(a) \ge a and g(b) \le b. So h(a) \ge 0 and h(b) \le 0. The Intermediate Value Theorem gives some p \in [a,b] with h(p) = 0, i.e. g(p) = p.

Uniqueness. Suppose p and q were both fixed points in [a,b], with p \ne q. The Mean Value Theorem gives some \xi between them with g(p) - g(q) = g'(\xi)(p - q) \quad\Longrightarrow\quad p - q = g'(\xi)(p-q) (using g(p)=p, g(q)=q). Since p \ne q, dividing both sides by p-q gives 1 = g'(\xi). That contradicts |g'(\xi)| \le k < 1. So p is the only fixed point.

Convergence. Because g maps [a,b] into itself, every iterate x_n stays in [a,b]. So the Mean Value Theorem applies at every step. For some \xi_n between x_n and p, x_{n+1} - p = g(x_n) - g(p) = g'(\xi_n)(x_n - p) \quad\Longrightarrow\quad |x_{n+1}-p| \le k\,|x_n - p|. Apply this bound n times in a row. That gives |x_n - p| \le k^n |x_0 - p| \le k^n \max(x_0-a,\,b-x_0). Since 0 \le k < 1, k^n \to 0. That forces x_n \to p.

Since f'(p) \ne 0 and f' is continuous, f' stays bounded away from zero on some interval around p. That interval is where \delta lives. Shrink it further if needed, so that f'' (also continuous, by hypothesis) is bounded there too. Expand f in a first-order Taylor polynomial about x_n, evaluated at p, using the Lagrange remainder (taylors-theorem-lagrange-remainder, ch.00): f(p) = f(x_n) + f'(x_n)(p - x_n) + \frac{f''(\xi_n)}{2}(p-x_n)^2 for some \xi_n between x_n and p. The left side is 0, since f(p) = 0. Divide through by f'(x_n). This is legal, since f'(x_n) \ne 0 once x_n is in the interval where f' stays bounded away from zero: 0 = \frac{f(x_n)}{f'(x_n)} + (p - x_n) + \frac{f''(\xi_n)}{2f'(x_n)}(p-x_n)^2. Rearrange, using x_{n+1} = x_n - f(x_n)/f'(x_n) to replace the first two terms: x_{n+1} - p = -\frac{f''(\xi_n)}{2f'(x_n)}(x_n - p)^2. This is the claimed error identity. f'' and 1/f' are both bounded near p. So the coefficient -f''(\xi_n)/(2f'(x_n)) stays bounded by some constant M as x_n \to p. This gives |x_{n+1}-p| \le M|x_n-p|^2. Once |x_n - p| is small enough that M|x_n-p| < 1, each step's error is a shrinking multiple of the SQUARE of the last one. That forces x_n \to p, and makes the convergence quadratic.

Existence, by direct construction. Write the unknown polynomial as P(x) = a_0 + a_1x + a_2x^2 + \cdots + a_nx^n. The n+1 interpolation conditions P(x_i) = y_i become a linear system in the unknown coefficients a_0, \ldots, a_n: \begin{bmatrix} 1 & x_0 & x_0^2 & \cdots & x_0^n \\ 1 & x_1 & x_1^2 & \cdots & x_1^n \\ \vdots & & & & \vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^n \end{bmatrix}\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{bmatrix}=\begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix}. The coefficient matrix is a Vandermonde matrix, whose determinant is the classical identity \prod_{0\le i<j\le n}(x_j - x_i) — a product of factors that are each nonzero exactly because the nodes were assumed distinct. A nonzero determinant means the matrix is invertible, so this system has a solution (a_0, \ldots, a_n) for any right-hand side y_0, \ldots, y_n whatsoever. That solution's coefficients define a polynomial of degree at most n meeting every interpolation condition, so a solution exists.

Uniqueness. Suppose P and Q are both polynomials of degree at most n with P(x_i) = Q(x_i) = y_i for every i = 0, \ldots, n. The difference D = P - Q is then a polynomial of degree at most n with D(x_i) = 0 at all n+1 of the distinct points x_0, \ldots, x_n. A nonzero polynomial of degree at most n can have at most n roots — so a polynomial of that degree bound with n+1 distinct roots must be the zero polynomial. Hence D \equiv 0, i.e. P \equiv Q: the two polynomials are identical, not merely equal at the nodes.

Each L_i has degree exactly n (a product of n linear factors), so P, as a sum of degree-n polynomials, has degree at most n — the right degree bound for n+1 nodes.

The whole formula rests on one property of each L_i: it equals 1 at its own node and 0 at every other node. At x=x_i, every factor in L_i's product is of the form \frac{x_i-x_j}{x_i-x_j}=1, so L_i(x_i)=1. At x=x_k for some k\ne i, the product defining L_i includes the factor \frac{x-x_k}{x_i-x_k} (since k \ne i, this factor is genuinely present), which is 0 when x=x_k — so the whole product L_i(x_k)=0.

Evaluating P at a node x_k, every term y_iL_i(x_k) with i\ne k vanishes (since L_i(x_k)=0), leaving only y_kL_k(x_k)=y_k\cdot1=y_k. So P(x_k)=y_k for every k=0,\ldots,n: P meets every interpolation condition by construction, which is existence — the interpolation existence-uniqueness theorem then guarantees this is the interpolating polynomial, not merely one candidate.

Let P_{0..k} denote the (unique) degree-at-most-k polynomial interpolating f at x_0,\ldots,x_k, and, for the induction below, P_{1..k} the degree-at-most-(k-1) polynomial interpolating f at x_1,\ldots,x_k. It suffices to prove, by induction on k, that f[x_0,\ldots,x_k] is the coefficient of x^k in P_{0..k} — the Newton form's own structure (each new term adds exactly one new node and one new factor) then builds P_{0..k} from P_{0..k-1} automatically, since P_{0..k}-P_{0..k-1} has degree at most k and vanishes at x_0,\ldots,x_{k-1}, forcing it to equal (leading coefficient of P_{0..k})\cdot\prod_{i=0}^{k-1}(x-x_i).

Base case, k=0: P_{0..0}(x)=f(x_0) is the constant polynomial, and its "coefficient of x^0" is f(x_0)=f[x_0].

Inductive step. Assume f[x_0,\ldots,x_{k-1}] is P_{0..k-1}'s leading coefficient and f[x_1,\ldots,x_k] is P_{1..k}'s leading coefficient. Define Q(x) = \frac{(x-x_0)\,P_{1..k}(x) - (x-x_k)\,P_{0..k-1}(x)}{x_k-x_0}. Check Q interpolates f at every one of x_0,\ldots,x_k: at an interior node x_j (0<j<k), both P_{1..k}(x_j)=f(x_j) and P_{0..k-1}(x_j)=f(x_j), so Q(x_j)=\frac{(x_j-x_0)f(x_j)-(x_j-x_k)f(x_j)}{x_k-x_0}=\frac{(x_k-x_0)f(x_j)}{x_k-x_0}=f(x_j). At x_0, the second term's factor (x-x_k)P_{0..k-1}(x) is unaffected by the first term vanishing... directly: Q(x_0)=\frac{0-(x_0-x_k)P_{0..k-1}(x_0)}{x_k-x_0}=\frac{(x_k-x_0)f(x_0)}{x_k-x_0}=f(x_0) since P_{0..k-1}(x_0)=f(x_0). Symmetrically Q(x_k)=f(x_k) using P_{1..k}(x_k)=f(x_k). So Q is a degree-at-most-k polynomial matching f at all k+1 nodes x_0,\ldots,x_k — by the interpolation existence-uniqueness theorem, Q=P_{0..k}.

Q's coefficient of x^k comes only from the x^{k-1} leading terms of P_{1..k} and P_{0..k-1} (each multiplied by the linear factor x-x_0 or x-x_k), giving [\text{leading coeff of } P_{0..k}] = \frac{f[x_1,\ldots,x_k]-f[x_0,\ldots,x_{k-1}]}{x_k-x_0} = f[x_0,\ldots,x_k], exactly the divided-difference recursion applied to the full node range — completing the induction.

This is the identity proved in the previous section (Newton's divided-difference formula) in general form, reindexed. That proof showed, for the specific block of nodes x_0,\ldots,x_k, that Q(x) = \frac{(x-x_0)\,P_{1..k}(x) - (x-x_k)\,P_{0..k-1}(x)}{x_k-x_0} interpolates f at every node x_0,\ldots,x_k. Here P_{1..k} and P_{0..k-1} are the two length-k sub-interpolants obtained by dropping one endpoint each. The identical argument applies verbatim to any contiguous block of nodes x_{i-j},\ldots,x_i in place of x_0,\ldots,x_k. It works by checking agreement at the shared interior nodes, then at each of the two "new" endpoints in turn, using nothing but P_{1..k}'s and P_{0..k-1}'s own interpolation conditions. Here Q_{i,j-1} plays the role of P_{1..k} (dropping the leftmost node x_{i-j}), and Q_{i-1,j-1} plays the role of P_{0..k-1} (dropping the rightmost node x_i). By the interpolation existence-uniqueness theorem, Q_{i,j} is therefore the unique degree-at-most-j interpolant of that node block, evaluated at x. Taking i=n, j=n covers all n+1 nodes. So Q_{n,n}=P(x).

Fix x \in [a,b]. If x equals one of the nodes x_i, both sides are 0: the product on the right vanishes, and f(x_i) - P_n(x_i) = 0 by the interpolation conditions. So the equality holds trivially; pick any \xi(x) in the interval. Assume from here that x is not a node.

Let W(t) = \prod_{i=0}^{n}(t - x_i), a monic polynomial of degree n+1. Define g(t) = f(t) - P_n(t) - \frac{f(x) - P_n(x)}{W(x)}\,W(t). W(x) \ne 0 since x is not a node. So g is well-defined. By construction, g(x_i) = f(x_i) - P_n(x_i) - (\cdots)W(x_i) = 0 - 0 = 0 for each of the n+1 nodes. Also g(x) = f(x) - P_n(x) - \frac{f(x)-P_n(x)}{W(x)}W(x) = 0. So g has at least n+2 distinct roots in [a,b]: the n+1 nodes, together with x itself.

Apply Rolle's theorem repeatedly. Between each pair of consecutive roots of g, g' has a root. That gives g' at least n+1 roots. Between each pair of consecutive roots of g', g'' has a root. That gives g'' at least n roots. Continuing this way, g^{(n+1)} has at least (n+2)-(n+1) = 1 root in [a,b]. Call it \xi = \xi(x). It lies strictly between the smallest and largest of the n+2 points that produced it, exactly as claimed.

Differentiate g itself n+1 times. P_n has degree at most n, so P_n^{(n+1)} \equiv 0. W is monic of degree n+1, so its (n+1)-th derivative is the constant (n+1)! (the same fact used for any monic degree-(n+1) polynomial: differentiating t^{n+1} that many times leaves (n+1)!, and every lower-degree term of W vanishes after n+1 derivatives). So g^{(n+1)}(t) = f^{(n+1)}(t) - \frac{f(x)-P_n(x)}{W(x)}(n+1)!. Evaluate at t = \xi, where g^{(n+1)}(\xi) = 0. Solving for f(x) - P_n(x) gives f(x) - P_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}W(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^{n}(x-x_i), which is the claimed formula.

The formula comes from a geometric construction, not an algebraic derivation: place n+1 points equally spaced in angle around the upper half of the unit circle, at angles \theta_k=\dfrac{(2k+1)\pi}{2(n+1)} for k=0,\ldots,n — these angles are themselves equally spaced between 0 and \pi, each one centered in its own equal slice of the semicircle. Project each point straight down onto the horizontal diameter: the point at angle \theta_k projects to x_k=\cos\theta_k, giving exactly the formula above.

This projection is what produces the clustering: cosine changes slowly near \theta=0 and \theta=\pi (near the semicircle's two ends) and quickly near \theta=\pi/2 (the top of the arc) — so equally spaced angles project to x-values that are close together near \pm1 and spread apart near 0, the opposite of equally spaced x-values. Nodes distinct: since \theta\mapsto\cos\theta is strictly decreasing on [0,\pi] and the n+1 angles are themselves distinct, the n+1 nodes x_k are guaranteed distinct too.

The deeper reason this specific spacing controls the Runge phenomenon — that it minimizes the largest value \max_{x\in[-1,1]}\left|\prod_{i=0}^n(x-x_i)\right| over every possible choice of n+1 nodes — is a genuine theorem (Chebyshev's equioscillation property of the monic polynomials 2^{-n}T_{n+1}(x)), whose proof belongs to approximation theory and is outside this course's scope. What this course verifies directly instead, numerically, is that the minimized node-product term actually keeps the interpolation error controlled in practice — see the next section's own worked example.

The one new idea is the substitution itself. An ordinary divided difference f[x_i,x_i']=\dfrac{f(x_i')-f(x_i)}{x_i'-x_i} is exactly a difference quotient. As x_i'\to x_i, that quotient has a limit. By the definition of the derivative, that limit is f'(x_i). So setting f[z_{2i},z_{2i+1}]:=f'(x_i) is not an arbitrary rule. It is the value the ordinary divided-difference formula would converge to. That holds if the two repeated nodes were allowed to merge continuously, instead of being equal outright. So the divided-difference table stays a meaningful, well-defined object. This holds even where the standard formula's denominator would otherwise vanish.

Make that one substitution at each repeated pair. Then every remaining entry in the table is built by the same recursion Section 3 already proved correct, applied to the 2n+2-node list z_0,\ldots,z_{2n+1} (treating repeated values as distinct list POSITIONS, not distinct numbers). So the resulting Newton form is a well-defined degree-at-most-(2n+1) polynomial, by that same argument. It also meets the derivative conditions H'(x_i)=f'(x_i), not just the value conditions H(x_i)=f(x_i). That follows from differentiating the Newton form term by term, and confirming each derivative condition against the construction. It is verified directly, on real numbers, in this section's own worked example, rather than carried out here in full symbolic generality.

Fix x \in [a,b]. If x equals one of the nodes x_i, both sides are 0: the product on the right vanishes, and f(x_i) - P_n(x_i) = 0 by the interpolation conditions. So the equality holds trivially; pick any \xi(x) in the interval. Assume from here that x is not a node.

Let W(t) = \prod_{i=0}^{n}(t - x_i), a monic polynomial of degree n+1. Define g(t) = f(t) - P_n(t) - \frac{f(x) - P_n(x)}{W(x)}\,W(t). W(x) \ne 0 since x is not a node. So g is well-defined. By construction, g(x_i) = f(x_i) - P_n(x_i) - (\cdots)W(x_i) = 0 - 0 = 0 for each of the n+1 nodes. Also g(x) = f(x) - P_n(x) - \frac{f(x)-P_n(x)}{W(x)}W(x) = 0. So g has at least n+2 distinct roots in [a,b]: the n+1 nodes, together with x itself.

Apply Rolle's theorem repeatedly. Between each pair of consecutive roots of g, g' has a root. That gives g' at least n+1 roots. Between each pair of consecutive roots of g', g'' has a root. That gives g'' at least n roots. Continuing this way, g^{(n+1)} has at least (n+2)-(n+1) = 1 root in [a,b]. Call it \xi = \xi(x). It lies strictly between the smallest and largest of the n+2 points that produced it, exactly as claimed.

Differentiate g itself n+1 times. P_n has degree at most n, so P_n^{(n+1)} \equiv 0. W is monic of degree n+1, so its (n+1)-th derivative is the constant (n+1)! (the same fact used for any monic degree-(n+1) polynomial: differentiating t^{n+1} that many times leaves (n+1)!, and every lower-degree term of W vanishes after n+1 derivatives). So g^{(n+1)}(t) = f^{(n+1)}(t) - \frac{f(x)-P_n(x)}{W(x)}(n+1)!. Evaluate at t = \xi, where g^{(n+1)}(\xi) = 0. Solving for f(x) - P_n(x) gives f(x) - P_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}W(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^{n}(x-x_i), which is the claimed formula.

By construction. Writing each piece as S_j(x)=a_j+b_j(x-x_j)+c_j(x-x_j)^2+d_j(x-x_j)^3 on [x_j,x_{j+1}], the interpolation conditions fix every a_j=y_j, and matching S' and S'' at each interior node reduces b_j and d_j to explicit formulas in c_j and c_{j+1} alone (exactly the algebra this chapter's own construction carries out). What remains is showing the resulting tridiagonal system for c_0,\ldots,c_n — with c_0=c_n=0 fixed by the natural condition, and h_{j-1}c_{j-1}+2(h_{j-1}+h_j)c_j+h_jc_{j+1} =\alpha_j for each interior node j=1,\ldots,n-1 — has one and only one solution. That's a claim about n-1 equations in n-1 unknowns (c_1,\ldots,c_{n-1}), so it suffices to show the forward-elimination step this chapter's own algorithm performs never divides by zero, since a square linear system reduced by elimination to a triangular system with every pivot nonzero has exactly one solution (Chapter 0's own review of Gaussian elimination).

The proof is by induction on the elimination's own pivots l_1,\ldots, l_{n-1}, defined by l_1=2(h_0+h_1) and, for i\ge2, l_i=2(h_{i-1}+h_i)-h_{i-1}\mu_{i-1} where \mu_{i-1}=h_{i-1}/l_{i-1} — exactly the running quantities the construction computes along the way. Claim: l_i>h_i>0 for every i=1,\ldots,n-1.

Base case, i=1: l_1=2(h_0+h_1)>h_1 since 2h_0+h_1>0, and every h_j=x_{j+1}-x_j>0 because the nodes are distinct and increasing.

Inductive step. Suppose l_{i-1}>h_{i-1}>0 for some i\ge2. Then 0<\mu_{i-1}=h_{i-1}/l_{i-1}<1, so h_{i-1}\mu_{i-1}<h_{i-1}, and l_i = 2(h_{i-1}+h_i) - h_{i-1}\mu_{i-1} > 2(h_{i-1}+h_i)-h_{i-1} = h_{i-1}+2h_i > h_i > 0. This closes the induction: every pivot the elimination produces is strictly positive, hence never zero, so the elimination completes and the tridiagonal system has exactly one solution for c_1,\ldots,c_{n-1}. Together with c_0=c_n=0, every c_j is determined; b_j and d_j then follow from the explicit formulas already derived from c_j and c_{j+1}, and a_j=y_j was fixed from the start. Every coefficient of every piece is therefore uniquely determined, which is exactly the claim.

Apply Taylor's theorem with the Lagrange remainder (ch.00, Section 1) to f at x, expanded to the cubic term, once in each direction: f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{6}f'''(\xi_1), \qquad \xi_1 \in (x, x+h), f(x-h) = f(x) - hf'(x) + \frac{h^2}{2}f''(x) - \frac{h^3}{6}f'''(\xi_2), \qquad \xi_2 \in (x-h, x). Subtracting the second expansion from the first, the f(x) and \frac{h^2}{2}f''(x) terms cancel exactly (they carry the same sign in both expansions), leaving f(x+h) - f(x-h) = 2hf'(x) + \frac{h^3}{6}\big(f'''(\xi_1) + f'''(\xi_2)\big). Dividing by 2h and solving for f'(x), f'(x) = \frac{f(x+h)-f(x-h)}{2h} - \frac{h^2}{12}\big(f'''(\xi_1) + f'''(\xi_2)\big). Since f''' is continuous on [x-h, x+h] and \xi_1, \xi_2 both lie in that interval, the average \frac{1}{2}\big(f'''(\xi_1)+f'''(\xi_2)\big) lies between the minimum and maximum of f''' on [x-h,x+h]. By the Intermediate Value Theorem, f''' attains that average value at some point \xi in [x-h, x+h] (strictly between the endpoints, since it is a genuine average of two values attained at two interior-or-boundary points \xi_1 \ne \xi_2 unless f''' is constant, in which case any \xi in the interval works). Substituting \frac{1}{2}(f'''(\xi_1)+f'''(\xi_2)) = f'''(\xi) gives f'(x) = \frac{f(x+h)-f(x-h)}{2h} - \frac{h^2}{6}f'''(\xi), which is the claimed formula.

The n=2 case is small enough to derive directly by the method of undetermined coefficients, and the pattern generalizes. Seek x_1, x_2, w_1, w_2 so that w_1f(x_1)+w_2f(x_2) exactly reproduces \int_{-1}^1 f(x)\,dx for every polynomial up to degree 3. That is four free parameters, matched against four conditions, one per basis polynomial 1,x,x^2,x^3: \begin{aligned} f(x)=1:&\quad w_1+w_2 = \int_{-1}^1 1\,dx = 2\\ f(x)=x:&\quad w_1x_1+w_2x_2 = \int_{-1}^1 x\,dx = 0\\ f(x)=x^2:&\quad w_1x_1^2+w_2x_2^2 = \int_{-1}^1 x^2\,dx = \tfrac{2}{3}\\ f(x)=x^3:&\quad w_1x_1^3+w_2x_2^3 = \int_{-1}^1 x^3\,dx = 0 \end{aligned} Symmetry (nothing distinguishes +x from -x on [-1,1]) suggests x_2=-x_1 and w_1=w_2. That satisfies the second and fourth equations immediately. The first equation then gives w_1=w_2=1. The third gives 2x_1^2=\tfrac23, so x_1=\tfrac{1}{\sqrt3}, x_2=-\tfrac{1}{\sqrt3}. This is exactly the n=2 row of the standard table. Four conditions are satisfied by only four unknowns, with no degree-4 condition imposed. So exactness up to degree 3=2(2)-1 is automatic, and exactness at degree 4 is not guaranteed (and generally fails). This is the same 2n-1 pattern that appears for every n.

For general n, the same idea applies: matching 2n moment conditions with 2n unknowns. This has a solution if and only if the nodes are chosen to be orthogonal to every polynomial of degree <n, under the standard inner product \langle p,q\rangle=\int_{-1}^1 p(x)q(x)\,dx. That is exactly what characterizes the roots of P_n(x), the degree-n member of the Legendre polynomial family (the orthogonal polynomials for this exact inner product). The full argument is a standard result in approximation theory and is not re-derived here. The n=2 case above already shows where the 2n-1 degree-of-precision figure comes from.

Separate the defining sum by parity of n: write n=2r for the even terms and n=2r+1 for the odd terms, each with r = 0,\dots,\tfrac N2 - 1: X_k = \sum_{r=0}^{N/2-1} x_{2r}\,e^{-i2\pi k(2r)/N} + \sum_{r=0}^{N/2-1} x_{2r+1}\,e^{-i2\pi k(2r+1)/N}. In the first sum, e^{-i2\pi k(2r)/N} = e^{-i2\pi kr/(N/2)} — exactly the kernel of a length-(N/2) DFT in r, evaluated at frequency index k. In the second sum, factor out the part not depending on r: e^{-i2\pi k(2r+1)/N} = e^{-i2\pi k/N}\cdot e^{-i2\pi kr/(N/2)}. So X_k = \underbrace{\sum_{r=0}^{N/2-1} x_{2r}\,e^{-i2\pi kr/(N/2)}}_{E_k} \;+\; e^{-i2\pi k/N}\underbrace{\sum_{r=0}^{N/2-1} x_{2r+1}\,e^{-i2\pi kr/(N/2)}}_{O_k} = E_k + W_N^{\,k}\,O_k. As length-(N/2) DFTs, E_k and O_k are periodic in k with period N/2: E_{k+N/2} = E_k and O_{k+N/2} = O_k. The twiddle factor, however, flips sign every half-period: W_N^{\,k+N/2} = e^{-i2\pi(k+N/2)/N} = e^{-i2\pi k/N}\cdot e^{-i\pi} = -W_N^{\,k}. So evaluating the same split sum at index k + N/2 (for k = 0,\dots,\tfrac N2-1) gives X_{k+N/2} = E_{k+N/2} + W_N^{\,k+N/2}O_{k+N/2} = E_k - W_N^{\,k}O_k — the two outputs X_k and X_{k+N/2} share the same pair (E_k, O_k) and differ only in the sign attached to W_N^{\,k}O_k, which is exactly one "butterfly": two inputs (E_k, O_k) producing two outputs (X_k, X_{k+N/2}) with a single complex multiplication (W_N^{\,k}O_k) and two additions.

Complexity. Let T(N) be the operation count for a length-N DFT computed this way. Computing E_k and O_k costs 2T(N/2), and combining them into all N outputs via N/2 butterflies costs O(N) more, giving the recurrence T(N) = 2T(N/2) + O(N). Unrolling this recursion \log_2 N times (until the base case T(1) = O(1) is reached) accumulates O(N) work at each of \log_2 N levels, for a total of T(N) = O(N\log N) — asymptotically far below the O(N^2) of evaluating the defining sum directly for every k.

Track the elimination as a sequence of matrix multiplications rather than as row-by-row bookkeeping. Each row replacement step — "add m times row j to row i" — is itself left-multiplication by an elementary matrix E: the identity matrix with m dropped into the (i,j) slot. If reducing A to U takes k such steps with elementary matrices E_1, E_2, \dots, E_k applied in that order, then

Every E_t here is unit lower triangular (1's on the diagonal, a single nonzero entry below it, from a row replacement that only ever adds a multiple of a higher row to a lower one — no swaps means row i only ever absorbs rows above it). The set of unit lower triangular matrices is closed under multiplication and under inversion, so the product M = E_k \cdots E_2 E_1 is itself unit lower triangular, and so is its inverse L = M^{-1}. Solving MA = U for A gives

It remains to identify the entries of L. Each E_t^{-1} undoes its row operation by flipping the sign of its one off-diagonal entry: if E_t adds m times row j to row i, then E_t^{-1} subtracts it back out, so E_t^{-1} carries -m in slot (i,j). Write \ell_{ij} = -m for that entry, i.e. \ell_{ij} is the negative of the multiplier used to clear (i,j) — equivalently, the multiplier used to reconstruct it from U's reduction of A. Because no swaps occur, the elimination order never revisits a cleared entry, so the individual E_t^{-1} matrices combine without their off-diagonal entries interfering with one another, and

comes out unit lower triangular with exactly \ell_{ij} sitting in each below-diagonal slot (i,j) — precisely the multiplier record from the forward elimination. This is why L can be written down directly from the row reduction as it happens, with no matrix inversion required in practice.

Sanity check: if A is already upper triangular, elimination needs no steps at all, so L = I and U = A — consistent with A = IA.

Since A\mathbf{x}=\mathbf{b},

Taking norms and using the induced-norm inequality \lVert M\mathbf{v}\rVert\le\lVert M\rVert\,\lVert\mathbf{v}\rVert (true of every matrix norm induced by a vector norm, by the very definition of \lVert M\rVert as the largest such ratio):

The same inequality applied to \mathbf{b}=A\mathbf{x} gives \lVert\mathbf{b}\rVert\le\lVert A\rVert\,\lVert\mathbf{x}\rVert, so

Dividing the first bound by this lower bound on \lVert\mathbf{x}\rVert:

Sanity check: if \hat{\mathbf{x}}=\mathbf{x} exactly, then \mathbf{r}=\mathbf{0} and both sides are 0 — the bound is tight in the one case where there's nothing left to bound.

Direct, via a contraction argument in the infinity norm (Chapter 7, Section 6). Write the Jacobi iteration as \mathbf{x}^{(k+1)}= T_J\mathbf{x}^{(k)}+\mathbf{c}, where T_J=-D^{-1}(L+U) has entries (T_J)_{ij}=-a_{ij}/a_{ii} for j\ne i and (T_J)_{ii}=0. The exact solution \mathbf{x}^* satisfies the SAME fixed-point equation, \mathbf{x}^*=T_J\mathbf{x}^*+\mathbf{c}, since \mathbf{x}^* solving A\mathbf{x}=\mathbf{b} is exactly what makes every row's own rearrangement (Section 2) come out true. Subtracting the two equations gives \mathbf{x}^{(k+1)}-\mathbf{x}^*=T_J\bigl(\mathbf{x}^{(k)}- \mathbf{x}^*\bigr).

The infinity norm of T_J is its largest absolute row sum: \lVert T_J\rVert_\infty=\max_i\frac{\sum_{j\ne i}\lvert a_{ij}\rvert} {\lvert a_{ii}\rvert}. Strict diagonal dominance means every one of these ratios is strictly less than 1, so \lVert T_J\rVert_\infty=\delta for some \delta<1. Taking infinity norms of the subtracted equation and using that an induced matrix norm is submultiplicative (\lVert T_J\mathbf{v}\rVert_ \infty\le\lVert T_J\rVert_\infty\lVert\mathbf{v}\rVert_\infty, Chapter 7, Section 6): \bigl\lVert\mathbf{x}^{(k+1)}-\mathbf{x}^*\bigr\rVert_\infty\le\delta \bigl\lVert\mathbf{x}^{(k)}-\mathbf{x}^*\bigr\rVert_\infty. Applying this bound repeatedly, by induction on k, \bigl\lVert\mathbf{x}^{(k)}-\mathbf{x}^*\bigr\rVert_\infty\le\delta^k \bigl\lVert\mathbf{x}^{(0)}-\mathbf{x}^*\bigr\rVert_\infty. Since 0\le\delta<1, \delta^k\to0 as k\to\infty, regardless of how far \mathbf{x}^{(0)} started from \mathbf{x}^* -- so \mathbf{x}^{(k)}\to\mathbf{x}^* from ANY starting guess. This is exactly the "geometric-decay to zero" pattern Chapter 2's convergence-order vocabulary already named linear convergence, with rate at most \delta.

Gauss-Seidel's iteration matrix T_{GS}=-(D+L)^{-1}U satisfies the same bound, \lVert T_{GS}\rVert_\infty\le\delta, by an analogous but more delicate row-by-row argument: each row's bound on the already-updated entries (using \delta<1 inductively on earlier rows) combines with the same diagonal-dominance margin on the not-yet-updated entries to keep the overall ratio under 1 as well. The convergence conclusion for Gauss-Seidel then follows by the identical contraction argument above, applied to T_{GS} in place of T_J.

A least-squares solution \hat{\mathbf x} is, by definition, a minimizer of \lVert A\mathbf x - \mathbf b \rVert over all \mathbf x — equivalently, a vector for which A\hat{\mathbf x} is the point of \operatorname{Col}(A) closest to \mathbf b. The closest-point characterization of a projection says exactly that the residual \mathbf b - A\hat{\mathbf x} must be orthogonal to \operatorname{Col}(A): no vector in that column space can be nearer to \mathbf b than its orthogonal projection, so if the residual had any component lying in \operatorname{Col}(A), subtracting that component would move A\hat{\mathbf x} closer to \mathbf b still, contradicting minimality.

Orthogonality to every vector in \operatorname{Col}(A) means orthogonality to every column of A in particular, since those columns span \operatorname{Col}(A). Writing \mathbf a_1, \dots, \mathbf a_n for the columns of A, this is the system of dot-product conditions

Stacking these n scalar equations into one vector equation is precisely what A^{\mathsf T} does: the j-th row of A^{\mathsf T} is \mathbf a_j^{\mathsf T}, so left-multiplying by A^{\mathsf T} produces a vector whose j-th entry is \mathbf a_j \cdot (\mathbf b - A\hat{\mathbf x}). Since every one of those entries is 0,

Distributing A^{\mathsf T} across the difference and moving the A^{\mathsf T} A \hat{\mathbf x} term to the other side gives

Sanity check: if \mathbf b already lies in \operatorname{Col}(A), the system A\mathbf x = \mathbf b is consistent and any exact solution \mathbf x^* leaves residual \mathbf b - A\mathbf x^* = \mathbf 0, which is trivially orthogonal to every column of A — so \mathbf x^* satisfies the normal equations too, and least squares reduces to ordinary solving exactly when there's no error left to minimize.

Write \mathbf{x} for a candidate vector and split the target \mathbf{b} against the subspace \operatorname{Col}(A): every \mathbf{b} decomposes uniquely as \mathbf{b} = \hat{\mathbf{b}} + \mathbf{b}^{\perp}, where \hat{\mathbf{b}} = \operatorname{proj}_{\operatorname{Col}(A)}\mathbf{b} is the closest point to \mathbf{b} inside \operatorname{Col}(A) and \mathbf{b}^{\perp} is orthogonal to every vector in \operatorname{Col}(A).

For any \mathbf{x}, the vector A\mathbf{x} lies in \operatorname{Col}(A), so the residual splits as A\mathbf{x} - \mathbf{b} = (A\mathbf{x} - \hat{\mathbf{b}}) - \mathbf{b}^{\perp}, a sum of a vector inside \operatorname{Col}(A) and a vector orthogonal to it. The Pythagorean theorem for orthogonal vectors then gives \lVert A\mathbf{x} - \mathbf{b} \rVert^2 = \lVert A\mathbf{x} - \hat{\mathbf{b}} \rVert^2 + \lVert \mathbf{b}^{\perp} \rVert^2. The second term does not depend on \mathbf{x} at all, so this quantity is smallest exactly when the first term is 0 — that is, exactly when A\mathbf{x} = \hat{\mathbf{b}}. Every minimizer \hat{\mathbf{x}} therefore satisfies A\hat{\mathbf{x}} = \hat{\mathbf{b}}, which is equivalent to \mathbf{b} - A\hat{\mathbf{x}} = \mathbf{b} - \hat{\mathbf{b}} = \mathbf{b}^{\perp} \in \operatorname{Col}(A)^{\perp}.

A vector lies in \operatorname{Col}(A)^{\perp} exactly when it is orthogonal to every column of A, i.e. when A^{\mathsf T} sends it to \mathbf{0} — because each entry of A^{\mathsf T}\mathbf{v} is one such dot product. So the residual condition above rewrites as A^{\mathsf T}\left(\mathbf{b} - A\hat{\mathbf{x}}\right) = \mathbf{0} \quad\Longleftrightarrow\quad A^{\mathsf T} A \hat{\mathbf{x}} = A^{\mathsf T} \mathbf{b}, the normal equations. Every step above is reversible, so a vector \hat{\mathbf{x}} satisfies the normal equations if and only if it minimizes \lVert A\mathbf{x} - \mathbf{b} \rVert.

For uniqueness under linear independence: if the columns of A are linearly independent then A\mathbf{v} = \mathbf{0} forces \mathbf{v} = \mathbf{0}, so for any nonzero \mathbf{v}, \mathbf{v}^{\mathsf T} A^{\mathsf T} A \mathbf{v} = (A\mathbf{v})^{\mathsf T}(A\mathbf{v}) = \lVert A\mathbf{v} \rVert^2 > 0, which shows A^{\mathsf T} A \mathbf{v} \neq \mathbf{0} for every nonzero \mathbf{v} — so A^{\mathsf T} A has trivial null space and, being square, is invertible. Multiplying both sides of the normal equations by \left(A^{\mathsf T} A\right)^{-1} then isolates the unique solution \hat{\mathbf{x}} = \left(A^{\mathsf T} A\right)^{-1} A^{\mathsf T} \mathbf{b}.

Run the Gram–Schmidt process on the columns \mathbf{a}_1, \dots, \mathbf{a}_n, producing orthonormal vectors \mathbf{q}_1, \dots, \mathbf{q}_n one at a time. At each step, build an intermediate vector \mathbf{w}_k by stripping off \mathbf{a}_k's component along every \mathbf{q}_i already found, \mathbf{w}_k = \mathbf{a}_k - \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\,\mathbf{q}_i, \qquad \mathbf{q}_k = \frac{\mathbf{w}_k}{\lVert\mathbf{w}_k\rVert}. Linear independence of \mathbf{a}_1, \dots, \mathbf{a}_n guarantees \mathbf{w}_k \neq \mathbf{0} at every step (a zero \mathbf{w}_k would make \mathbf{a}_k a combination of \mathbf{a}_1, \dots, \mathbf{a}_{k-1}), so each division is legal and \lVert\mathbf{w}_k\rVert > 0.

Solve the defining relation for \mathbf{a}_k itself by moving the sum to the other side and substituting \mathbf{w}_k = \lVert\mathbf{w}_k\rVert\, \mathbf{q}_k: \mathbf{a}_k = \sum_{i=1}^{k-1} (\mathbf{a}_k \cdot \mathbf{q}_i)\, \mathbf{q}_i \;+\; \lVert\mathbf{w}_k\rVert\,\mathbf{q}_k. Every vector on the right is one of \mathbf{q}_1, \dots, \mathbf{q}_k, so this says each column of A is a linear combination of the \mathbf{q}_i's with an explicit, computable set of coefficients — exactly the data needed to assemble a matrix equation.

Let Q be the m \times n matrix with columns \mathbf{q}_1, \dots, \mathbf{q}_n, and define R entry by entry from the coefficients just found: r_{ik} = \mathbf{a}_k \cdot \mathbf{q}_i \ \ (i < k), \qquad r_{kk} = \lVert\mathbf{w}_k\rVert, \qquad r_{ik} = 0 \ \ (i > k). The zero pattern is exactly upper triangular, since \mathbf{a}_k was only ever built from \mathbf{q}_1, \dots, \mathbf{q}_k — it carries no component along \mathbf{q}_{k+1}, \dots, \mathbf{q}_n, which is what r_{ik} = 0 for i > k records. With this R, column k of the product QR is \sum_{i=1}^{n} r_{ik}\,\mathbf{q}_i = \sum_{i=1}^{k} r_{ik}\,\mathbf{q}_i, which is precisely the expression for \mathbf{a}_k derived above. Since this holds for every column, A = QR.

Finally, Q has orthonormal columns by construction of Gram–Schmidt, so Q^{\mathsf T}Q = I_n, and R is upper triangular with diagonal entries r_{kk} = \lVert\mathbf{w}_k\rVert > 0, as required.

Sanity check: take n = 1, so A is a single nonzero column \mathbf{a}_1. Then \mathbf{w}_1 = \mathbf{a}_1, \mathbf{q}_1 = \mathbf{a}_1 / \lVert\mathbf{a}_1\rVert, and R = [\lVert\mathbf{a}_1\rVert] — a single positive number. Indeed QR = \mathbf{q}_1 \cdot \lVert\mathbf{a}_1\rVert = \mathbf{a}_1, recovering A exactly, as the factorization must in the smallest nontrivial case.

Since the columns of A are linearly independent, A^{\mathsf T}A is invertible and the least-squares solution is the unique solution of the normal equations A^{\mathsf T}A\,\mathbf{x} = A^{\mathsf T}\mathbf{b}.

Substitute A = QR into both sides. On the left, A^{\mathsf T}A = (QR)^{\mathsf T}(QR) = R^{\mathsf T}Q^{\mathsf T}QR. The columns of Q are orthonormal, which is exactly the statement Q^{\mathsf T}Q = I_n, so this collapses to A^{\mathsf T}A = R^{\mathsf T}R. On the right, A^{\mathsf T}\mathbf{b} = (QR)^{\mathsf T}\mathbf{b} = R^{\mathsf T}Q^{\mathsf T}\mathbf{b}. The normal equations become R^{\mathsf T}R\,\mathbf{x} = R^{\mathsf T}Q^{\mathsf T}\mathbf{b}.

R is invertible — it is a factor of the full-column-rank matrix A in a factorization where Q has full column rank — so R^{\mathsf T} is also invertible, and left-multiplying both sides by \left(R^{\mathsf T}\right)^{-1} cancels it entirely: R\mathbf{x} = Q^{\mathsf T}\mathbf{b}. Since R is invertible this system has the unique solution \hat{\mathbf{x}} = R^{-1}Q^{\mathsf T}\mathbf{b}, which must therefore be the least-squares solution guaranteed by the normal equations.

Why this is the practical version of the normal equations. Forming A^{\mathsf T}A squares the conditioning of the problem — round-off error that was small in A gets amplified when A is multiplied against itself. The QR route never forms A^{\mathsf T}A at all: it computes Q^{\mathsf T}\mathbf{b} directly and then solves the triangular system R\hat{\mathbf x} = Q^{\mathsf T}\mathbf{b} by back substitution, which is both cheaper and numerically safer. The two formulas agree exactly — R^{\mathsf T}R above is nothing but a triangular factorization of A^{\mathsf T}A itself — they just reach the same \hat{\mathbf{x}} by different arithmetic paths.

Write \mathbf x_0 in the eigenbasis: \mathbf x_0=c_1\mathbf v_1+c_2\mathbf v_2+\cdots+c_n\mathbf v_n with c_1\ne0 by hypothesis. A\mathbf v_i=\lambda_i\mathbf v_i, so applying A repeatedly gives A^k\mathbf x_0=c_1\lambda_1^k\mathbf v_1+c_2\lambda_2^k\mathbf v_2+\cdots+c_n\lambda_n^k\mathbf v_n =\lambda_1^k\left(c_1\mathbf v_1+\sum_{i=2}^n c_i\left(\frac{\lambda_i}{\lambda_1}\right)^k\mathbf v_i\right). Every ratio \lvert\lambda_i/\lambda_1\rvert<1 for i\ge2. That is exactly what "strictly dominant" means. So (\lambda_i/\lambda_1)^k\to0 as k\to\infty. The term in parentheses converges to c_1\mathbf v_1. The leading scalar factor \lambda_1^k changes only the LENGTH of A^k\mathbf x_0, never its direction. Normalizing at every step removes that factor. (That means dividing by \lVert A^k\mathbf x_0\rVert at each step, exactly what the algorithm does.) This leaves the direction converging to \mathbf v_1 itself, up to an overall sign. (\lambda_1^k can be negative when \lambda_1<0 and k is odd.) The slowest-decaying error term is the one with the largest surviving ratio, \lvert\lambda_2/ \lambda_1\rvert. That ratio governs the overall rate.

The normalized iterate \mathbf x_k converges to a unit eigenvector for \lambda_1. The Rayleigh quotient is continuous in \mathbf x, and exactly equals \lambda_1 when evaluated AT an exact unit eigenvector for \lambda_1 (\mathbf v_1^{\mathsf T}A\mathbf v_1/\mathbf v_1^{\mathsf T}\mathbf v_1=\mathbf v_1^{\mathsf T}(\lambda_1\mathbf v_1)/1=\lambda_1). So continuity carries the vector's own convergence over to the Rayleigh quotient's convergence to \lambda_1, at the same \lvert\lambda_2/\lambda_1\rvert^k rate.

A^{\mathsf T}A is an n\times n symmetric matrix, so the real spectral theorem gives an orthonormal eigenbasis \{\mathbf v_1,\dots,\mathbf v_n\} of \mathbb R^n with A^{\mathsf T}A\mathbf v_i = \lambda_i\mathbf v_i. Every eigenvalue is nonnegative: \lambda_i \lVert\mathbf v_i\rVert^2 = \mathbf v_i^{\mathsf T}A^{\mathsf T}A\mathbf v_i = \lVert A\mathbf v_i\rVert^2 \ge 0. Order them so \lambda_1\ge \cdots\ge\lambda_n\ge0 and set \sigma_i = \sqrt{\lambda_i}.

Let r = \text{rank}(A) (equivalently, the number of nonzero \lambda_i, since A\mathbf v_i=\mathbf 0 exactly when \lambda_i=0). For i\le r, define \mathbf u_i = \frac{1}{\sigma_i}A\mathbf v_i. These are orthonormal: \mathbf u_i \cdot \mathbf u_j = \frac{1}{\sigma_i\sigma_j}(A\mathbf v_i)^{\mathsf T}(A\mathbf v_j) = \frac{1}{\sigma_i\sigma_j}\mathbf v_i^{\mathsf T}A^{\mathsf T}A\mathbf v_j = \frac{\lambda_j}{\sigma_i\sigma_j}\mathbf v_i\cdot\mathbf v_j, which is 1 when i=j (since \lambda_i=\sigma_i^2) and 0 otherwise (orthonormality of the \mathbf v_i). Extend \{\mathbf u_1,\dots,\mathbf u_r\} to a full orthonormal basis \{\mathbf u_1,\dots,\mathbf u_m\} of \mathbb R^m (possible since any orthonormal set extends to a full orthonormal basis).

Let V have columns \mathbf v_1,\dots,\mathbf v_n and U have columns \mathbf u_1,\dots,\mathbf u_m — both orthogonal by construction. By definition A\mathbf v_i = \sigma_i\mathbf u_i for i\le r, and A\mathbf v_i=\mathbf 0=\sigma_i\mathbf u_i (since \sigma_i=0) for i>r — so AV = U\Sigma for the matrix \Sigma with \sigma_i down its diagonal. Since V is orthogonal, V^{-1} = V^{\mathsf T}, giving A = U\Sigma V^{\mathsf T} exactly. \blacksquare

This result is proved by constructing a sequence of successive approximations y_0(x), y_1(x), y_2(x), \dots (Picard iterates), each defined from the last by y_{n+1}(x) = y_0 + \int_{x_0}^{x} f(t, y_n(t))\,dt, and showing this sequence converges to a genuine solution — the continuity of f alone guarantees a solution exists (this weaker half is Peano's theorem). The extra hypothesis, continuity of \partial f/\partial y, is exactly what makes the map y_n \mapsto y_{n+1} a contraction on a small enough interval: two different starting guesses get pulled closer together with each iteration, at a rate bounded by a Lipschitz constant derived from \partial f/\partial y. A genuine second solution through the same point would be a second fixed point of that contraction, which a contraction mapping can never have. Neither half of this argument produces a formula — both are pure existence arguments, which is exactly why this theorem answers "is there a solution, and only one" without ever constructing it explicitly.

A method built by matching a Taylor expansion through the h^k term has local truncation error O(h^{k+1}) — one power of h higher than what it matched, since the first term it fails to match is the source of the error. Reaching a fixed endpoint from a fixed starting point takes n = (\text{interval length})/h steps, so the global error is (local error per step) \times (number of steps) \sim h^{k+1}\cdot(1/h) = h^k — one power of h lower than the local error. This is why Euler's method (matches through h^1, local error O(h^2)) is global order 1, while Heun's method (matches through h^2, local error O(h^3)) is global order 2, and RK4 (matches through h^4, local error O(h^5)) is global order 4: reading a method's global order directly off its local error requires subtracting exactly one power of h.

k_1 is the slope at the left endpoint; k_2 and k_3 are two independent estimates of the slope at the midpoint, the second refining the first by using k_2's own predicted y-value instead of k_1's; k_4 is the slope at the right endpoint, predicted using the midpoint slope k_3. Weighting these 1:2:2:1 and dividing by 6 reproduces Simpson's rule's own weighting for approximating an integral by samples at both endpoints and (doubly) the midpoint. Carrying this weighted average through a Taylor expansion of the true solution shows it matches every term through h^4, giving a local truncation error of O(h^5) per step and a global error of O(h^4) — two full orders better than Euler's method, at the cost of four function evaluations per step instead of one.

Both formulas are constructed the way any Runge-Kutta method is. Each set of weights is chosen so the resulting y_{n+1} matches the true solution's Taylor expansion through the stated order. What makes RKF45 different: Fehlberg found one set of six stage evaluations k_1,\ldots,k_6 rich enough to support two different weighted combinations at once. One matches the Taylor expansion through h^4 (local truncation error O(h^5)). The other matches it through h^5 (local truncation error O(h^6)). Neither needs a second, independent set of stages for the higher-order estimate. Both estimates are built from the exact same six function evaluations, so their difference is available essentially for free. The 5th-order estimate's own error is one full power of h smaller than the 4th-order estimate's. So that difference is dominated by the 4th-order formula's own local truncation error. That gives an estimate of the error a step is making, obtained without ever knowing the true solution.

Both formulas come from replacing f, over the interval [t_n,t_{n+1}], by the polynomial that interpolates it at the known grid points. Then integrate \int_{t_n}^{t_{n+1}} f\,dt \approx y_{n+1}-y_n against that polynomial's own Lagrange form. The Adams-Bashforth predictor interpolates f at the four points already computed, t_{n-3},\ldots,t_n. Every node lies AT OR BEFORE t_n. So integrating over [t_n,t_{n+1}] is an extrapolation past the last known point. The Adams-Moulton corrector interpolates f at t_{n-2},\ldots,t_{n+1} instead, one step further along. So the same integral now falls INSIDE the node range [t_{n-2},t_{n+1}]: a genuine interpolation, not an extrapolation. This is exactly why the corrector is more accurate, and more stable, than a predictor built from the same number of points. Interpolation is intrinsically better-conditioned than extrapolation. The coefficients above are precisely what integrating each interpolating cubic against its own Lagrange basis produces.

Apply any explicit Runge-Kutta method to y'=\lambda y. Every stage k_i comes out as \lambda times a fixed linear combination of y_n and the earlier stages, since f(t,y)=\lambda y is itself linear. So y_{n+1} works out to R(h\lambda)y_n for some polynomial R built entirely from the method's own weights. For Euler, R(z)=1+z. For RK4, R(z)=1+z+\tfrac{z^2}{2}+\tfrac{z^3}{6}+\tfrac{z^4}{24}: the Taylor series of e^z truncated at the method's own order. That is exactly what a consistent method has to reproduce, since the test equation's true solution is y=y_0e^{\lambda t}. Iterate y_{n+1}=R(z)y_n for n steps. That gives y_n=R(z)^n y_0, which stays bounded as n grows exactly when |R(z)|\le1, and grows without bound otherwise. This is the numerical analogue of the true solution's own boundedness for \operatorname{Re}(\lambda)\le0. The exact solution respects that bound at every step size, but a numerical method only respects it for h\lambda lying inside its own region.