Apps

You know the values of a function \(f\) at \(n+1\) different points \(x_0\lt x_1\lt \dots\lt x_n\), the nodes, say \(y_j=f(x_j)\). The values may come from a measurement, or from a formula that is expensive to evaluate. Polynomial interpolation looks for a polynomial \(p\) of degree at most \(n\) that takes the same values at the nodes:

Why polynomials? They are the simplest functions there are: evaluating one takes only additions and multiplications, and it can be differentiated and integrated exactly. Weierstrass' theorem says that every continuous function on a closed interval is a uniform limit of polynomials. And interpolation is the engine inside many other methods: Simpson's rule and all quadrature formulas, finite differences, the Adams methods for differential equations, the secant method and inverse interpolation for roots, Richardson extrapolation and Romberg integration.

The questions of this course are the natural ones. Does such a \(p\) exist, and is it unique? How do we compute it, and which way is best? How large is \(f-p\) between the nodes? Can the nodes be chosen well? What happens when \(n\) grows?

Theorem. If the nodes \(x_0,\dots,x_n\) are distinct, then for any values \(y_0,\dots,y_n\) there is exactly one polynomial \(p\) of degree at most \(n\) with \(p(x_j)=y_j\) for all \(j\).

Uniqueness. If \(p\) and \(q\) both interpolate, then \(p-q\) has degree at most \(n\) and vanishes at the \(n+1\) distinct nodes. A nonzero polynomial of degree at most \(n\) has at most \(n\) roots, so \(p-q=0\).

Existence. Write \(p(x)=a_0+a_1x+\dots+a_nx^n\). The conditions form a linear system \(Va=y\) with the Vandermonde matrix

Uniqueness says that \(Va=0\) only for \(a=0\), so \(V\) is invertible and the system has a solution. (The next group builds the solution explicitly, which is a second proof.)

Consequences. If the data come from a polynomial of degree at most \(n\), the interpolant is that very polynomial: choose “polynomial x⁵ − 2x³ + x/2” in the app and raise \(n\) from 4 to 5; the error jumps from 0.11 to about \(10^{-15}\). The degree of \(p\) can be smaller than \(n\) (all \(y_j\) equal gives a constant). The nodes must be distinct; with repeated nodes one must prescribe derivatives instead (Hermite interpolation, group 6).

Fix the nodes and call \(\mathcal I_nf\) the interpolating polynomial of \(f\). Three facts organize everything that follows.

• \(\mathcal I_n\) is linear: \(\mathcal I_n(\alpha f+\beta g)=\alpha\,\mathcal I_nf+\beta\,\mathcal I_ng\), because the system \(Va=y\) is linear in \(y\).
• \(\mathcal I_n\) is a projection: \(\mathcal I_np=p\) for every polynomial \(p\) of degree at most \(n\), hence \(\mathcal I_n^2=\mathcal I_n\).
• \(\mathcal I_nf\) depends on \(f\) only through the \(n+1\) numbers \(f(x_j)\), so it knows nothing about \(f\) between the nodes. The size of \(f-\mathcal I_nf\) there must come from smoothness (group 3), and the sensitivity to errors in the \(y_j\) is the norm of \(\mathcal I_n\), the Lebesgue constant (group 4).

The polynomial itself does not depend on how we write it. The five algorithms of the next group (Lagrange, Newton, Neville, monomial, barycentric) produce the same function, with very different costs and rounding behavior.

Fix the nodes \(x_0,\dots,x_n\) (distinct). For \(j=0,\dots,n\) call \(e_j\) the table that takes the value \(1\) at \(x_j\) and \(0\) at the other \(n\) nodes: \(e_j(x_k)=\delta_{jk}\). Each such table has an interpolating polynomial, and it is unique; call it \(\ell_j\). Every table \(y\) satisfies \(y=\sum_jy_je_j\) (compare the two sides node by node), and the map table \(\mapsto\) interpolating polynomial is linear, so

The point that cannot be overstated is the first one: the Lagrange polynomials are not just a formula that happens to work, they are the interpolating polynomials of the simplest possible tables, built on the same nodes. The general problem has been reduced to these \(n+1\) special problems, and the special problems have an explicit answer.

Look for a polynomial \(q\) of degree at most \(n\) with \(q(x_k)=0\) for \(k\ne j\) and \(q(x_j)=1\). Having the \(n\) roots \(x_k\) (\(k\ne j\)) and degree at most \(n\), it must be \(q(x)=c\prod_{k\ne j}(x-x_k)\) for some constant \(c\): dividing out those \(n\) linear factors leaves a polynomial of degree at most \(0\). The condition \(q(x_j)=1\) gives \(c=1/\prod_{k\ne j}(x_j-x_k)\), that is

So the formula is not guessed, it is derived from the special table. The same root-counting argument proves uniqueness in general (two interpolants differ by a polynomial of degree \(\le n\) with \(n+1\) roots), so it is the only polynomial that works. Try it in the card above: choose the special table \(e_j\) and compare its stems (the values \(0,\dots,1,\dots,0\)) with the curve; the status line below the plot states the largest value of \(|\ell_j|\), which can be much bigger than \(1\) between the nodes.

With \(n\le6\) the card shows the table of the pieces \(y_je_j\): column \(j\) is zero except in row \(j\), where it holds \(y_j\), and the last column, the sum of each row, is the original table \(y\). On the right of the card the curves \(y_j\ell_j\) are the interpolating polynomials of those columns, and the slider adds them one by one: the partial sums are the interpolating polynomials of the partial tables \(y_0e_0+\dots+y_me_m\), and after the last piece the curve passes through every node. In “your data” mode, drag one point up and watch exactly one piece scale, and the sum follow.

Three identities are the same idea read with other tables. The table of all ones is \(e_0+\dots+e_n\) and is interpolated by the constant \(1\), so \(\sum_j\ell_j=1\). The table \((x_0^k,\dots,x_n^k)\) is interpolated by \(x^k\) when \(k\le n\), so \(\sum_jx_j^k\ell_j(x)=x^k\). And since \(\ell_j\) has degree \(n\) while the \(n+1\) functions \(\ell_0,\dots,\ell_n\) are independent (look at their values at the nodes), they form a basis of the polynomials of degree at most \(n\), whose coordinate vector of \(p\) is just the table \((y_0,\dots,y_n)\).

Choose a basis \(b_0,\dots,b_n\) of the polynomials of degree at most \(n\) and write \(p=\sum_ka_kb_k\). The interpolation conditions become the linear system \(Ba=y\) with \(B_{ik}=b_k(x_i)\): the matrix of the basis evaluated at the nodes. The three classical forms are three choices:

The more structure the basis has with respect to the nodes, the less work is left. Lagrange's basis is the extreme case. It also tells you what \(V^{-1}\) is: since \(\ell_j(x)=\sum_ka_{kj}x^k\) is the polynomial of the table \(e_j\), the \(j\)-th column of \(V^{-1}\) holds the monomial coefficients of \(\ell_j\), and the interpolation operator, as a map from tables to coefficient vectors, is multiplication by \(V^{-1}\). Its size in the sup norm is the Lebesgue constant of group 4: \(|p(x)|\le\Lambda_n\max_j|y_j|\) because \(|p(x)|\le\sum_j|y_j||\ell_j(x)|\le(\max_j|y_j|)\sum_j|\ell_j(x)|\).

The name suggests that Lagrange invented the formula, but it is older. Edward Waring published it in 1779, and Euler derived the same formula a few years later (around 1783). Lagrange's own presentation came in his lectures at the École Normale around 1795. What the name does honour is a way of seeing: with Lagrange the formula appears inside a general theory of \(\ell_j\) as the building blocks, the viewpoint of this card, and that is the form in which it spread. A generation later Gauss, in a lecture of 1812 that was published only afterwards, showed that the formulas of Newton, Euler and Waring–Lagrange are one and the same polynomial written in different ways, which is the uniqueness theorem of group 1 in disguise. Newton's own form, older than all of these, is the subject of the next app of this section.

Sources: J. Meijering, A chronology of interpolation: from ancient astronomy to modern signal and image processing, Proc. IEEE 90 (2002) 319–342, for the dates; the exact year of Lagrange's lectures is the one given there, as far as I could check.

Instead of solving a system, build \(p\) directly. For each \(j\) let

Each \(\ell_j\) is a product of \(n\) linear factors, so it has degree \(n\). At a node \(x_i\ne x_j\) the numerator contains the factor \(x_i-x_i=0\); at \(x_j\) every fraction equals 1. Hence \(\ell_j(x_i)=\delta_{ij}\) (the Kronecker delta), and \(p(x_i)=\sum_jy_j\delta_{ij}=y_i\). This is the explicit proof of existence.

The conceptual reading, that \(\ell_j\) is the interpolant of the special table \(e_j\) and that linearity assembles everything else, is the subject of group 1b.

Example. Nodes \(-1,0,1\): \(\ell_0=\tfrac12x(x-1)\), \(\ell_1=1-x^2\), \(\ell_2=\tfrac12x(x+1)\). For the values \(y=(1,0,1)\) we get \(p=\ell_0+\ell_2=x^2\).

Two identities. The constant 1 interpolates the values \((1,\dots,1)\), so by uniqueness \(\sum_j\ell_j(x)=1\) for all \(x\) (the dashed line in the basis card). More generally \(\sum_jx_j^k\ell_j(x)=x^k\) for \(k\le n\).

Assessment. The Lagrange form is wonderful for theory, since the basis \(\ell_j\) is tied to the nodes and the coefficients are the data themselves. As an algorithm it is poor: evaluating \(p(x)\) costs \(O(n^2)\), and adding a node changes every \(\ell_j\).

Newton writes \(p\) in the basis \(1,\ (x-x_0),\ (x-x_0)(x-x_1),\dots\):

The coefficients are the divided differences, defined by

They fill a triangular table (see the Forms card; the coefficients are the highlighted top row). Example. Nodes \(0,1,2,3\), values \(1,2,0,5\): the columns are \(1,2,0,5\); then \(1,-2,5\); then \(-\tfrac32,\tfrac72\); then \(\tfrac53\). So

and indeed \(p(0)=1,\ p(1)=2,\ p(2)=0,\ p(3)=5\) and \(p(1.5)=0.75\).

Why it is useful. The coefficient \(c_k\) uses only the first \(k+1\) nodes, so adding a node adds one new term and leaves the old ones alone (Lagrange cannot do that). Evaluation is nested like Horner's rule, \(p(x)=c_0+(x-x_0)\bigl(c_1+(x-x_1)(c_2+\cdots)\bigr)\), with \(n\) multiplications. The table costs about \(n^2/2\) divisions. Divided differences are symmetric in their arguments, and for smooth \(f\) one has \(f[x_0,\dots,x_k]=f^{(k)}(\xi)/k!\) for some \(\xi\) (a consequence of the error formula of group 3); when all the nodes merge the Newton form becomes the Taylor polynomial.

If you need \(p\) at a single point \(t\), you can skip the coefficients altogether. Let \(P_{i,\dots,i+k}(t)\) be the value at \(t\) of the polynomial that interpolates only the nodes \(x_i,\dots,x_{i+k}\). Starting from \(P_i=y_i\),

Proof. Both sides are polynomials in \(t\) of degree at most \(k\). At \(t=x_i\) the right side gives \(P_{i,\dots,i+k-1}(x_i)=y_i\); at \(t=x_{i+k}\) it gives \(P_{i+1,\dots,i+k}(x_{i+k})=y_{i+k}\); at an inner node both pieces equal \(y\) and the weights \((t-x_i)-(t-x_{i+k})=x_{i+k}-x_i\) add up to the denominator. By uniqueness the two sides coincide.

The cost is \(O(n^2)\) per point and only a vector of storage. The successive columns are the interpolants of increasing degree, which is exactly what Richardson extrapolation and Romberg integration use: extrapolating to \(h=0\) means evaluating an interpolating polynomial at 0.

The most obvious algorithm solves \(Va=y\) by Gaussian elimination. It costs \(O(n^3)\) and, worse, \(V\) is badly conditioned. The condition number \(\kappa_\infty(V)\) (nodes in \([-1,1]\)) grows geometrically:

With 16 digits of arithmetic, a condition number \(\kappa\) costs you about \(\log_{10}\kappa\) digits in the coefficients. The “Barycentric vs Vandermonde” view of the Error card shows the effect: two computations of the same polynomial, whose difference should be zero, differ by \(3.6\cdot10^{-14}\) (n = 10), \(2.9\cdot10^{-11}\) (n = 20) and \(3.1\cdot10^{-8}\) (n = 30) for the Runge function on Chebyshev nodes. For \(e^x\) on equispaced nodes the difference is \(8\cdot10^{-13}\) at n = 20, \(5\cdot10^{-10}\) at n = 30, \(3\cdot10^{-7}\) at n = 40 and \(10^{-4}\) at n = 50.

(The app already helps the Vandermonde route: it rescales the interval to \([-1,1]\) and pivots. On a raw interval like \([0,100]\) it is catastrophically worse.) The moral: the monomial basis is perfect for proofs and for tiny \(n\), and a poor basis for computing.

Put \(\omega(x)=\prod_k(x-x_k)\) and \(w_j=1/\prod_{k\ne j}(x_j-x_k)\) (which equals \(1/\omega'(x_j)\)). Then \(\ell_j(x)=\omega(x)\,\dfrac{w_j}{x-x_j}\), and Lagrange's form becomes the first barycentric form

Interpolating the constant 1 gives \(1=\omega(x)\sum_jw_j/(x-x_j)\). Dividing the first form by this identity eliminates \(\omega\) and gives the second (true) barycentric form

At \(x=x_j\) the formula is \(0/0\) and one simply returns \(y_j\). Since the weights appear in numerator and denominator, any common factor is irrelevant, which is why simple closed forms exist:

• equispaced nodes: \(w_j=(-1)^j\binom nj\);
• Chebyshev points (zeros of \(T_{n+1}\)): \(w_j=(-1)^j\sin\frac{(2j+1)\pi}{2n+2}\);
• Chebyshev–Lobatto points (extrema of \(T_n\), endpoints included): \(w_j=(-1)^j\), halved for \(j=0\) and \(j=n\).

For other nodes the weights cost \(O(n^2)\) once. The app uses the closed forms and checks them against the general formula (agreement to \(5\cdot10^{-15}\)).

Properties. Evaluating costs \(O(n)\) per point; adding a node costs \(O(n)\); the formulas work for any nodes. Higham (2004) proved that the second form is forward stable when the Lebesgue constant of the nodes is small, as for Chebyshev nodes, which makes it the method of choice (it is what Chebfun uses). With equispaced nodes and large \(n\) the interpolation problem itself is ill-conditioned (group 4), and no algorithm can fix that.

The Forms card evaluates the same \(p(x)\) with five algorithms. For small \(n\) the five values agree to 15 or 16 digits; push \(n\) up with equispaced nodes and the Vandermonde value is the first to go wrong.

Theorem. Let \(f\in C^{n+1}[a,b]\), let \(x_0,\dots,x_n\in[a,b]\) be distinct and let \(p\) be the interpolating polynomial. For every \(x\in[a,b]\) there is a point \(\xi\) in the smallest interval containing \(x,x_0,\dots,x_n\) such that

Proof. If \(x\) is a node both sides are zero. Otherwise fix \(x\), and set \(c=\dfrac{f(x)-p(x)}{\omega(x)}\) and \(g(t)=f(t)-p(t)-c\,\omega(t)\). Then \(g\) vanishes at the \(n+1\) nodes (there \(f-p=0\) and \(\omega=0\)) and at \(t=x\) by the choice of \(c\): that is \(n+2\) distinct zeros. By Rolle's theorem \(g'\) has at least \(n+1\) zeros, \(g''\) at least \(n\), and so on until \(g^{(n+1)}\) has at least one zero \(\xi\). But \(p^{(n+1)}=0\) and \(\omega^{(n+1)}=(n+1)!\), so

Corollaries. \(\|f-p\|_\infty\le\dfrac{\|f^{(n+1)}\|_\infty}{(n+1)!}\,\|\omega\|_\infty\). If \(f\) is a polynomial of degree at most \(n\), the error is zero. For linear interpolation (\(n=1\), spacing \(h\)) the maximum of \(|\omega|\) is \(h^2/4\), so the error is at most \(h^2\|f''\|/8\). There is also a version that needs no derivatives: \(f(x)-p(x)=\omega(x)\,f[x_0,\dots,x_n,x]\), which is the Newton form with one more node.

The error is a product of two things with different owners. The factor \(f^{(n+1)}(\xi)/(n+1)!\) belongs to the function; we cannot change it. The factor \(\omega(x)\) belongs to the nodes; we choose it. Use the “Node polynomial” view of the Error card.

On \([-1,1]\) the largest \(|\omega|\) is

At \(n=20\) the Chebyshev \(\omega\) is about 245 times smaller. For equispaced nodes \(\omega\) is tiny in the middle and large near the two ends, so the error piles up there; for Chebyshev nodes the \(n+2\) lobes of \(|\omega|\) all have the same height.

The derivative factor can also be dangerous. For \(1/(1+25x^2)\) the Taylor coefficients at 0 are \(\pm25^k\), so \(f^{(n+1)}(0)/(n+1)!\) is about \(5^{n+1}\) in size. The factor \(\omega\) of equispaced nodes shrinks only by about 0.7 per step (8.5·10⁻³ at \(n=10\), 2.3·10⁻⁴ at \(n=20\)), so the bound grows with \(n\), and so does the actual error (group 3, Runge).

Take \(f=e^x\) on \([-1,1]\), where \(\|f^{(n+1)}\|_\infty=e\). With Chebyshev nodes \(\|\omega\|=2^{-n}\), so

The bound holds and is sharp within a factor of about 2.4 for the Chebyshev nodes. Use it to choose \(n\): to guarantee error below \(10^{-6}\) we need \(e/((n+1)!2^n)\lt 10^{-6}\), which fails for \(n=6\) (\(8.4\cdot10^{-6}\)) and holds for \(n=7\) (\(5.3\cdot10^{-7}\)). Since \(e^x\) is entire, the error falls faster than any geometric rate until it reaches rounding level, about \(10^{-15}\) from \(n=15\) on.

Choose “Runge” with equispaced nodes. At \(n=5\) the polynomial is a poor but sane approximation; at \(n=10\) it swings far above the function near the ends; at \(n=20\) it is off the chart. The maximum errors on \([-1,1]\) for \(1/(1+25x^2)\):

At \(x=0.95\) and \(n=10\) the polynomial has the value 1.92 while \(f(0.95)=0.042\). More points made it worse. Runge found this in 1901: the interpolants converge to \(f\) only for \(|x|\lt 0.726\) and diverge for \(|x|>0.726\) (approximately).

The reason lies in the complex plane: \(1/(1+25z^2)\) has poles at \(\pm i/5\), close to the interval. Equispaced interpolation converges only for functions analytic in a lens-shaped region around \([-1,1]\) that crosses the imaginary axis at about \(\pm0.53\,i\); the poles at \(\pm0.2\,i\) are inside it. The same story for “1/(1 + x²) on [−5, 5]”: it is the same function after rescaling, and the numbers are identical, which shows that only the relative position of the singularity matters. Chebyshev nodes cure it (next lesson, and group 5).

Since \(\omega\) is the part we control, choose the nodes to make \(\|\omega\|_\infty\) as small as possible. Theorem (Chebyshev). Among all monic polynomials of degree \(n+1\), the smallest sup-norm on \([-1,1]\) belongs to \(T_{n+1}(x)/2^n\), where \(T_{n+1}(x)=\cos((n+1)\arccos x)\); its norm is \(2^{-n}\). So the best nodes are the zeros of \(T_{n+1}\),

Proof. \(T_{n+1}\) has leading coefficient \(2^n\), so \(\omega^\ast=T_{n+1}/2^n\) is monic, and it alternates between \(+2^{-n}\) and \(-2^{-n}\) at the \(n+2\) points \(\cos(k\pi/(n+1))\). If a monic \(q\) had \(\|q\|_\infty\lt 2^{-n}\), then \(r=\omega^\ast-q\) would have degree at most \(n\) and, at those \(n+2\) points, the signs of \(\omega^\ast\): \(n+1\) sign changes, hence \(n+1\) zeros, hence \(r=0\), contradicting \(\|q\|\lt 2^{-n}\). \(\blacksquare\)

On \([a,b]\) the nodes are mapped affinely and the bound of the previous lesson becomes

The Chebyshev–Lobatto points (extrema of \(T_n\), endpoints included) give an \(\omega\) of nearly the same size and are often more convenient. Measured on the Runge function, maximum errors for Chebyshev nodes at \(n=10,20,40\) are \(0.109,\ 0.0153,\ 2.9\cdot10^{-4}\); for Lobatto \(0.132,\ 0.0177,\ 3.4\cdot10^{-4}\). One curiosity: the Chebyshev error at \(n=15\) (0.083) is larger than at \(n=12\) (0.069). For odd \(n\) there is an even number of nodes and none sits at \(x=0\), exactly where the peak of the function is.

Now perturb the data. If every \(y_j\) changes by at most \(\delta\), the interpolant changes by \(\sum_j\delta_j\ell_j(x)\), whose size is at most \(\delta\lambda_n(x)\), where

The bound is attained by taking \(\delta_j=\delta\,\mathrm{sign}\,\ell_j(x)\). So \(\Lambda_n\), the Lebesgue constant, is exactly the worst-case amplification of data errors: the condition number of interpolation. It depends only on the nodes. See the Lebesgue function in the basis card. Measured on \([-1,1]\):

In the app, take \(\sin(2\pi x)\), \(n=20\), and add noise \(10^{-3}\) to the values: with equispaced nodes the maximum error becomes about 1.8, with Chebyshev nodes about \(1.2\cdot10^{-3}\), the size of the noise itself.

Let \(E_n(f)=\min_{q\in\mathcal P_n}\|f-q\|_\infty\) be the error of the best polynomial of degree at most \(n\) (it exists and is unique, by Chebyshev's equioscillation theorem). Interpolation is never much worse:

Proof. Let \(q\) be any polynomial of degree at most \(n\). Because \(\mathcal I_nq=q\) and \(\mathcal I_n\) is linear, \(f-\mathcal I_nf=(f-q)-\mathcal I_n(f-q)\), so \(\|f-\mathcal I_nf\|\le(1+\|\mathcal I_n\|)\|f-q\|\). Take \(q\) to be the best polynomial. \(\blacksquare\)

So two things control interpolation: how well polynomials can approximate \(f\) at all (\(E_n(f)\), a question of smoothness, next group) and how good the nodes are (\(\Lambda_n\)). With Chebyshev nodes \(1+\Lambda_n\) is 3.9 at \(n=40\) and about 5.4 even at \(n=400\): you lose at most a small factor against the very best polynomial, and you get it with an \(O(n)\) formula and no iteration. With equispaced nodes the factor \(1+\Lambda_n\) is astronomical and the guarantee is worthless.

The asymptotics, all visible in the “Lebesgue constant against n” plot:

• equispaced nodes: \(\Lambda_n\approx\dfrac{2^{n+1}}{e\,n\ln n}\) (about \(1.1\cdot10^{4}\) measured at \(n=20\), \(4.7\cdot10^{9}\) at \(n=40\)): exponential growth, roughly doubling at each step;
• Chebyshev–Lobatto nodes: \(\Lambda_n\approx\dfrac2\pi\ln n+0.9625\), checked numerically to four digits at \(n=100\) (3.894) and \(n=400\) (4.777); Chebyshev zeros are the same up to a constant: logarithmic growth.

Is there a perfect set of nodes with bounded \(\Lambda_n\)? No. Erdős (1961) proved that for every choice of nodes \(\Lambda_n\ge\frac2\pi\ln(n+1)-C\) for a constant \(C\). So Chebyshev nodes are optimal up to a bounded additive constant, the best one can hope for.

Look at the Lebesgue function \(\lambda_n(x)\) for equispaced nodes at \(n=20\): it is of order 1 to 10 in the middle of the interval and shoots up to about \(10^4\) in the first and last subintervals. Nodes too far apart near the ends are what makes a polynomial through them swing wildly there; Chebyshev nodes cluster near the ends with density \(1/(\pi\sqrt{1-x^2})\) to prevent exactly this.

Even if \(f\) is evaluated perfectly, each \(y_j\) is stored with relative error up to about \(\varepsilon\approx1.1\cdot10^{-16}\). That is a perturbation of the data, and interpolation can multiply it by \(\Lambda_n\). Two measurements:

• \(f=e^x\), equispaced, \(n=40\): \(\Lambda_n=4.7\cdot10^9\), expected damage \(\sim5\cdot10^{-7}\), measured error \(3.3\cdot10^{-7}\), although the theoretical interpolation error is far below \(10^{-30}\) (the function is entire).
• The polynomial \(x^5-2x^3+x/2\) is reproduced exactly for every \(n\ge5\) in exact arithmetic. In floating point, equispaced nodes give errors \(1.1\cdot10^{-15}\) (n = 10), \(5.5\cdot10^{-13}\) (n = 20), \(2.4\cdot10^{-10}\) (n = 30) and \(1.6\cdot10^{-7}\) (n = 40), while Chebyshev nodes stay below \(10^{-15}\) throughout.

This is a property of the problem, not of the formula: it afflicts barycentric, Newton and Lagrange alike. For large \(n\), equispaced interpolation is therefore doubly wrong, inaccurate for functions like Runge's and noisy even for polynomials.

The tempting guess is that \(\mathcal I_nf\to f\) uniformly for every continuous \(f\) if the nodes are decent. It is false. Faber's theorem (1914): for every triangular array of nodes there is a continuous function \(f\) for which \(\|f-\mathcal I_nf\|_\infty\) does not tend to 0. The reason is the uniform boundedness principle: \(\Lambda_n=\|\mathcal I_n\|\to\infty\) for every array (Erdős's bound of group 4), so some \(f\) is hit by the growth.

For equispaced nodes things are worse than “some function”: Bernstein (1918) showed that for \(f(x)=|x|\) the interpolants diverge at every point of \([-1,1]\) except \(-1,0,1\). Try “|x|” with equispaced nodes: the maximum errors are 0.66, 95 and \(1.5\cdot10^7\) for \(n=10,20,40\). With Chebyshev nodes they are 0.055, 0.029 and 0.015: convergence, slowly.

The positive side: with Chebyshev nodes \(\Lambda_n=O(\ln n)\), so \(\|f-\mathcal I_nf\|\le(1+\Lambda_n)E_n(f)\) tends to 0 whenever \(E_n(f)\ln n\to0\), in particular for every Lipschitz function, and every smoother one. (Marcinkiewicz showed that for each continuous \(f\) some node array gives convergence, but the array depends on \(f\), which is of no practical use.)

By the near-best property, the speed is that of the best approximation \(E_n(f)\), up to a factor \(1+\Lambda_n\). And \(E_n(f)\) is governed by smoothness:

• (Jackson) if \(f\in C^k[-1,1]\) then \(E_n(f)\le C_k\,\|f^{(k)}\|_\infty/n^k\): a power law, the faster the smoother;
• if \(f\) is analytic in a neighborhood of \([-1,1]\), then \(E_n(f)\) decays geometrically, \(O(\rho^{-n})\) (next lesson);
• a corner costs one power: for \(f=|x|\), \(nE_n(f)\to0.2802\) (Bernstein's constant). Chebyshev interpolation gives \(0.0546,\ 0.0285,\ 0.0146\) at \(n=10,20,40\), about twice \(0.28/n\), as the factor \(1+\Lambda_n\approx3\) allows.

The conclusion is practical. The accuracy you can reach with a polynomial of degree \(n\) is a property of the function, and Chebyshev interpolation gets within a small factor of it. Nothing smarter can beat \(E_n\); it can only avoid the factor \(1+\Lambda_n\).

Let \(E_\rho\) be the ellipse with foci \(\pm1\) and semi-axes \((\rho+\rho^{-1})/2\) and \((\rho-\rho^{-1})/2\). Theorem (Bernstein; see Trefethen). If \(f\) is analytic inside \(E_\rho\) with \(|f|\le M\) there, then the interpolant at \(n+1\) Chebyshev points satisfies

Geometric convergence, at a rate fixed by the nearest singularity. Check with Runge's function. Its poles at \(\pm i/5\) lie on \(E_\rho\) when \((\rho-\rho^{-1})/2=1/5\), that is \(\rho=\frac{1}{5}+\sqrt{1+\frac1{25}}=1.2198\). The prediction: the error shrinks by \(\rho^{-10}=0.137\) every 10 steps of \(n\) and by \(\rho^{-20}=0.0188\) every 20. Measured with Chebyshev nodes: from \(n=10\) to 20 the factor is 0.140, from \(n=20\) to 40 it is 0.0189. The theory predicts the slope.

The rate shows up as a straight line on the logarithmic convergence plot. For \(\sin(2\pi x)\) and \(e^x\), entire functions with no singularity at all, the line bends downward faster than any straight line (superexponential convergence) until it meets the rounding floor near \(10^{-15}\). The same theory explains equispaced failure: those nodes behave as if singularities near the imaginary axis (the poles of Runge's function) were far closer than they are.

Maximum error with Chebyshev nodes, measured in the app (choose the function and set n):

Reading the rows: for the jump, nothing helps; for the corner and the endpoint singularity the error halves when \(n\) doubles (\(O(1/n)\)); for the bump it falls faster than any power of \(1/n\) but slower than any geometric rate; for Runge it falls geometrically at rate \(1.22^{-n}\); for entire functions it falls superexponentially until rounding. The “Maximum error against n” plot shows exactly these shapes.

Put \(x=\cos\theta\). The Chebyshev points are \(x_j=\cos\theta_j\) with equally spaced angles \(\theta_j=(2j+1)\pi/(2n+2)\): they are the shadows, on the diameter, of equally spaced points on the unit circle. A polynomial of degree \(n\) in \(x\) is the same thing as an even trigonometric polynomial \(\sum_{k\le n}a_k\cos k\theta\), because \(T_k(\cos\theta)=\cos k\theta\).

So interpolating in Chebyshev points is interpolating a trigonometric polynomial at equispaced angles, and that is the discrete Fourier transform: orthogonal, perfectly conditioned, and the source of the \(\frac2\pi\ln n\) growth of \(\Lambda_n\). Points equispaced in \(x\) are not equispaced in \(\theta\), and that is the whole difference between Chebyshev and Runge.

Practical bonus: the Chebyshev coefficients \(a_k\) of the interpolant can be computed in \(O(n\log n)\) by the fast Fourier transform, which is the basis of spectral methods and of the Chebfun system.

Replace \(f\) by \(p\) and integrate: \(\displaystyle\int_a^bf\approx\int_a^bp=\sum_jW_jy_j\) with \(W_j=\int_a^b\ell_j\). This is how every interpolatory quadrature rule is born. For equispaced nodes on \([0,n]\) (spacing 1):

For \(n=8\) the weights are \(0.279,\ 1.662,\ -0.262,\ 2.962,\ -1.281,\ 2.962,\ -0.262,\ 1.662,\ 0.279\): some are negative, so \(\sum|W_j|>b-a\) and rounding errors in \(y_j\) are amplified, the same disease as the Lebesgue constant. That is why high-order Newton–Cotes rules are abandoned. Remedies: composite rules (low degree on many small panels), Clenshaw–Curtis (integrate the Chebyshev interpolant), and Gauss quadrature, where the nodes are chosen to make the rule exact up to degree \(2n+1\).

Differentiating the error formula gives rules for \(f'\): from \(n=2\) one gets the central difference \((f(x+h)-f(x-h))/(2h)\) with error \(-h^2f'''(\xi)/6\). Differentiation is intrinsically ill-posed: the norm of the Chebyshev differentiation matrix grows like \(n^2\), so noise is amplified much more than in interpolation. Adams–Bashforth multistep methods for ODEs integrate the interpolant of \(f\) through the last few steps.

Prescribe values and first derivatives at \(n+1\) nodes: \(2n+2\) conditions, and exactly one polynomial \(H\) of degree at most \(2n+1\) satisfies them. The proof of the error formula goes through verbatim (each node is a double zero of \(f-H\)) and gives

Hermite interpolation is computed with the Newton table by repeating nodes and putting \(f[x_j,x_j]=f'(x_j)\) in the slots where a division by zero would occur. For \(n=1\) it is the cubic Hermite interpolant, the building block of splines in computer graphics and of monotone “PCHIP” interpolation.

The extreme case is a single node repeated \(n+1\) times: the Newton form becomes the Taylor polynomial \(\sum_kf^{(k)}(x_0)(x-x_0)^k/k!\), and the error formula becomes Lagrange's remainder \(f^{(n+1)}(\xi)(x-x_0)^{n+1}/(n+1)!\). So Taylor's theorem is the special case of this lesson in which all the nodes collide, and interpolation is its natural generalization to separated nodes, which is why it is far better behaved on a whole interval.

A single polynomial of high degree is a global object: a bump at one end can shake the whole curve. The alternative is low degree on many small pieces.

• Piecewise linear interpolation is continuous with error at most \(h^2\|f''\|_\infty/8\), from the error formula with \(n=1\). Stable (\(\Lambda=1\)), but only first-order smooth.
• Cubic splines are piecewise cubics with two continuous derivatives. For smooth \(f\) the error is \(O(h^4)\), and the Lebesgue constant stays bounded as the mesh is refined, unlike global polynomial interpolation.

A working rule: if \(f\) is analytic (or very smooth) on the whole interval and you may choose the nodes, global interpolation at Chebyshev points converges geometrically and wins. If the data are given at arbitrary points, or \(f\) has corners and jumps, or you only need a few digits, low-degree pieces are the safe choice. In the free-data mode, put twelve points on a gentle zigzag: the global polynomial oscillates between and beyond them, which a spline would not do.

If the \(y_j\) are measurements with errors of size \(\sigma\), interpolation reproduces the errors faithfully, and can amplify them up to \(\Lambda_n\sigma\) in between. With equispaced data and \(n\) around 20 this is a catastrophe (use the noise control with “equispaced”, then switch the node family to see the contrast).

What one does in practice is least-squares fitting: choose a polynomial of low degree \(m\ll n\) minimizing \(\sum_j(p(x_j)-y_j)^2\), or a smoothing spline. The fit does not pass through the points, and that is a feature: it averages the noise instead of chasing it. Interpolation is the right tool when the data are exact (values of a known function, or nodes you choose) and the right question is “what is between the nodes?”.

1. Choose “Runge: 1/(1 + 25x²)”, equispaced nodes, and press Grow n. The maximum error first stays near 0.4 (0.44 at \(n=4\), 0.43 at \(n=5\)) and then explodes: 0.43 at \(n=5\), 1.92 at \(n=10\), 59.8 at \(n=20\), \(2.4\cdot10^3\) at \(n=30\), \(1.05\cdot10^5\) at \(n=40\).
2. Tick “Fit the window to p” to see the whole polynomial, and open the Lebesgue-function view of the basis card.
3. Switch to Chebyshev: \(0.109\), \(0.0153\), \(2.9\cdot10^{-4}\) at \(n=10,20,40\). In the convergence card the Chebyshev curve is a straight line of slope \(\log_{10}1.22^{-1}\) per unit of \(n\).
4. Choose “1/(1 + x²) on [−5, 5]” and compare: same function after rescaling, identical numbers.

Question: why does the Chebyshev error at \(n=15\) exceed the one at \(n=12\)? (Lesson on Chebyshev nodes.)

With Chebyshev nodes, run through “|x|”, “√(x + 1)”, “smooth bump”, “Runge” and “sin(2πx)” at \(n=10,20,40\) and compare with the table of the catalogue lesson. In the convergence card, notice the shapes: a flat-ish line for the corner (\(O(1/n)\) is a slow straight line only on a log–log plot, here it looks like a gentle bend), a straight line for Runge, a curve falling away for \(\sin\). Then switch to equispaced nodes for \(|x|\): the curve goes up after \(n\approx5\), which is Bernstein's divergence.

Question: the table suggests doubling \(n\) halves the error for \(|x|\). Which part of Jackson's theorem is that?

Choose “polynomial x⁵ − 2x³ + x/2”. For \(n\ge5\) the error is rounding only (about \(10^{-15}\)): uniqueness. For \(n=4\) it is not zero, and it can be predicted exactly. The function is monic of degree 5, so \(f^{(5)}/5!=1\) and the error formula says \(f-p=\omega(x)\) with no \(\xi\) left over. Therefore the maximum error equals \(\|\omega\|_\infty\): for equispaced nodes 0.113, and for Chebyshev nodes exactly \(2^{-4}=0.0625\), which is what the app reports. Check both against the “Node polynomial” view of the Error card.

Choose \(\sin(2\pi x)\), \(n=20\), noise \(10^{-3}\). Equispaced nodes: maximum error about 1.8, a thousand times the noise, close to the \(\Lambda_{20}\approx1.1\cdot10^4\) allowance. Chebyshev nodes: about \(1.2\cdot10^{-3}\), the size of the noise. Look at the polynomial with “Fit the window to p”: with equispaced nodes the damage sits in the first and last subintervals, exactly where the Lebesgue function is large. Try noise \(10^{-6}\) and \(10^{-2}\) and verify that the damage scales linearly with the noise.

Choose “Your data”. Drag a point and watch the whole polynomial move: it is a global object. Add points to the right of the last one: the polynomial far to the left changes too. Put two points very close together with different heights: the polynomial needs a huge slope between them and \(\Lambda_n\) leaps. Press Add a random point repeatedly and watch the Vandermonde condition number in the metrics. In the Forms card (for \(n\le8\)) follow the divided-difference table as you drag: the top row, the Newton coefficients, changes only from the dragged point on to the right.

1. Choose “Your data” and place four points. In the card Lagrange's idea pick \(e_0\) and look at the left plot: the polynomial is \(1\) at the first point and \(0\) at the other three. Do the same with \(e_1,e_2,e_3\).
2. Drag the slider on the right from \(0\) to the last piece. Which piece is largest? Drag the corresponding point upward and watch that piece, and only that piece, grow.
3. Choose “Your function f(x)”, type sin(3x) on \([-1,1]\) with \(n=3\) equispaced nodes, and compare the sum of the four pieces with \(f\).
4. With the Python script below, change NODES to "cheb" and compare the pictures lagrange_pieces.png for the two families: the pieces \(y_j\ell_j\) are much smaller with Chebyshev nodes, which is the Lebesgue constant at work.

1. Find the polynomial through \((-1,1),(0,0),(1,1)\) by Lagrange's formula and by uniqueness. Answer: \(\ell_0+\ell_2=x^2\); and \(x^2\) has degree 2 and fits, so by uniqueness it is the one.
2. Show that \(\sum_j\ell_j(x)\equiv1\) and \(\sum_jx_j^k\ell_j(x)=x^k\) for \(k\le n\). Answer: \(1\) and \(x^k\) are polynomials of degree at most \(n\), so each equals its own interpolant, whose Lagrange form is the stated sum.
3. Show that \(w_j=1/\omega'(x_j)\). Answer: write \(\omega=(x-x_j)\prod_{k\ne j}(x-x_k)\); then \(\omega'(x_j)=\prod_{k\ne j}(x_j-x_k)\).
4. Bound the error of the linear interpolation of \(\ln x\) at the nodes 1 and 2, at \(x=1.5\). Answer: \(|f''|=1/x^2\le1\), \(|\omega(1.5)|=0.25\), so the error is at most \(\tfrac12\cdot0.25=0.125\). The actual error is \(\ln1.5-\tfrac12\ln2=0.0589\).
5. How many Chebyshev nodes guarantee error below \(10^{-6}\) for \(e^x\) on \([-1,1]\)? Answer: \(e/((n+1)!2^n)\lt 10^{-6}\) first holds at \(n=7\), that is 8 nodes.
6. Derive Simpson's rule from \(\int\ell_j\). Answer: on nodes \(0,1,2\): \(\ell_0=\tfrac12(x-1)(x-2)\), \(\ell_1=x(2-x)\), \(\ell_2=\tfrac12x(x-1)\), with integrals over \([0,2]\) equal to \(\tfrac13,\tfrac43,\tfrac13\). So \(\int_0^2f\approx\tfrac13(f_0+4f_1+f_2)\).
7. Prove that the Newton coefficient \(c_k\) depends only on \(x_0,\dots,x_k\) and \(y_0,\dots,y_k\). Answer: \(c_k\) is the leading coefficient of the polynomial of degree at most \(k\) interpolating the first \(k+1\) data, since the Newton terms \(c_j\prod_{i\lt j}(x-x_i)\) with \(j\lt k\) have degree below \(k\).

To integrate \(f\) on \([a,b]\), choose nodes \(x_0,\dots,x_n\), let \(p_n\) be the polynomial of degree \(\le n\) through the points \((x_j,f(x_j))\), and take \(\int_a^b p_n\) as the answer. In Lagrange form \(p_n=\sum f(x_j)\ell_j\), so

The weights depend only on the nodes, never on \(f\). A closed Newton–Cotes rule takes equally spaced nodes \(x_j=a+jh\), \(h=(b-a)/n\), endpoints included. Substituting \(x=a+th\) gives \(w_j=h\,c_j\) with \(c_j=\int_0^n\prod_{k\ne j}\frac{t-k}{j-k}\,dt\), a number you compute once. The panel “Closed Newton–Cotes weights” lists them: \(n=1\) gives \(\tfrac12,\tfrac12\) (trapezoid), \(n=2\) gives \(\tfrac13,\tfrac43,\tfrac13\) (Simpson 1/3), \(n=3\) gives \(\tfrac38,\tfrac98,\tfrac98,\tfrac38\) (Simpson 3/8).

Connection with the previous app: this is exactly “integrate the Lagrange interpolant”. The error of the rule is the integral of the interpolation error, and that is how every error term below is obtained.

A rule has degree of precision \(d\) if it is exact for every polynomial of degree \(\le d\) and not for some polynomial of degree \(d+1\). Testing on monomials \(1,x,x^2,\dots\) is enough, by linearity. The panel “Degree of precision” does this on \([0,1]\).

An interpolatory rule with \(n+1\) nodes has \(d\ge n\), because \(p_n=f\) whenever \(\deg f\le n\). For even \(n\) there is one more degree: put \(m\) at the midpoint of \([a,b]\); the nodes are symmetric about \(m\) and so are the weights, so both \(\int (x-m)^{n+1}\) and the rule applied to it vanish (odd function, symmetric sum). Hence

trapezoid \(n=1\): \(d=1\); Simpson 1/3 \(n=2\): \(d=3\), not 2; Simpson 3/8 \(n=3\): \(d=3\).

In the app, choose \(x^3\) on \([0,1]\), Simpson 1/3, simple: the error is zero (below \(10^{-16}\)), even though the parabola through three points of a cubic is not the cubic. For \(x^4\) the simple errors are \(-1/120\approx-8.3\cdot10^{-3}\) (Simpson 1/3) and \(-1/270\approx-3.7\cdot10^{-3}\) (Simpson 3/8); they will fix the constants of the error terms.

With \(n=1\) the interpolant is the chord, and its integral is the area of a trapezoid:

Choose “Trapezoid”, “Simple”, function \(\sin x\) on \([0,\pi]\): the chord is the \(x\)-axis, the rule gives \(0\), the integral is \(2\), the error is \(2\). A single panel is a crude tool; the point of the next lessons is to quantify how crude, and to cure it by cutting \([a,b]\) into many panels.

Let \(f\in C^2[a,b]\). The interpolation error through \(a,b\) is \(f(x)-p_1(x)=\tfrac12 f''(\xi_x)(x-a)(x-b)\). Integrate:

The factor \((x-a)(x-b)\) is \(\le0\) on all of \([a,b]\): it does not change sign. So the mean value theorem for integrals pulls \(f''\) out, at some \(\xi\in(a,b)\):

Reading it: convex \(f\) (\(f''>0\)) gives \(E<0\), the chord lies above the graph and the rule overestimates. Example in the app: \(e^x\) on \([0,1]\), simple trapezoid. The rule gives \(1.8591\), \(E=-0.14086\); the bound is \(e/12=0.2265\); solving \(-\tfrac1{12}e^{\xi}=E\) gives \(\xi=0.525\), inside \((0,1)\), as the theorem promises. The “Error term” card recomputes all this for the case you choose.

Cut \([a,b]\) into \(n\) panels of width \(h=(b-a)/n\) and apply the simple rule on each. The interior nodes are counted twice and the ends once:

Each panel contributes \(-\tfrac{h^3}{12}f''(\xi_i)\). Summing, \(-\tfrac{h^3}{12}\sum_{i=1}^n f''(\xi_i)\). The average \(\frac1n\sum f''(\xi_i)\) lies between the min and the max of the continuous \(f''\), so by the intermediate value theorem it equals \(f''(\xi)\) for some \(\xi\), and with \(nh=b-a\):

The error is of order 2: halving \(h\) divides it by 4. Check in the app with \(\sin x\) on \([0,\pi]\): for \(n=6,12,24,48\) the errors are \(4.59\cdot10^{-2}\), \(1.144\cdot10^{-2}\), \(2.857\cdot10^{-3}\), \(7.14\cdot10^{-4}\), and the convergence table shows order \(2.00\). On the log–log plot it is a straight line of slope \(-2\).

The mean-value form hides the true behavior. For smooth \(f\) the error has an expansion in even powers of \(h\) (Euler–Maclaurin):

Only the endpoints appear. Two consequences. First, the leading term is \(-\tfrac{h^2}{12}(f'(b)-f'(a))\), which is the integral of \(f''\) in disguise, so the true error is usually a bit smaller than the bound with \(\max|f''|\). For \(\sin x\) on \([0,\pi]\) the leading term is \(h^2/6\), which is the observed error to three digits.

Second, if \(f\) is periodic (or just has matching odd derivatives at \(a\) and \(b\)) every term vanishes and the trapezoid rule converges faster than any power of \(h\). For \(e^{\cos x}\) on \([0,2\pi]\) the error is \(3.4\cdot10^{-2}\) for \(n=4\), \(2.8\cdot10^{-4}\) for \(n=6\), \(1.3\cdot10^{-6}\) for \(n=8\), \(6.5\cdot10^{-12}\) for \(n=12\): no straight line of slope \(-2\) here. This is why the trapezoid rule is the method of choice for periodic integrands, and it prepares the ground for Gaussian quadrature, the next topic.

Since \(I-T_n=c_2h^2+c_4h^4+\cdots\) with the same constants for \(h\) and \(h/2\), the combination that kills the \(h^2\) term is

Writing it out, the weights on the \(2n+1\) points are \(\tfrac13(1,4,2,4,\dots,4,1)\) times the finer step: this is the composite Simpson 1/3 rule. Numerically, for \(\sin x\) on \([0,\pi]\): \(T_6\) and \(T_{12}\) give \((4T_{12}-T_6)/3=2.0000526243\), identical to Simpson 1/3 with \(n=12\). Repeating the idea on the \(h^4\) term gives Romberg integration. The lesson: the trapezoid rule and Simpson's rule are one family, and the extra order of Simpson comes from cancelling the leading term of the Euler–Maclaurin series.

Take three nodes \(m-h,\,m,\,m+h\) (with \(m=\tfrac{a+b}2\), \(h=\tfrac{b-a}2\)) and the parabola through them. Instead of integrating Lagrange polynomials, use the degree of precision: the rule \(w_0f_0+w_1f_1+w_2f_2\) must be exact on \(1\) and \((x-m)^2\), and by symmetry \(w_0=w_2\). Then

So

On \(\sin x\) over \([0,\pi]\) the simple rule gives \(2.0944\) (error \(-0.0944\)) against \(0\) for the simple trapezoid: three points and a huge gain. Geometrically, we integrate a parabola through the midpoint instead of a chord.

The Lagrange error gives \(f'''\), but the rule has degree 3, so we expect \(f^{(4)}\). The trick: let \(q\) be the cubic with \(q=f\) at \(m-h,m,m+h\) and \(q'(m)=f'(m)\) (Hermite). The rule only sees the three values, so \(Q(f)=Q(q)=\int q\), because Simpson is exact on cubics. Hence \(E=\int(f-q)\), and the Hermite error is

With \(t=x-m\), the weight \((t^2-h^2)t^2\) is \(\le0\) on \([-h,h]\): constant sign again, so

Check against the \(x^4\) test: \(f^{(4)}=24\), \(h=\tfrac12\): \(-\tfrac{1}{32\cdot90}\cdot24=-\tfrac1{120}\), the value in the exactness table. For \(e^x\) on \([0,1]\) the bound is \(\tfrac{e}{90\cdot32}=9.4\cdot10^{-4}\) and the actual error is \(-5.8\cdot10^{-4}\).

Take \(n\) even and apply the simple rule on each of the \(n/2\) pairs of subintervals (each of width \(2h\)):

Each pair has error \(-\tfrac{h^5}{90}f^{(4)}(\xi_i)\); adding \(n/2\) of them and using the intermediate value theorem as before:

Order 4: halving \(h\) divides the error by 16. In the app, \(e^{-x^2}\) on \([0,2]\) with \(n=6,12\): \(4.98\cdot10^{-5}\) then \(3.14\cdot10^{-6}\), ratio \(15.9\). And the cost comparison is striking: for \(\sin x\) on \([0,\pi]\) with \(n=12\), trapezoid and Simpson use the same 13 points, but the errors are \(1.1\cdot10^{-2}\) and \(5.3\cdot10^{-5}\), a factor of 200.

Four nodes, spacing \(h=\tfrac{b-a}3\), cubic interpolant. Symmetry gives \(w_0=w_3\), \(w_1=w_2\). With \(m\) the midpoint the nodes are \(m\pm\tfrac{3h}2,\ m\pm\tfrac h2\). Exactness on \(1\) and \((x-m)^2\):

Its degree of precision is 3 (not 4: the symmetry bonus applies only to even \(n\)), so \(E=K\,h^5f^{(4)}(\xi)\) with a constant of sign \(-\) (the Peano kernel of the rule has constant sign, a fact we quote). The constant is fixed by the test \(f=x^4\): on \([0,1]\) the table gives \(E=-\tfrac1{270}\), with \(h=\tfrac13\), \(f^{(4)}=24\), so \(K\,\tfrac1{243}\cdot24=-\tfrac1{270}\), which gives

The composite rule needs \(n\) a multiple of 3, has coefficients \(1,3,3,2,3,3,2,\dots,3,3,1\), and error \(-\tfrac{(b-a)h^4}{80}f^{(4)}(\xi)\) (there are \(n/3\) groups, each contributing \(-\tfrac{3h^5}{80}f^{(4)}\)).

Both rules are of order 4; the constants of the composite errors are \(\tfrac1{180}\) and \(\tfrac1{80}\). At the same \(n\) (hence the same \(h\) and the same number of evaluations) the 3/8 error is \(\tfrac{180}{80}=2.25\) times larger. Observed with \(\sin x\) on \([0,\pi]\), \(n=12\): \(5.26\cdot10^{-5}\) for 1/3 against \(1.19\cdot10^{-4}\) for 3/8, ratio \(2.27\).

So when \(n\) is even, use 1/3. The 3/8 rule earns its keep when \(n\) is odd: take \(n-3\) (even) subintervals with Simpson 1/3 and the last three with Simpson 3/8; the combined rule has the same order. Example: \(n=9\) is one 3/8 plus three 1/3 panels. The two rules are also the first members of a pattern: Boole's rule (\(n=4\), degree 5) is the next with even \(n\).

To guarantee \(|E|\le\varepsilon\) solve the bound for \(n\). Take \(\int_0^\pi\sin x\,dx\) with \(\varepsilon=10^{-8}\) and \(\max|f''|=\max|f^{(4)}|=1\):

Direct experiments show how pessimistic the bound is: the trapezoid rule reaches \(10^{-8}\) with about \(12\,830\) subintervals (the Euler–Maclaurin term \(h^2/6\) predicts \(\pi/\sqrt{6\cdot10^{-8}}\approx12\,826\)) and Simpson with \(102\). The ratio \(16\,075/114\) already says why Simpson is the default for smooth integrands. In practice one also does the reverse: double \(n\) until two successive values agree to the tolerance.

The error terms assume \(f\in C^2\) or \(C^4\). Test the app on two functions that violate this.

\(\sqrt x\) on \([0,1]\): \(f''\) blows up at 0, so the bound is infinite. Measured orders: \(1.46\) to \(1.49\) for the trapezoid, \(1.50\) for both Simpson rules. The singularity \(x^{1/2}\) at the end limits every rule to order \(3/2\), and Simpson's extra accuracy is wasted. Remedy: substitute \(x=t^2\) to remove the singularity, or split the interval.

\(|x-1/\pi|\) on \([0,1]\): a corner inside the interval. The error is irregular because it depends on where the corner falls relative to the nodes (a corner exactly on a node would give a deceptively exact answer, which is why the corner is placed at the irrational \(1/\pi\)). The convergence table shows orders that jump around and can even turn negative. The lesson: a high-order rule needs smoothness on the whole interval; put a node at every kink, i.e. integrate piece by piece.

The orders 2 and 4 are what you get for a generic smooth \(f\); special functions can do better. Take \(1/(1+x^2)\) on \([0,1]\), whose integral is \(\pi/4\). The trapezoid rule has order \(2.00\), as expected (\(f'(1)=-\tfrac12\ne0\)). But Simpson 1/3 and 3/8 both show order \(6.0\): errors \(2.2\cdot10^{-7}\), \(3.3\cdot10^{-9}\), \(5.2\cdot10^{-11}\) for \(n=6,12,24\).

The reason is in the Euler–Maclaurin series for Simpson, whose \(h^4\) coefficient is proportional to \(f'''(b)-f'''(a)\). Here \(f'''(x)=-24x(1-x^2)/(1+x^2)^4\) vanishes at both \(0\) and \(1\), so the \(h^4\) term is absent and \(h^6\) leads. The error bound with \(f^{(4)}\) is still true, just far from sharp. Practical consequence: estimate the order from the data (the convergence table), do not trust the textbook number blindly. As a bonus, \(4\,S_6\) with this integrand is \(3.14159178\), six correct digits of \(\pi\) from seven function values.

Why stop at \(n=3\) panels? Because higher closed rules inherit the weakness of equally spaced interpolation (the Runge phenomenon of the previous app). The weight table shows it. For \(n=6\) (Boole's neighbor, degree 7) all weights are positive. From \(n=8\) the weights change sign:

and for \(n=10\) the largest weight is \(7.14\) and the quantity \(\sum|c_j|/n\) has grown from \(1\) to \(3.06\) (it is \(1.45\) already at \(n=8\)). Since \(\sum c_j=n\) always, large positive and negative weights cancel, and any rounding error in \(f\) is amplified by that factor: the rule becomes unstable. The stable alternative is what we did: keep a low-degree rule (trapezoid, Simpson) and refine \(h\). The only way to get high degree without losing stability is to give up equally spaced nodes, which is the subject of the next class, Gaussian quadrature.

Class 1: from interpolation to the trapezoid rule.

• Lessons 1.1–1.2: interpolatory quadrature, weights \(w_j=\int\ell_j\), degree of precision, symmetry bonus (open the weights and exactness panels).
• Lessons 2.1–2.3: simple trapezoid, error proof with the mean value theorem for integrals, composite rule and order 2. In the app: sin, trapezoid, simple then composite, read the order in the convergence table.
• If time remains, 2.4 (Euler–Maclaurin and periodic functions) as a remark.

Class 2: Simpson and practice.

• Lessons 3.1–3.3: Simpson 1/3, the Hermite proof of \(-h^5f^{(4)}/90\), composite order 4. Lesson 2.5 (Richardson) connects it to the first class.
• Lessons 3.4–3.5: Simpson 3/8, comparison at equal cost, the odd-\(n\) combination.
• Lessons 4.1–4.4 as discussion: choosing \(n\), rough functions (\(\sqrt x\), the corner), the lucky \(1/(1+x^2)\), negative weights as the bridge to Gaussian quadrature.

1. Derive the weights \(\tfrac38,\tfrac98,\tfrac98,\tfrac38\) by solving the moment equations for \(1,x,x^2,x^3\) on \([0,3]\). Answer: the system is the one in the weights panel; the solution is \(c=(\tfrac38,\tfrac98,\tfrac98,\tfrac38)\), in units of \(h\), summing to \(3\).
2. For \(\int_0^1e^x\,dx\) with the composite trapezoid and \(n=10\), give the bound and compare with the app. Answer: bound \(e/1200=2.27\cdot10^{-3}\); actual \(-1.43\cdot10^{-3}\). With Simpson 1/3, \(n=10\): bound \(e/(180\cdot10^4)=1.51\cdot10^{-6}\); actual \(-9.5\cdot10^{-7}\).
3. Show that Simpson 1/3 integrates \(x^3\) exactly on any \([a,b]\). Answer: write \(x^3=(x-m)^3+3m(x-m)^2+3m^2(x-m)+m^3\); the odd powers of \(x-m\) cancel by symmetry and the rest has degree \(\le2\).
4. Estimate \(\pi=4\int_0^1\frac{dx}{1+x^2}\) with \(n=6\). Answer: trapezoid \(3.13696\) (error \(4.6\cdot10^{-3}\)); Simpson 1/3 \(3.14159178\) (error \(8.7\cdot10^{-7}\)).
5. Integrate \(e^x\) on \([0,1]\) with \(n=9\) using Simpson 1/3 on the first six subintervals and Simpson 3/8 on the last three, then the other way round. Answer: errors \(-2.27\cdot10^{-6}\) and \(-1.87\cdot10^{-6}\); the pure composite 3/8 with \(n=9\) gives \(-3.26\cdot10^{-6}\). Placing the (larger-error) 3/8 block where \(f^{(4)}\) is smaller, near \(x=0\), is better.
6. Explain why the trapezoid rule with \(n=8\) on \(e^{\cos x}\), \([0,2\pi]\) is accurate to about \(10^{-6}\) while on \(e^x\), \([0,1]\) it needs \(n\approx 380\) (trapezoid) for the same. Answer: lesson 2.4, periodicity.

A rule \(Q(f)=\sum_{j=1}^n w_jf(x_j)\) for \(\int f\,w\) has \(2n\) free numbers: \(n\) weights and \(n\) nodes. Newton–Cotes fixes the nodes and uses only the weights, so it reaches degree \(n-1\) (or \(n\)). If the nodes are free too, we can ask for exactness on \(1,x,\dots,x^{2n-1}\), which is \(2n\) equations:

n = 1 (Legendre, \([-1,1]\)): \(w_1=2\) and \(w_1x_1=0\), so \(x_1=0\): the midpoint rule, one evaluation, degree 1. n = 2: by symmetry \(w_1=w_2=1\), and \(2x^2=\mu_2=\tfrac23\) gives \(x=\pm\tfrac1{\sqrt3}=\pm0.5774\). Two evaluations and degree 3, as many as Simpson's rule uses with three. The system is nonlinear in the nodes, which is why the structure of the next lessons is needed.

No rule can do better than \(2n-1\). Take \(q(x)=\prod_j(x-x_j)^2\), of degree \(2n\): \(q\ge0\), so \(\int q\,w>0\), while \(Q(q)=0\).

Theorem. Let \(x_1,\dots,x_n\) be distinct, \(\omega(x)=\prod(x-x_j)\), and let the weights be those of interpolatory quadrature (\(w_j=\int\ell_jw\)). The rule has degree of precision \(\ge2n-1\) if and only if \(\int\omega\,q\,w=0\) for every polynomial \(q\) of degree \(\le n-1\).

Proof. (⇐) Let \(\deg f\le2n-1\). Divide: \(f=q\,\omega+r\) with \(\deg q,\deg r\le n-1\). Then \(\int fw=\int q\omega w+\int rw=\int rw\) by orthogonality. On the other side, \(f(x_j)=r(x_j)\) because \(\omega(x_j)=0\), so \(Q(f)=Q(r)=\int rw\), the last step because the rule is interpolatory and \(\deg r\le n-1\). (⇒) Take \(f=q\,\omega\) with \(\deg q\le n-1\), so \(\deg f\le2n-1\): exactness gives \(\int q\omega w=Q(q\omega)=0\) since \(\omega\) vanishes at the nodes. ∎

So the problem of choosing the nodes becomes: find the monic polynomial \(\omega\) of degree \(n\) orthogonal to all lower degrees, and take its roots. That polynomial is \(p_n\), the \(n\)-th orthogonal polynomial for the weight \(w\).

For a weight \(w>0\) on an interval (with finite moments) put \(\langle f,g\rangle=\int fg\,w\). Applying Gram–Schmidt to \(1,x,x^2,\dots\) gives a unique sequence of monic polynomials \(p_k\), \(\deg p_k=k\), with \(\langle p_j,p_k\rangle=0\) for \(j\ne k\). Equivalently, \(p_k\perp\) all polynomials of degree \(\lt k\) (they span the same space as \(p_0,\dots,p_{k-1}\)).

The Gram table in the app checks this numerically: for Legendre the diagonal reads \(2,\ 0.6667,\ 0.4,\ 0.2857,\ 0.2222\), i.e. \(2/(2k+1)\), and every off-diagonal entry is zero to rounding.

Orthogonal polynomials satisfy, for monic \(p_k\),

Proof. \(xp_k-p_{k+1}\) has degree \(\le k\), so \(xp_k=p_{k+1}+\sum_{j\le k}c_jp_j\) with \(c_j=\langle xp_k,p_j\rangle/h_j\). For \(j\lt k-1\): \(\langle xp_k,p_j\rangle=\langle p_k,xp_j\rangle=0\), because \(\deg(xp_j)=j+1\lt k\). Only \(j=k\) and \(j=k-1\) survive. ∎

The coefficients are explicit: \(\alpha_k=\langle xp_k,p_k\rangle/h_k\) and \(\beta_k=h_k/h_{k-1}\). In the usual normalizations the recurrences are, for Legendre \((k+1)P_{k+1}=(2k+1)xP_k-kP_{k-1}\); Chebyshev \(T_{k+1}=2xT_k-T_{k-1}\); Laguerre \((k+1)L_{k+1}=(2k+1-x)L_k-kL_{k-1}\); Hermite \(H_{k+1}=2xH_k-2kH_{k-1}\). The code in this page evaluates all four with one loop; the recurrence is also the best way to compute \(p_n(x)\) (it is cheap and stable).

Claim. \(p_n\) has \(n\) distinct real roots, all in the open interval \((a,b)\).

Proof. Let \(t_1,\dots,t_m\) be the points of \((a,b)\) where \(p_n\) changes sign (roots of odd multiplicity); \(m\le n\). Put \(q=\prod_{i\le m}(x-t_i)\). Then \(p_nq\) does not change sign on \((a,b)\), so \(\int p_nq\,w\ne0\). If \(m\lt n\), orthogonality of \(p_n\) to polynomials of degree \(\lt n\) would give \(0\): contradiction. Hence \(m=n\): all \(n\) roots are real, simple and in \((a,b)\). ∎

This is what makes the rule usable: the nodes are real, inside the domain of integration, never repeated. Look at the plot: for Legendre with \(n=4\) the roots are at \(\pm0.3400\) and \(\pm0.8611\), and the nodes crowd toward the ends of \([-1,1]\) as \(n\) grows (try \(n=30\)); for Laguerre and Hermite the roots spread over \([0,\approx4n]\) and \([-2\sqrt n,2\sqrt n]\).

Between two consecutive roots of \(p_n\) there is exactly one root of \(p_{n-1}\), and \(p_{n-1}\) has none outside \([x_1,x_n]\): the roots interlace. The plot shows it: the small dots (roots of \(p_{n-1}\)) sit in the gaps between the big dots (roots of \(p_n\)). It follows from the Christoffel–Darboux identity, or from the recurrence by induction: \(p_{n+1}(t_i)=-\beta_np_{n-1}(t_i)\) at a root \(t_i\) of \(p_n\), so \(p_{n+1}\) and \(p_{n-1}\) have opposite signs at the roots of \(p_n\), which forces a sign change of \(p_{n+1}\) in each gap.

Practical use: to find all roots of \(p_n\) one needs a bracket for each. The gaps between the roots of \(p_{n-1}\), plus the two end intervals, are \(n\) disjoint brackets each containing exactly one root. The code computes the rows of roots for \(p_1,p_2,\dots,p_n\) successively and in each gap runs Newton's method kept inside the bracket (bisection when a Newton step escapes). No initial guess is ever needed and no root can be missed or found twice, which is the usual failure of general polynomial root-finders at degree 20 or 30.

With the nodes fixed, the weights are those of interpolatory quadrature, \(w_j=\int\ell_j\,w\), where \(\ell_j\) is the Lagrange basis polynomial. Two facts:

Positivity. \(\ell_j^2\) has degree \(2n-2\le2n-1\), so the rule integrates it exactly: \(0\lt \int\ell_j^2w=\sum_kw_k\ell_j(x_k)^2=w_j\), since \(\ell_j(x_k)=\delta_{jk}\). Every weight is positive, and \(\sum w_j=\int w=\mu_0\).

Christoffel's formula. \(w_j=\dfrac{k_n\,h_{n-1}}{k_{n-1}\,p_n'(x_j)\,p_{n-1}(x_j)}\), with \(k_m\) the leading coefficient of \(p_m\). This is what the code uses, with \(p_n'\) obtained from the recurrence differentiated. Check in the app (Nodes and weights card): Legendre \(n=3\) gives weights \(\tfrac59,\tfrac89,\tfrac59\) at \(0,\pm0.7746\); Chebyshev gives \(\pi/n\) at every node.

Positive weights are the deep difference with high-order Newton–Cotes (whose weights, from \(n=8\), change sign): no cancellation, so the rule is numerically stable for every \(n\).

Let \(f\in C^{2n}\) and let \(q\) be the Hermite interpolant: degree \(\le2n-1\), \(q(x_j)=f(x_j)\), \(q'(x_j)=f'(x_j)\). The rule only sees the values, so \(Q(f)=Q(q)=\int qw\) (exact, degree \(2n-1\)). Thus \(E=\int(f-q)\,w\), and the Hermite error formula is

Now \(\omega^2\ge0\): constant sign, as in the Simpson proof. By the mean value theorem for integrals

Closed forms: Legendre \(\frac{(b-a)^{2n+1}(n!)^4}{(2n+1)((2n)!)^3}\); Chebyshev \(\frac{2\pi}{2^{2n}(2n)!}\); Laguerre \(\frac{(n!)^2}{(2n)!}\); Hermite \(\frac{n!\sqrt\pi}{2^n(2n)!}\). For \(n=2\), Legendre: \(E=f^{(4)}(\xi)/135\) on \([-1,1]\). Check with \(f=x^4\) (\(f^{(4)}=24\)): the rule gives \(2\cdot\tfrac19=\tfrac29\), the integral is \(\tfrac25\), the difference \(\tfrac8{45}=\tfrac{24}{135}\). ✓. With \(e^x\): the constant for \(n=2\) is \(1/135\), the bound \(e/135=2.0\cdot10^{-2}\), the actual error \(7.7\cdot10^{-3}\); for \(n=3\) the constant is \(6.35\cdot10^{-5}\) and the error \(6.5\cdot10^{-5}\); the app prints both, and the implied value \(f^{(6)}(\xi)=1.03\in[1,e]\).

Because the weights are positive and \(\sum w_j=\mu_0\), for continuous \(f\) on a finite interval and any polynomial \(r\) of degree \(\le2n-1\): \(|I-Q|\le|\int(f-r)w|+|Q(r-f)|\le2\mu_0\max|f-r|\). Choosing the best \(r\) (Weierstrass), Gauss quadrature converges for every continuous \(f\), and as fast as the best polynomial approximation of degree \(2n-1\): geometrically for analytic functions. Newton–Cotes cannot promise this (negative weights).

Numbers (Legendre, \(e^x\) on \([-1,1]\), same number of function evaluations): 3 points: Gauss \(6.5\cdot10^{-5}\), Simpson \(1.2\cdot10^{-2}\); 5 points: Gauss \(8.2\cdot10^{-10}\), Simpson \(7.9\cdot10^{-4}\); with \(n=6\) Gauss is already at \(1.6\cdot10^{-12}\), while Simpson with 11 points is at \(2.1\cdot10^{-5}\). Another instance from the Newton–Cotes app: \(\int_0^\pi\sin x\,dx\) with three evaluations: Simpson \(2.0944\) (error \(-0.094\)), Gauss–Legendre \(2.00139\) (error \(-1.4\cdot10^{-3}\)). Open the convergence card: the Gauss curve bends down faster than any straight line, while Simpson's decays like \(n^{-4}\).

Three ways to get the rule for a given \(n\). (1) Tables, fine for \(n\le5\) by hand. (2) The method of this page: roots of \(p_n\) by safeguarded Newton inside the interlacing brackets, weights by Christoffel's formula. (3) Golub–Welsch (1969): the recurrence says \(xp_k=p_{k+1}+\alpha_kp_k+\beta_kp_{k-1}\), i.e. the nodes are the eigenvalues of the symmetric tridiagonal Jacobi matrix with diagonal \(\alpha_k\) and off-diagonal \(\sqrt{\beta_k}\), and \(w_j=\mu_0\,v_{j,1}^2\) where \(v_{j,1}\) is the first component of the normalized eigenvector. This is the standard in libraries.

Double precision limits. The nodes and weights are accurate to near rounding for all \(n\le30\) here (Legendre \(n=20\): \(\sum w_j=2\) to \(4\cdot10^{-14}\)). The weights of Laguerre and Hermite decay super-exponentially: for \(n=10\) the smallest Laguerre weight is \(9.9\cdot10^{-13}\) at the node \(29.9\), and the smallest Hermite weight is \(7.6\cdot10^{-6}\). That is why the “degree of precision” table is shown only for \(n\le10\): the monomials \(x^k\), \(k\) near \(2n\), are a poor basis and the test itself, not the rule, loses digits.

Weight \(1\) on \([-1,1]\); roots of the Legendre polynomial \(P_n\). Explicit small cases:

For \([a,b]\): \(x=\frac{b-a}2t+\frac{a+b}2\), \(\int_a^bf=\frac{b-a}2\sum w_jf(x(t_j))\), and the error constant becomes \(\frac{(b-a)^{2n+1}(n!)^4}{(2n+1)((2n)!)^3}\) (for \(n=2\): \((b-a)^5/4320\)). This is the all-purpose rule for smooth integrands on a finite interval. In the app, Legendre, \(e^x\), \(n=1,2,3,4,5\): errors \(3.5\cdot10^{-1}\), \(7.7\cdot10^{-3}\), \(6.5\cdot10^{-5}\), \(2.9\cdot10^{-7}\), \(8.2\cdot10^{-10}\): each node multiplies the accuracy by a factor between \(45\) and \(360\). The test function \(x^6\) is exact for \(n\ge4\) (degree \(7\ge6\)) and has error \(4.6\cdot10^{-2}\) for \(n=3\), which equals \(720\cdot6.35\cdot10^{-5}\), exactly the error term since \(f^{(6)}=720\).

Weight \((1-x^2)^{-1/2}\) on \([-1,1]\). Here the roots of \(T_n\) are explicit, because \(T_n(\cos\theta)=\cos n\theta\):

Why equal weights: with \(x=\cos\theta\), \(\int_{-1}^1\frac{f(x)}{\sqrt{1-x^2}}dx=\int_0^\pi f(\cos\theta)\,d\theta\), and the rule is the midpoint rule in \(\theta\) for the even, periodic function \(f(\cos\theta)\), the case where the trapezoid/midpoint rules are spectrally accurate (Euler–Maclaurin lesson of the Newton–Cotes app). Use it for integrals with the factor \(1/\sqrt{1-x^2}\), which makes endpoint singularities harmless. Example in the app: \(\int_{-1}^1\frac{e^x}{\sqrt{1-x^2}}dx=\pi I_0(1)=3.97746326\); \(n=4\) gives \(3.97746264\), error \(6.3\cdot10^{-7}\), below the bound \(6.09\cdot10^{-7}\cdot e=1.7\cdot10^{-6}\); for \(n=2\): nodes \(\pm0.7071\), value \(\pi\cosh0.7071=3.9603\), error \(1.7\cdot10^{-2}\).

Weight \(e^{-x}\) on \([0,\infty)\): \(\int_0^\infty e^{-x}f(x)\,dx\approx\sum w_jf(x_j)\), \(E=\frac{(n!)^2}{(2n)!}f^{(2n)}(\xi)\). For \(n=2\): nodes \(2\mp\sqrt2=0.5858,\ 3.4142\), weights \(\frac{2\pm\sqrt2}4=0.8536,\ 0.1464\). The rule has degree 3: \(\int_0^\infty x^3e^{-x}=6\) is reproduced (\(5.99999\ldots\)) while for \(x^4\) it gives \(20\) instead of \(24\). The test function \(x^5\) is exact for \(n\ge3\).

Use it when the integrand decays like \(e^{-x}\) times something smooth. To integrate \(\int_0^\infty g\) take \(f=ge^{x}\), which is only sensible if \(g\) already decays at least like \(e^{-x}\). In the app: \(\int_0^\infty e^{-x}\sin x\,dx=\tfrac12\): the error oscillates in sign and falls to \(-1.9\cdot10^{-8}\) for \(n=12\) and \(1.5\cdot10^{-11}\) for \(n=16\). With a function that has a singularity at the endpoint, \(\sqrt x\) (exact \(\sqrt\pi/2\)) or \(\ln x\) (exact \(-\gamma\)), the convergence is only algebraic: for \(\sqrt x\) the error is \(1.3\cdot10^{-2},\ 4.6\cdot10^{-3},\ 1.6\cdot10^{-3}\) for \(n=4,8,16\), a factor \(2^{3/2}\) per doubling.

Weight \(e^{-x^2}\) on \((-\infty,\infty)\), \(E=\frac{n!\sqrt\pi}{2^n(2n)!}f^{(2n)}(\xi)\). Small cases:

The nodes are symmetric, so odd functions are integrated exactly for any \(n\). With \(n=2\): \(\int x^2e^{-x^2}=\tfrac{\sqrt\pi}2\) is exact (degree 3), \(x^4\) gives \(0.4431\) against \(\tfrac{3\sqrt\pi}4=1.3293\). Test \(x^6\) (exact \(\tfrac{15\sqrt\pi}8\)): error \(1.33\) for \(n=3\), zero for \(n\ge4\). Example with a transcendental function: \(\int e^{-x^2}\cos x\,dx=\sqrt\pi\,e^{-1/4}=1.38038845\); the errors for \(n=4,6,8\) are \(5.9\cdot10^{-5}\), \(3.7\cdot10^{-8}\), \(1.2\cdot10^{-11}\), and for \(n=12\) rounding level. Use Gauss–Hermite for expectations under a normal distribution (with \(x\mapsto\sqrt2\,x\) rescaling), integrals in quantum mechanics, and any integrand with a Gaussian factor. For \(1/(1+x^2)\) (exact \(1.3433\)) the convergence is slow and algebraic: poles at \(\pm i\) are too close to the real axis.

• Smoothness decides. For \(|x|\) on \([-1,1]\) (a corner at the middle) the errors for \(n=2,4,8,16,32\) are \(1.5\cdot10^{-1}\), \(4.3\cdot10^{-2}\), \(1.2\cdot10^{-2}\), \(3.0\cdot10^{-3}\), \(7.8\cdot10^{-4}\): order \(n^{-2}\) only, the same as the composite trapezoid rule. Split the interval at the kink and use Gauss on each piece. For \(1/(1+25x^2)\) the errors are \(3.4\cdot10^{-1}\), \(1.8\cdot10^{-1}\), \(4.1\cdot10^{-2}\), \(1.8\cdot10^{-3}\), \(3.1\cdot10^{-6}\): Gauss converges here (interpolation at equispaced nodes does not), but late, since the poles at \(\pm i/5\) are near the interval.
• Oscillation needs nodes. For \(\cos10x\) on \([-1,1]\): \(n=2,4\) are garbage (\(-1.85\), \(1.63\)), then \(5.3\cdot10^{-3}\) at \(n=8\), \(3.8\cdot10^{-5}\) at \(n=10\), \(1.6\cdot10^{-13}\) at \(n=16\). Rule of thumb: about \(\pi\) nodes per wavelength before the geometric regime starts.
• Composite Gauss. For a long interval or a rough function, cut into panels and use a fixed \(n\) (3 to 7) on each, like composite Simpson but with order \(2n\). Gauss nodes are not nested, so refining \(n\to2n\) reuses no evaluations; Gauss–Kronrod and Clenshaw–Curtis rules fix this and are what adaptive libraries use.
• Choose the family by the weight, not by the interval. If the integrand is \(e^{-x^2}g(x)\) with \(g\) smooth, Gauss–Hermite with \(f=g\) is exponentially better than truncating the integral and using Legendre.

Class 1: why the nodes are roots of orthogonal polynomials.

• 1.1–1.2: counting parameters, degree \(2n-1\), the orthogonality theorem (the proof is short and is the heart of the subject).
• 2.1–2.4: inner product, Gram–Schmidt, three-term recurrence, real simple roots, interlacing. In the app: change \(n\) and watch the roots and the small dots; read the Gram table.
• 3.1 (weights, positivity) and 4.1 (Gauss–Legendre, \(n=2,3\) by hand, then the app).

Class 2: error, the other three families, practice.

• 3.2 (Hermite-interpolation proof of the error term; reuse the Simpson proof from the previous class), 3.3 (convergence and Gauss against Simpson in the convergence card).
• 4.2–4.4: Chebyshev (explicit nodes), Laguerre, Hermite: for each, the weight, one explicit \(n=2\) computation, and one experiment in the app.
• 3.4 and 5.1 as discussion: Golub–Welsch, composite Gauss, rough functions.

1. Derive the Gauss–Legendre rule with \(n=2\) from the four moment equations. Answer: \(w_1=w_2=1\), \(x=\pm1/\sqrt3\) (lesson 1.1).
2. Show that it is not exact for \(x^4\) and compute the error two ways. Answer: \(\int x^4=\tfrac25\), \(Q=\tfrac29\), \(E=\tfrac8{45}\); the formula gives \(24/135=\tfrac8{45}\).
3. Estimate \(\int_0^\pi\sin x\,dx\) with Gauss–Legendre \(n=2\) and \(n=3\) (change variables to \([-1,1]\)). Answer: \(1.93582\) (error \(6.4\cdot10^{-2}\)) and \(2.00139\) (error \(-1.4\cdot10^{-3}\)).
4. With the Chebyshev rule, \(n=2\), estimate \(\int_{-1}^1\frac{e^x}{\sqrt{1-x^2}}dx\). Answer: \(\pi\cosh(1/\sqrt2)=3.9603\); the exact value is \(3.97746\), error \(1.7\cdot10^{-2}\).
5. Check that the Laguerre rule with \(n=2\) is exact for \(x^3\) but not for \(x^4\). Answer: \(0.8536\cdot0.5858^3+0.1464\cdot3.4142^3=6\); for \(x^4\) it gives \(20\neq24\).
6. Show that Gauss–Hermite with \(n=2\) integrates \(x^2e^{-x^2}\) exactly. Answer: \(2\cdot\tfrac{\sqrt\pi}2\cdot\tfrac12=\tfrac{\sqrt\pi}2\).
7. Prove that no \(n\)-node quadrature rule is exact for all polynomials of degree \(2n\). Answer: \(q=\prod(x-x_j)^2\) (lesson 1.1).

The Tacoma Narrows Bridge in Washington State had a main span of 2,800 ft (\(853\,\mathrm{m}\)) and a deck only about 39 ft (\(12\,\mathrm{m}\)) wide, slender for its length. It opened on 1 July 1940 and from the first days it bounced in moderate wind. On 7 November 1940 the bouncing turned into a large torsional motion, the two edges of the deck going up and down in opposition, and within about an hour or so the main span fell. Nobody died; a dog in a car was the only casualty. Ambient wind was about 40 mph (64 km/h). Sources differ on some figures (the day is sometimes given as 10 November in older accounts), so treat the numbers as approximate.

Not in the textbook sense. The popular story is that the wind pushed at the bridge's natural frequency, like a child on a swing. Billah and Scanlan (1991) showed that physics textbooks repeated this and that it is not what the evidence supports. The vortices shed by the deck had a frequency far from the torsional frequency that appeared. The explanation most engineers accept is aeroelastic flutter: when the deck twists, the airflow around it changes in a way that pushes it further, a self-excited feedback that behaves like negative damping above a critical wind speed. A resonance needs an external rhythm; flutter creates its own.

Other mechanisms have been proposed, and published surveys state that no explanation of every detail is universally accepted. This app is about one of them, not about the verdict.

Because it is a clean example of a nonlinear mechanism that couples two motions that the linear theory keeps apart. Lazer and McKenna (1990) pointed out that a cable stretches but does not compress, so a hanger is a spring that turns off when slack. McKenna and Tuama (2001) showed in a simple model that vertical forcing can then create a torsional response. Critics, for example Green and Unruh, argue that the Tacoma hangers did not in fact go slack and that the parameters needed are unrealistic. So treat the card as a laboratory for the mechanism and for the numerics, not as a reconstruction of 7 November 1940.

The deck section is a rigid bar of half-width \(\ell\), mass \(m\) per metre. Its centre moves down by \(y\) and it rotates by \(\theta\), so the left and right edges move down by \(y\mp\ell\sin\theta\). Each hanger is stretched by \(f=y_0+y\mp\ell\sin\theta\), where \(y_0=mg/(2K)\) is the static stretch, and pulls with \(K f^{+}\) with \(f^+=\max(f,0)\). Newton for the centre gives \(m\ddot y=-K(f_L^++f_R^+-2y_0)\) plus damping and forcing. For the rotation, the torque of a vertical force \(F\) at distance \(\ell\cos\theta\) from the centre is \(F\ell\cos\theta\), and the moment of inertia of a bar about its centre is \(\int_{-\ell}^{\ell}(m/2\ell)\,x^2\,dx=m\ell^2/3\). That gives the second equation in the summary card.

Then \(f^+=f\), and for small \(\theta\) the equations become \(\ddot y+\delta\dot y+\tfrac{2K}{m}y=\varepsilon\sin\mu t\) and \(\ddot\theta+\delta\dot\theta+\tfrac{6K}{m}\theta=0\): two independent oscillators. Vertical forcing drives only \(y\); the twist \(\theta\) has no forcing and any initial twist simply dies out. With \(K=8000\,\mathrm{N/m^2}\) and \(m=2500\,\mathrm{kg/m}\), \(\omega_y=\sqrt{2K/m}=2.53\) rad/s (period \(2.48\,\mathrm{s}\)) and \(\omega_\theta=\sqrt3\,\omega_y=4.38\) rad/s (period \(1.43\,\mathrm{s}\)). Run “Same gust, hangers kept taut”: the twist stays at about 0.06°, the seed we gave it. The ratio \(\sqrt3\) is not a coincidence of the numbers: it comes only from the geometry, inertia \(m\ell^2/3\) against two springs at distance \(\ell\).

At rest the two hangers carry the weight: \(2Ky_0=mg\), so \(y_0=mg/(2K)=2500\cdot9.8/16000=1.53\) m. This number matters for slackness: a hanger goes slack when its stretch \(y_0+y\mp\ell\sin\theta\) becomes negative, that is, when the deck center rises more than \(y_0\) above its rest position (for \(\theta=0\)), or earlier when the deck is tilted. A stiffer hanger (larger \(K\)) gives a smaller \(y_0\), so slackness is reached for a smaller motion, but \(\omega_y\) is also higher. Earlier trials with \(K=1000\) gave an unphysical \(12\,\mathrm{m}\) sag and twists of tens of radians, so the default is \(K=8000\).

In the taut model the torque is \(-2K\ell^2\sin\theta\cos\theta\) and does not depend on \(y\). With slackness it becomes \(K\ell\cos\theta\,(f_L^+-f_R^+)\). If the deck is lifted high enough, both hangers are slack and the torque vanishes, so the twist is freed from its restoring spring. When one side goes slack before the other, which happens as soon as \(\theta\neq0\), the torque is not symmetric and depends on the heave. In this way the heave motion changes the stiffness of the torsion mode over each cycle, and a stiffness that varies with the right rhythm can pump energy into an oscillation (parametric excitation, as when a child pumps a swing by standing and squatting). That is a plausible reading, not a proof: the experiments show what happens, the mathematics of when is hard.

Run “Gust at the heave frequency, slack hangers” and watch the edge-lift curve: a tiny twist of 0.06° is amplified to the full 5.7° threshold and beyond.

With RK4 and \(h=0.02\) s the scan gives, for example, \(\varepsilon_{crit}\approx6.6\,\mathrm{m/s^2}\) at \(\mu/\omega_y=0.6\), 3.2 at 0.8, 0.78 at 1.0, 1.4 at 1.04, 7.4 at 1.08, and no twisting up to \(8\,\mathrm{m/s^2}\) from 1.16 on. Three lessons. First, there is a valley near the heave frequency, as resonance would suggest. Second, it is lopsided: weaker forcing is enough on the low-frequency side than on the high side, which fits the picture of a spring that softens when hangers go slack (a softening spring resonates at lower frequency when driven hard); this is a plausible reading, not something the app proves. Third, the curve is jagged: with a nonlinear coupling the response is not monotone in \(\varepsilon\), so “the smallest \(\varepsilon\) that twists” depends on how you search. We scan upward on a geometric grid, bracket the first crossing and bisect inside it; a threshold hidden between two scan points could be missed.

Euler's method turns \(y''=-\omega^2y\) into a map whose amplitude grows each step by \(\sqrt{1+z^2}\), \(z=\omega h\): an undamped oscillation spirals outward whatever \(h\) you choose. For the torsion mode at \(h=0.02\) s, \(z=0.088\) and the amplitude grows by 0.4% per step, about 21% per simulated second, and the energy by about 47% per second at the start. Try “Free fall, no wind: Euler”: the energy curve climbs and the status line warns that energy grew with no wind. The same trajectory with RK4 and the same \(h\) stays within \(3\cdot 10^{-4}\) of the initial energy over \(30\,\mathrm{s}\), and Heun within \(5\cdot 10^{-3}\) (maximum relative drift over \(30\,\mathrm{s}\): Euler 4.9, Heun \(4.7\cdot 10^{-3}\), RK4 \(2.9\cdot 10^{-4}\) at \(h=0.02\); at \(h=0.005\): 1.1, \(1.3\cdot 10^{-4}\), \(1.8\cdot 10^{-5}\).

The consequence for the critical wind: with Euler the scan returns thresholds near zero everywhere, because the method itself drives the twist.

Euler has local error \(O(h^2)\) and global error \(O(h)\); Heun (the explicit trapezoid rule) global \(O(h^2)\); the classical Runge–Kutta method global \(O(h^4)\). In the step-size card, on a smooth trajectory (hangers always taut), the measured slopes are 0.95, 2.02 and 4.02. When the hangers switch on and off the right-hand side has corners, so the Taylor expansion behind the order proofs fails near each switch: Euler and Heun keep roughly their order (0.99 and 2.04), but RK4 falls to about 2.4 and its error curve is irregular, depending on where the corners fall between steps. Moral: a high-order method does not buy you high accuracy on a non-smooth problem; to recover it you would have to locate each switching time and restart the integrator there (event handling).

For \(y''=-\omega^2y\) one RK4 step multiplies the state by \(R(z)=1-\tfrac{z^2}{2}+\tfrac{z^4}{24}+i\bigl(z-\tfrac{z^3}{6}\bigr)\), \(z=\omega h\), and \(|R|\le1\) only for \(z\le2\sqrt2\approx2.83\). With \(h=0.7\) s the torsion mode has \(z=4.38\cdot0.7=3.07\), just past the limit, and \(|R|=1.75\): the oscillation grows by 75% per step. Run “RK4 with a huge step”: the energy rises by a factor of several hundred although damping is on and no wind blows. In the nonlinear model the twist does not grow without bound (the \(\sin\theta\cos\theta\) torque saturates), so you see an erratic motion at large amplitude instead of an overflow: the energy curve is the reliable alarm. Below the limit RK4 is stable and slightly damps high frequencies (at \(z=2\) the amplitude loses about 25% per step), so too large a step can also make a real oscillation die from numerical damping alone.

You never know for sure; you collect evidence. (1) Check a case with a known answer: the taut model has periods \(2\pi/\omega_y=2.484\) s and \(2\pi/\omega_\theta=1.434\) s, and the simulation reproduces them to four digits. (2) Check an invariant: with no wind and no damping the energy must be flat. (3) Halve the step and see whether the answer changes by about the expected \(2^{-p}\); for RK4 on smooth parts, a factor 16. (4) Compare two methods: Heun and RK4 give the same critical thresholds to three digits (0.776 against 0.78 at \(\mu/\omega_y=1\)) while Euler does not, which tells you which to distrust. Agreement is not proof, but disagreement is a certain warning.

For teaching: with a fixed step, the formulas for the amplification factor and the order apply directly, and you can see how the error depends on \(h\). A production code would use an adaptive embedded pair, such as Dormand–Prince, to choose \(h\) from an error estimate, and an event finder to land exactly on each slack/taut switch. Adaptive steps shrink automatically near the switches, which is where the fixed-step RK4 loses its order.

Start with “Gust at the heave frequency”. Keep the strength at 1.2 and move the frequency slider from 0.90 to 1.10 in steps of 0.02. In our runs the deck twists only for \(0.94\lesssim\mu/\omega_y\lesssim1.04\), reaching 12° to 22° within two minutes, and it takes about \(25\,\mathrm{s}\) to cross the 5.7° threshold at 1.00. Outside that band the twist stays at its seed. Then lower the strength to 0.5: no frequency twists the deck at all, and the lowest threshold of the whole scan is \(0.78\,\mathrm{m/s^2}\). Compare with the critical-wind curve (run the scan first): the band where it dips below 1.2 is the band you found.

Change \(K\) from 8000 to 2000 and then to 30000 with the gust frequency fixed at 1.0 times \(\omega_y\). Predict first: how do \(\omega_y\), \(\omega_\theta\) and \(y_0\) scale with \(K\)? (Answers: \(\omega\propto\sqrt K\), halving to 1.26 and \(2.53\,\mathrm{rad/s}\) when \(K\) is divided by 4, and \(y_0\propto1/K\), growing to \(6.1\,\mathrm{m}\).) Then check the readout beside the slider. Which setting twists more easily? Notice that the preset menu goes to “Custom” and the simulation restarts, because changing \(K\) changes what “rest” means.

Write \(y''=-\omega^2y\) as a first-order system \(u'=Au\), \(u=(y,y')\), and apply one step of each method. Since \(A\) has eigenvalues \(\pm i\omega\), it is enough to apply the method to the scalar test equation \(u'=\lambda u\), \(\lambda=i\omega\). Euler: \(R=1+\lambda h=1+iz\), so \(|R|=\sqrt{1+z^2}\). Heun: \(R=1+\lambda h+\tfrac12(\lambda h)^2=1-\tfrac{z^2}2+iz\), and \(|R|^2=(1-\tfrac{z^2}{2})^2+z^2=1+\tfrac{z^4}{4}\). RK4: \(R\) is the Taylor polynomial of \(e^{\lambda h}\) up to order 4: \(1+\lambda h+\tfrac{(\lambda h)^2}{2}+\tfrac{(\lambda h)^3}{6}+\tfrac{(\lambda h)^4}{24}\). Show that \(|R|=1\) is reached at \(z=2\sqrt2\) (there \(R=-\tfrac13-\tfrac{2\sqrt2}{3}i\)) and that near \(z=0\), \(|R|^2=1-z^6/72+\ldots\) so RK4 damps slightly. Check against the plot.

List what the model leaves out and say which would matter most. A partial list: the wind force depends on the angle of attack and therefore on \(\theta\) and \(\dot\theta\) (this is flutter); the cables and towers have their own motion; the bridge is a long structure with many modes, not one section; the hangers have different lengths and tensions; the damping is not a single viscous constant. Then decide which new term would turn the damping \(\delta\) negative above a critical wind speed, which is the signature of flutter. (Hint: a force proportional to \(+\dot\theta\) in the torque equation.)

Walking sideways pushes on the deck with a force that alternates left and right once per two steps, at roughly \(1\,\mathrm{Hz}\), near the lateral frequencies of many footbridges. If the pushes of different people are random they largely cancel, and the bridge only feels a force of order \(\sqrt N\). If the people are in step with each other and with the sway, the forces add up to order \(N\) and push the bridge in phase with its own velocity, which is exactly what negative damping means.

It is a feedback loop. More sway makes more walkers adjust their pace to it, more walkers in step make more sway. Whether the loop amplifies or dies depends on whether the push from the crowd beats the damping of the bridge. Below the threshold any small sway is damped out and the walkers stay at their own paces; above it the loop runs away until it saturates, because walkers whose natural pace is far from the bridge's cannot lock. Try the “below critical” preset (a push dies out) and “above critical” (a tiny seed grows over a couple of minutes).

The critical crowd is proportional to the damping coefficient: \(N_c=4\zeta K\gamma/(GC\sin\alpha)\propto B\). Doubling the damping doubles the threshold. Stiffness enters only through \(K\), so tripling the threshold by stiffening would need three times the stiffness, moving the frequency by a factor \(\sqrt3\), and a stiffer bridge is still excited by a crowd at its own frequency. Adding damping raises the threshold directly and cheaply. In the app the “Fit dampers” button raises \(\zeta\) to 3 %, so that \(N_c\) goes from 143 to 428 and a wobbling crowd of 250 has its sway brought below 1 cm in about \(35\,\mathrm{s}\).

One lateral mode of the bridge is an oscillator \(M\ddot X+B\dot X+KX=G\sum\sin\Theta_i\), where \(\Theta_i\) is the phase of the \(i\)-th walker's sway cycle and \(G\) the lateral force amplitude of one walker. Each walker obeys \(\dot\Theta_i=\omega_i+CA\sin(\Psi-\Theta_i+\alpha)\): its own pace \(\omega_i\) plus a pull toward the bridge phase \(\Psi\), proportional to the bridge amplitude \(A\). With \(X=A\sin\Psi\), \(\dot X=\Omega_0A\cos\Psi\) the amplitude and phase are computed from \(X\) and \(\dot X\) at each step. With \(\alpha=\pi/2\) the walkers' force runs a quarter turn ahead of the displacement and so pushes in phase with the velocity. The \(N+2\) equations are integrated together. This is a Kuramoto model in which the coupling is not fixed but is carried by the bridge itself.

Look for a steady state with sway amplitude \(A\). In the frame turning with the bridge, a walker with detuning \(\nu=\omega_i-\Omega_0\) obeys \(\dot\varphi=\nu-CA\sin\varphi\). Those with \(|\nu|\lt CA\) lock at \(\sin\varphi=\nu/CA\); the others drift and on average contribute nothing to the force. The locked walkers give an order parameter \(r=\int_{-CA}^{CA}P(\nu)\sqrt{1-\nu^2/(CA)^2}\,d\nu\), and the bridge, driven at resonance, has its damping balanced by the part of the crowd force in phase with its velocity: \(B\Omega_0A=GNr\sin\alpha\). For small \(A\), \(r\approx\tfrac\pi2P(0)\,CA\); inserting it gives a nonzero \(A\) only if \(N>N_c=\dfrac{2B\Omega_0}{\pi\,G\,C\,P(0)\sin\alpha}\). With a Lorentzian, \(P(0)=1/(\pi\gamma)\) and \(N_c=2\gamma B\Omega_0/(GC\sin\alpha)\). Equivalently \(N_c=4\zeta K\gamma/(GC\sin\alpha)\), proportional to the damping ratio \(\zeta\).

For a Lorentzian the integral can be done in closed form: \(r=(\sqrt{a^2+\gamma^2}-\gamma)/a\) with \(a=CA\). The balance \(B\Omega_0A=GNr\sin\alpha\) then becomes \(\sqrt{a^2+\gamma^2}-\gamma=q\,a^2\), whose solution is \(a^2=(1-2\gamma q)/q^2\); hence \(A_\ast=\dfrac{G\sin\alpha}{B\Omega_0}\sqrt{N(N-N_c)}\) and \(r_\ast=\sqrt{1-N_c/N}\). A true Lorentzian has tails reaching negative walking frequencies, so the app cuts it at \(\pm T=\pm\min(5\gamma,0.8\Omega_0)\) and renormalizes; the kept share is \(Z=\tfrac2\pi\arctan(T/\gamma)=0.874\) for the default settings. Everything above then holds with \(N_c\) multiplied by \(Z\), \(A_\ast\) and \(r_\ast\) divided by \(Z\), as long as \(CA_\ast\le T\).

With the default settings the predicted critical crowd is \(N_c=142.6\). Simulating \(1500\,\mathrm{s}\) and averaging the last \(150\,\mathrm{s}\) gives, for \(N=150,175,200,250,300,400\), sways of 2.15, 5.15, 7.40, 11.36, 15.09 and 22.31 cm against the predicted 2.33, 5.28, 7.51, 11.48, 15.23 and 22.49 cm, within 1 to 8 %; the discrepancy is largest closest to threshold, where growth is slow. The order parameter agrees within 0.005 for \(N\ge200\). Fitting \(A^2/N\) against \(N\) to the simulated sways for \(N=175\) to \(300\) and reading off the intercept gives \(N_c=143.7\), within 1 % of the formula. Walkers are sampled at quantiles, so there is no randomness: the runs are exactly reproducible.

The state has \(N+2\) numbers: bridge displacement and velocity and \(N\) phases. The fast part is the oscillation of the bridge at \(\Omega_0\) and of the walkers' phases at their own paces, all within a few rad/s, and the slow part is the growth or decay of the amplitude over minutes. A fourth-order method lets one step cover a good fraction of a period, so that 25 minutes of simulated time cost under 20 000 steps. The cost per step is proportional to \(N\) (four evaluations, each one sine per walker), so one run with 250 walkers takes under a second.

Two limits. For the bridge alone, the amplification factor per step of each method is a function of \(z=\Omega_0h\): \(|R|=\sqrt{1+z^2}\) for Euler, \(\sqrt{1+z^4/4}\) for Heun, both larger than 1 for every \(z>0\), so they pump energy at the rate \(\Omega_0^2h/2\) (Euler) and \(\Omega_0^4h^3/8\) (Heun) per second. The true damping is \(\zeta\Omega_0\), so a method fakes a wobble when its pumping beats it: \(h>2\zeta/\Omega_0=0.004\) s for Euler and \(h>(8\zeta/\Omega_0^3)^{1/3}=0.086\) s for Heun. These agree with the table above (Euler grows already at \(0.004\,\mathrm{s}\), Heun at \(0.09\,\mathrm{s}\)). For the coupled system with RK4, the error in the final sway relative to \(h=0.01\) is \(-76\%\) at \(h=0.2\), \(-3.4\%\) at 0.1, \(-0.12\%\) at 0.05 and below 0.005 % at 0.02.

Random frequencies would make each run different, and a threshold measured on a few hundred walkers would carry noise of several percent. Placing walker \(i\) at the quantile \((i+\tfrac12)/N\) of the distribution makes the sum over walkers a midpoint rule for the integral that defines \(r\), so the run is deterministic and the finite-\(N\) error is smooth. The initial phases are a golden-ratio sequence \(2\pi\{i\varphi\}\), which spreads evenly (so the initial order parameter is small), and both sequences are shuffled so that the first \(n\) walkers, the ones who have entered when the crowd arrives gradually, still represent the whole distribution.

Pick the preset “below critical”, press “Push the bridge” a few times, and watch the sway die. Raise \(N\) by 10 each time and repeat. At which \(N\) does a 1 cm push stop dying and start to grow? Because growth near threshold is slow, wait a few minutes of simulated time (speed \(32\times\) and compare with the predicted \(N_c\) in the metrics row. Then double \(\zeta\) and see how the threshold moves.

Run “Walkers with nearly equal paces”. The critical crowd falls from 143 to 54, because \(N_c\) is proportional to the spread \(\gamma\): when everyone naturally walks at nearly the same pace, more of the crowd is within reach of the bridge frequency and a smaller crowd suffices. With the settings of that preset and \(N=100\), the predicted sway is 4.8 cm and the simulated one 4.7 cm. Now lower \(C\) to 8 and check that \(N_c\) doubles: it is inversely proportional to the sensitivity.

Run “No walkers, explicit Euler”. The bridge is empty and damped, yet the amplitude grows. Lower \(h\) to \(0.004\,\mathrm{s}\) and see the growth stop; lower it again and see the decay that should be there. Switch to Heun with \(h=0.1\): also growth. Then switch to RK4. What would a model of the Tacoma bridge, a real bridge whose behaviour depended on a small damping, have said if its author had used explicit Euler?

In general relativity gravity is the curvature of spacetime, and light, which has no mass to push around, simply follows the straightest possible paths in that curved spacetime, the null geodesics. Near a mass these paths are no longer straight lines in the space we draw. A Newtonian argument that treats light as particles of speed \(c\) also predicts a bend, but half as large: \(2M/b\) against Einstein's \(4M/b\). The other half comes from the curvature of space itself, not only from the time part of the metric.

The horizon \(r=2M\) is the point of no return: nothing inside can get out, light included. At \(r=3M\), the photon sphere, light can travel in a circle around the hole, but the orbit is unstable like a ball balanced on a hilltop: the slightest push sends the light in or out. A ray arriving from far away is captured if its impact parameter is below \(b_c=3\sqrt3\,M\approx5.196\,M\), which is not \(3M\) because of the bending on the way in. The dark disc that a distant observer sees, the shadow, has radius \(b_c\).

Imagine the ray coming from very far away along a straight line. The impact parameter \(b\) is the distance between that line and the parallel line through the hole: the height of the ray in the top view. It is the ratio \(L/E\) of two conserved quantities, angular momentum and energy, so it is the same at every point of the ray and it alone decides the whole orbit. That is why the app needs only one number to launch a ray: \(b=8\) always gives the same curve.

Write the conserved quantities \(E=(1-2M/r)\dot t\) and \(L=r^2\dot\varphi\) (dots are derivatives in the affine parameter) and the null condition \(0=-(1-2M/r)\dot t^2+\dot r^2/(1-2M/r)+r^2\dot\varphi^2\). Eliminating \(t\) gives \(\dot r^2=E^2-(L^2/r^2)(1-2M/r)\). With \(u=1/r\) and \(du/d\varphi=-\dot r/L\) this becomes \((du/d\varphi)^2=1/b^2-u^2+2Mu^3\), the first integral. Differentiating once in \(\varphi\) and dividing by \(2u'\) gives \(u''+u=3Mu^2\).

Only the last two lines are needed in this course; the derivation from the metric is a good exercise for a relativity class.

Put \(y_1=u\) and \(y_2=u'\). Then \(y_1'=y_2\) and \(y_2'=3y_1^2-y_1\) (with \(M=1\)): a first-order system \(\mathbf y'=\mathbf f(\mathbf y)\) in which the angle \(\psi\) does not appear on the right (an autonomous system). Euler: \(\mathbf y_{n+1}=\mathbf y_n+h\,\mathbf f(\mathbf y_n)\). Heun: average the slope at the start and at the Euler prediction. RK4: four slopes with weights \(1,2,2,1\) over 6. The file bh-core.js has the three in about twenty lines, and so has the Python script at the end of the page.

The horizon is at \(u=1/2\). After each step the program checks whether \(u\ge1/2\); if so the ray is captured, the bell rings, and the last step is shortened so that the path ends exactly on the horizon. The ray escapes when \(u'<0\) and \(u\) is back at its starting value \(u_0=1/r_0\). If neither happens after \(14\pi\) radians the ray is declared undecided, still circling. Locating the last step is a root-finding problem: a few Newton iterations on the fraction of the step, done with the same method, so the method keeps its order.

With no gravity the ray would sweep the polar angle from \(\psi_0=\arcsin(b/r_0)\) to \(\pi-\psi_0\). The bent ray sweeps a larger angle \(\psi_{\text{exit}}-\psi_0\), and the difference is the deflection \(\alpha=(\psi_{\text{exit}}-\psi_0)-\pi+2\psi_0\). The program starts at \(r_0=1000M\), far enough that the result is close to the value for a ray coming from infinity. To first order \(\alpha=4M/b\); the expansion continues as \(4M/b+15\pi M^2/(4b^2)+128M^3/(3b^3)+\dots\), and it fails near \(b_c\), where \(\alpha\) diverges.

Linearize the orbit equation about the circular photon orbit \(u=1/3\): with \(u=1/3+\delta\) we get \(\delta''=\delta+3\delta^2\), so small errors grow like \(e^{\psi}\) (a saddle, not a centre). A ray whose \(b-b_c\) is \(10^{-6}\) spends about \(\ln10^{6}\approx14\) radians, more than two turns, near the photon sphere before it decides. A numerical error of that size is amplified in the same way, so locating \(b_c\) to many digits needs a method whose error is far below the digits wanted. RK4 with \(h=0.02\) has it; Euler does not.

The equation is a nonlinear oscillator, \(u''=-u+3u^2\). Euler's method applied to an oscillator multiplies the amplitude by \(\sqrt{1+h^2}\) at every step: the numerical solution gains energy. Here the energy is the first integral \(1/b^2\), and gaining it means behaving like a ray with a smaller \(b\): Euler pulls in rays that should pass. Heun's amplification factor is \(\sqrt{1+h^4/4}\), smaller by two powers of \(h\), and RK4's is about \(1-h^6/144\) (it loses a little), so these two are far better. The same stability study is done in the Tacoma app.

In Schwarzschild coordinate time a ray approaching the horizon slows down and never gets there: \(dr/dt\to0\). A method marching in \(t\) would crawl and never see the capture. In the angle \(\psi\) the equation \(u''=3u^2-u\) is smooth everywhere and the horizon is simply the value \(u=1/2\). For the movie the app uses the affine parameter, \(d\psi/d\lambda=bu^2\), which far away is almost the ordinary time and keeps the motion lively near the hole.

When \(b\to0\) the ray goes almost straight in and the angle \(\psi\) hardly changes, so a fixed step in \(\psi\) covers nearly the whole fall: the slider stops at \(|b|=1\) for that reason. For radial rays one uses another formulation, for instance \(r(\lambda)\) with \((dr/d\lambda)^2=1-b^2(1-2M/r)/r^2\), or an adaptive step. Choosing the independent variable well is part of the modelling.

Choose “Euler with a big step lies” (\(b=5.3\), \(h=0.1\)). The exact answer is that the ray escapes, bent by 204°, but Euler makes it fall in, because its first-integral error moves the effective critical value up to \(b_c\approx5.356\). Press “Run the bisection” to see that number. Then lower \(h\) with the slider: below about \(h=0.07\) the ray escapes again, and the estimate of \(b_c\) tends to \(5.196\) from above, roughly in proportion to \(h\). Do the same with Heun and with RK4 at \(h=0.1\): both get the fate right.

Cast \(b=6\), \(5.3\) and \(5.2\) in turn with the experiment menu and read the bending in the status line: 99°, 204° and 390°. At \(b=5.2\), only \(0.004\) above \(b_c\), the ray goes more than once around the hole before leaving. By the linearization of the previous section each extra full turn needs the ray to be \(e^{2\pi}\approx535\) times closer to \(b_c\), so the slider's step of \(0.05\) is far too coarse for a second loop; use the bisection of the Python script. Does the logarithmic law \(\alpha\approx-\ln(b-b_c)+\text{const}\) fit the green curve near \(b_c\)?

1. Show from the first integral that capture happens exactly for \(b<3\sqrt3\,M\): the ray cannot turn back when \(F(u)=1/b^2-u^2+2u^3\) stays positive on \([0,1/2]\); \(F\) has its minimum there at \(u=1/3\), where \(F=1/b^2-1/27\).
2. Write Euler, Heun and RK4 for \(u''+u=3u^2\) in your own language and reproduce the table of the capture boundary.
3. Use \(\alpha\approx4M/b\) to estimate the bending of light at the rim of the Sun (\(M=1.48\,\mathrm{km}\), \(b=6.96\times10^{5}\,\mathrm{km}\)) in arcseconds, and compare with \(1.75''\).
4. Take \(r\) as the independent variable and the affine parameter as time, and try to integrate a ray with \(b\) just above \(b_c\). What happens at the turning point, where \(dr/d\lambda=0\)?
5. With the linearization about \(u=1/3\), estimate how many turns around the hole a ray with \(b-b_c=10^{-8}M\) makes.

The equation \(y^2=k\) is easy: \(y=\pm\sqrt k\). The trouble with \(x^2+px=q\) is that \(x\) appears twice, once squared and once alone. Completing the square rewrites the left side as one square, \((x+p/2)^2\), so that \(x\) appears once and the square root can peel it off. The geometric picture shows why it works: a square plus a strip of the same height is almost a bigger square, and the missing corner has area \((p/2)^2\).

The first real step divides by \(a\). If \(a=0\) the equation \(bx+c=0\) is linear, with the single solution \(x=-c/b\), and the formula would divide by zero. There is a pleasant limit: as \(a\to0\) with \(b>0\), the root \(\frac{-b+\sqrt{b^2-4ac}}{2a}\) tends to \(-c/b\) while the other root runs off to infinity. The stable form \(x=c/q\) shows this limit at once, which the naive form hides behind \(0/0\). Try the preset with a tiny \(a\).

Both signs give a solution: if \(y^2=k\) then \((-y)^2=k\) as well. So the formula is a short way to write two numbers, \(x_1=\frac{-b+\sqrt\Delta}{2a}\) and \(x_2=\frac{-b-\sqrt\Delta}{2a}\), which are mirror images of each other around the vertex \(x=-b/2a\). If \(\Delta=0\) they coincide (a double root, the parabola just touches the axis); if \(\Delta<0\) they are complex conjugates and the parabola does not reach the axis. In the animation the sign appears at the square-root step and survives to the end.

Dividing by \(a\ne0\), adding the same number to both sides and rewriting an expression as an equal one can each be undone, so every equation in the chain has exactly the same solutions as the previous one. Only the square-root step needs care: from \(y^2=k\) to \(y=\pm\sqrt k\) it keeps both solutions. A concrete check: for any \(x\), the difference between the two sides of each of the first seven frames equals \(f(x)/a\) with \(f(x)=ax^2+bx+c\) (or \(f(x)\) itself in the first frame), so all seven vanish at exactly the same \(x\). The Python script below verifies this with exact fractions.

\(\Delta=b^2-4ac\) is the number under the root, and the root is real exactly when \(\Delta\ge0\). Geometrically, the parabola \(y=ax^2+bx+c\) has its vertex at \(x=-b/2a\) with height \(-\Delta/4a\), so it crosses the axis twice when \(a\) and \(-\Delta\) have the same sign (\(\Delta>0\)), touches it once when \(\Delta=0\), and misses it when \(\Delta<0\). The distance between the two roots is \(\sqrt\Delta/|a|\). For integer coefficients \(\Delta\) is an integer, and the roots are rational exactly when \(\Delta\) is a perfect square.

Add and multiply the two roots, with \(s=\sqrt\Delta\): \(x_1+x_2=\frac{-b+s}{2a}+\frac{-b-s}{2a}=-\frac ba\) and \(x_1x_2=\frac{(-b)^2-s^2}{4a^2}=\frac{b^2-(b^2-4ac)}{4a^2}=\frac ca\). These are Viète's relations, and they say \(a(x-x_1)(x-x_2)=ax^2+bx+c\), which is the same as saying that \(x_1\) and \(x_2\) are the roots. The product relation is also the key to the stable form of the next section.

For \(a=1\), \(p>0\), \(q>0\) it is: a square of side \(x\) and a strip \(p\times x\) have area \(x^2+px=q\); cutting the strip in two halves \(\frac p2\times x\) and attaching them to two sides of the square gives an L-shape of area \(q\); adding the corner \((p/2)^2\) makes a square of side \(x+p/2\) and area \(q+(p/2)^2\), so \(x+\frac p2=\sqrt{q+p^2/4}\). That is exactly the formula for the positive root, \(x=-\frac p2+\sqrt{\frac{p^2}{4}+q}\). For other signs the pieces have negative lengths and the picture stops being a picture, while the algebra does not care.

When \(b^2\gg|4ac|\), \(\sqrt{b^2-4ac}\) is almost \(|b|\), and for one of the two signs the formula computes \(-b+\sqrt\Delta\) (for \(b>0\)) as the difference of two nearly equal numbers. Each of them carries a rounding error of about \(\varepsilon|b|\), with \(\varepsilon=2^{-53}\approx1.1\times10^{-16}\), and the subtraction cancels the leading digits and leaves that error in front: the absolute error of the result is still about \(\varepsilon|b|/(2|a|)\), while the root itself is only \(|c/b|\). The relative error is of order \(\varepsilon\,b^2/|ac|\). For \(x^2+10^8x+1\): the small root is \(-10^{-8}\), the formula in double precision returns \(-7.45\times10^{-9}\) (error 25 %), and for \(b\ge2^{27}\approx1.34\times10^8\) it returns exactly 0. The other root, \(-b-\sqrt\Delta\), has no cancellation and is fine.

Compute first the root that has no cancellation, by choosing the sign that adds magnitudes: \(q=-\tfrac12\bigl(b+\operatorname{sgn}(b)\sqrt\Delta\bigr)\), so that \(q\) has the size of \(b\) and no digits are lost, and then \(x_1=q/a\). The other root comes from Viète's product \(x_1x_2=c/a\): \(x_2=\frac{c/a}{q/a}=c/q\). Both are divisions and sums of numbers of equal sign, so each has a relative error of a few units of \(\varepsilon\). In the table above the stable column is correct to the last printed digit in every row. This is what libraries do: the standard routines for quadratics (and for solving small polynomial systems) use this form, not the one in the textbook.

First, overflow: \(b^2\) exceeds the largest double (about \(1.8\times10^{308}\)) for \(|b|>1.3\times10^{154}\) even if the roots are modest; scaling all coefficients by a power of two fixes it without changing the roots. Second, a discriminant that is nearly zero: \(b^2\) and \(4ac\) nearly equal, so \(\Delta\) loses digits, and near a double root the roots are ill-conditioned, a perturbation \(\delta\) of the coefficients moves them by about \(\sqrt\delta\), so \(10^{-16}\) turns into \(10^{-8}\) no matter how the formula is written. A compensated computation of \(b^2-4ac\) (with a fused multiply-add) recovers the lost digits when the coefficients themselves are exact. Third, \(a\) tiny: with \(a=10^{-9}\) in the preset, the formula for the root that tends to \(-c/b\) loses digits, and the stable form does not.

Set \(p=4\) and \(q=5\) in the area picture. The positive root is \(x=-2+\sqrt{4+5}=1\). Check it against the picture: the square has side 1, the two strips are \(2\times1\), the corner is \(2^2=4\), the big square has side \(x+2=3\) and area \(9=q+4\). Now take \(p=3\), \(q=10\) and compute \(x\) and the big square by hand before looking at the readout. Why does the readout never ask for \(q<0\)?

In the lab, type \(a=1\), \(c=1\) and \(b=10^6\), then \(10^7\), \(10^8\) and \(1.5\times10^8\). Compare the formula as written with the exact root and note where the relative error crosses 1 % and where it becomes exactly 1. Check the rule of thumb \(\varepsilon\,b^2/|ac|\): for \(b=10^6\) it predicts \(10^{-4}\); for \(b=10^8\) it predicts about \(1\). Then try \(a=1\), \(b=-10^8\), \(c=1\): which root is now the bad one, and why?

Take \(x^2+2x+5\), whose discriminant is \(-16\). The lab prints \(x=-1\pm2i\); check that \((x+1)^2=-4\) in the form of the fifth frame, and that \((-1+2i)+(-1-2i)=-2=-b/a\) and \((-1+2i)(-1-2i)=5=c/a\). What is the vertex of the parabola, and how does its height tell you the roots are not real?