Numerical analysis course of the University of Oxford (A7)

Source: Oxford Part A, A7 Numerical Analysis, A13759W1, 23 June 2023; parent PDF pages 135–137, paper pages 1–3. These original solutions cover all three questions and all eleven lettered subparts. In Question 1(d), the unrefined Gershgorin discs alone are inconclusive for the shift $5$; a further exact calculation resolves it. In Question 2(c), clearing denominators produces a linear system whose solution need not define a rational interpolant at every interpolation node. Question 3 uses the paper's normalized convention for the local truncation error.

For a complex $n\times n$ matrix $B=(b_{ij})$, putThe Gershgorin circle theorem states thatTo prove it, let $Bv=\lambda v$ with $v\ne 0$, and choose $i$ for which $|v_i|={\max }_j|v_j|>0$. The $i$th coordinate of the eigenvalue equation givesThe triangle inequality, followed by division by $|v_i|$, givesThus $\lambda \in D_i$, proving the assertion for every eigenvalue.

The second theorem, Gershgorin eigenvalue counting, says that if a union of $k$ indexed Gershgorin discs is disjoint from the remaining $n-k$ Gershgorin discs, that union contains exactly $k$ eigenvalues, counted with algebraic multiplicity. In particular,Discs that touch are not disjoint; overlapping Gershgorin discs cannot individually be assigned one eigenvalue each.

The power iteration starts from $v_0\ne 0$, normally scaled to unit Euclidean norm, and repeatsNormalization prevents the rapidly growing or decaying magnitudes of $A^k{v}_0$ from obscuring its direction. The iteration requires $w_{k+1}\ne 0$.

Because $A$ is a diagonalizable matrix, choose an eigenbasis $u_1,\dots ,u_n$, with corresponding eigenvalues ${\lambda }_1=\lambda ,{\lambda }_2,\dots ,{\lambda }_n$, and expand $v_0={\sum }_j{c}_j{u}_j$. First suppose the dominant eigenvalue is simple and $c_1\ne 0$. ThenForthe sum in parentheses after $u_1$ is $O(q^k)$ in Big O notation. Consequently the power iteration direction error is $O(q^k)$: the iterates approach the line spanned by $u_1$. For a real matrix, a uniquely dominant eigenvalue is real, since otherwise its complex conjugate would have the same modulus. If $\lambda >0$, $v_k$ approaches one unit eigenvector; if $\lambda <0$, successive $v_k$ alternate signs. In either case,This is the sign and phase in power iteration qualification. Also, the necessary assumption is a nonzero coefficient $c_1$ in the eigenbasis, rather than merely a nonzero ordinary inner product with $u_1$: those tests can differ for a nonsymmetric matrix. If $c_1=0$, the method cannot generate the missing component.

If the dominant value $\lambda$ has a higher-dimensional eigenspace, the same argument replaces $c_1{u}_1$ by the nonzero sum of all starting components in that eigenspace. The limiting direction is then that of this sum, which is still an eigenvector for $\lambda$. The usual algorithm assumes $\lambda \ne 0$; a strict dominance gap guarantees this whenever another eigenvalue is present.

An estimate retaining the sign of $\lambda$ is the Rayleigh quotientContinuity of this quotient at a nonzero eigenvector proves the limit even when $A$ is a nonsymmetric matrix. Alternatively, $(A{v}_k{)}_i/(v_k{)}_i\to \lambda$ for a coordinate whose limiting eigenvector component is nonzero. The scaling factor $\Vert A{v}_k{\Vert }_2$ alone tends only to $|\lambda |$.

Choose a shift $\mu$ which is not an eigenvalue. Inverse iteration applies power iteration to $B=(A-\mu I{)}^{-1}$, implemented by a linear system solve rather than by forming an explicit matrix inverse:A single LU decomposition of $A-\mu I$ can be reused for all the solves at this fixed shift. An eigenpair $A{u}_j={\lambda }_j{u}_j$ transforms intoThus an eigenvalue ${\lambda }_{\ast }$ uniquely closest to $\mu$ becomes the dominant eigenvalue of $B$. For a diagonalizable matrix, the argument in the preceding subpart gives convergence of diagonalizable fixed-shift inverse iteration, provided the starting coefficient in the target eigenspace is nonzero, with power iteration direction error $O(q_{\mu }^k)$, whereFor any prescribed simple eigenvalue, taking a sufficiently close but unequal shift makes the numerator small while the denominator stays close to its positive separation from the remaining eigenvalues. This makes $q_{\mu }$ small and gives very fast convergence per iteration. The Rayleigh quotient of $v_k$ estimates ${\lambda }_{\ast }$ directly. If ${\beta }_k$ estimates the associated eigenvalue of $B$, the corresponding estimate is $\mu +{\beta }_k^{-1}$.

The selection principle also applies when the other eigenvalues have nontrivial blocks in the Jordan normal form: their contributions acquire powers of $k$, multiplying strictly smaller geometric factors. A simple uniquely closest eigenvalue still attracts the direction when the Riesz projection of $v_0$ onto its generalized eigenspace is nonzero. For nonreal target eigenvalues, allow complex shifts and arithmetic; a real shift cannot distinguish a conjugate pair of equal distances. Setting $\mu ={\lambda }_{\ast }$ makes the solve singular, and an exactly zero starting target component prevents convergence to that target in exact arithmetic. When ${\lambda }_{\ast }<\mu$ is real, account for the sign behavior of fixed-shift inverse iteration.

The four row Gershgorin discs have centres $9,3,-4,-5$ and radii $1,2,2,1$, respectively. Their real-axis intersections areNone of the Gershgorin discs contains zero, so the Gershgorin circle theorem givesThe discs centred at $9$ and $3$ are each isolated. By Gershgorin eigenvalue counting, they contain exactly one eigenvalue apiece, and both eigenvalues are simple. Each must be real: the characteristic polynomial is real, so a nonreal root would bring a distinct conjugate root into the same disc. Denote these eigenvalues by $\beta \in [8,10]$ and $\alpha \in [1,5]$. The remaining two eigenvalues lie in the union of the negative-centred discs.

For the shift $4$, the distances satisfyfor either of the other eigenvalues $\gamma$, since its real part is at most $-2$. Thus $\alpha$ is uniquely closest to $4$. Once $4$ is checked not to be an eigenvalue, the natural full target-disc interval enclosure and a conservative asymptotic direction factor areAs always, inverse iteration requires a nonzero initial target component. The following calculation verifies the shift and gives sharper information.

For the tridiagonal matrix, let $p_j(t)$ be the determinant of the leading $j\times j$ block of $tI-A$. Expansion along the last row givesIn particular, $\det A=577$, $p_4(4)=-239$ and $p_4(5)=-583$. Both shifted linear systems are nonsingular. Further exact evaluations areThe intermediate value theorem gives a real root in each of $(-6,-5)$, $(-4,-3)$, $(3,4)$ and $(8,9)$. As $p_4$ has polynomial degree four, these are all its roots and each is simple. In particular, $A$ is a diagonalizable matrix, $\alpha \in (3,4)$ and $\beta \in (8,9)$. The closest competitor to $\alpha$ at the shift $4$ is $\beta$, so this refinement yieldsFor the shift $5$, the original Gershgorin discs only give $|\alpha -5|\le 4$ and $|\beta -5|\ge 3$; those bounds do not decide which is closer. The refined interval enclosures do:for the negative eigenvalues. Therefore the answer for this particular matrix is yes:Both shifts select the same eigenvector line and alternate signs because $\alpha <\mu$. The factor bounds concern asymptotic directional convergence, rather than a universal error reduction from the first iterate for a nonsymmetric matrix. The wider interval $[1,5]$ remains a valid full-disc enclosure for both shifts; $[3,4]$ records the stronger information obtained from the characteristic polynomial.

For the distinct interpolation nodes, define the Lagrange cardinal polynomialsEvery denominator is nonzero, each $L_i$ has polynomial degree at most $n$, and $L_i(x_j)={\delta }_{ij}$. Hence the Lagrange interpolation polynomialbelongs to ${\Pi }_n$ and has $p_n(x_j)=f_j$. This proves existence.

If $p_n$ and ${\tilde{p}}_n$ both fit the data, their difference has polynomial degree at most $n$ but vanishes at the $n+1$ distinct interpolation nodes. Repeated application of the factor theorem implies that ${\prod }_{j=0}^n(x-x_j)$ divides that difference. A nonzero polynomial of degree at most $n$ cannot have such a divisor of degree $n+1$. Thus $p_n-{\tilde{p}}_n=0$, proving uniqueness. When $n=0$, the empty product defining $L_0$ is $1$, so the argument includes the constant case.

The monomial $x^k$ lies in ${\Pi }_n$ and already fits the data. Uniqueness of polynomial interpolation therefore givesApply this result with a nonnegative integer $m\le n$. Its Lagrange interpolation representation is the polynomial identityThe numerator in each summand is a monic polynomial of degree $n$, so the coefficient of $x^n$ on the right is exactly ${\sum }_i{x}_i^m/{\prod }_{j\ne i}(x_i-x_j)$. The coefficient on the left is zero for $m<n$ and one for $m=n$. Comparing coefficients proves the Lagrange leading-coefficient identityThe intended exponent range is $m\ge 0$: arbitrary negative exponents would not give this polynomial identity. For $m=0$, $x_i^0$ denotes the value of the constant polynomial $1$, including at a zero interpolation node.

Clearing the denominators in the rational interpolation conditions givesThese are $n+1$ linear equations for $(m+1)+\ell =n+1$ unknowns. The $i$th row of the coefficient matrix, in the order $a_0,\dots ,a_m,b_0,\dots ,b_{\ell -1}$, isThis rational interpolation linearization is equivalent to the original fraction conditions only if $r(x_i)\ne 0$ at every interpolation node. If $r(x_i)=0$, the cleared equation merely forces $q(x_i)=0$, leaving the fraction undefined.

For $f_i=x_i^{m+1}$, the coefficient row becomesThis is the Vandermonde matrix, whose determinant is ${\prod }_{0\le i<j\le n}(x_j-x_i)\ne 0$. Equivalently, a vector in its null space would define a degree-at-most-$n$ polynomial vanishing at all $n+1$ distinct interpolation nodes, which is zero by the preceding uniqueness proof. Thus the cleared linear system has exactly one solution.

There is also a constructive description of that solution. Define the monic polynomialThe polynomial $x^{m+1}r(x)-q(x)$ is monic of degree $n+1$ and vanishes at every interpolation node. It must therefore equal $W(x)$. By polynomial division, write uniquelyThe quotient $R$ is monic of degree $\ell$, andis the unique coefficient solution. This does not assert that $R$ is nonzero at every interpolation node.

For a counterexample with $m=0$, $\ell =1$, choose $x_0=0,x_1=1$ and the specified data $f_i=x_i$, so $f_0=0,f_1=1$. The cleared equations reduce togiving the unique solution $q(x)=0$, $r(x)=x-1$. But at $x_1=1$ the quotient is $0/0$. Indeed, any genuine rational interpolant of this prescribed form would need $r(0)\ne 0$; then $q(0)/r(0)=0$ forces the constant numerator to be zero. At a nonzero denominator its value at $1$ would also be zero, never one. Cancelling the common factor in the cleared candidate would merely give the identically zero rational function, which still fails the data. Hence

Integrate the ordinary differential equation across one step, using the fundamental theorem of calculus:The midpoint quadrature rule approximates the integral by $hf(x_n+h/2,y(x_n+h/2))$. Approximate the unknown midpoint state by the midpoint value of linear interpolation between the endpoint states:Combining these approximations gives the implicit midpoint ruleIt is an implicit time-stepping method because its new state occurs inside the evaluation of $f$. An equivalent derivation uses a straight-line approximation on the step and imposes the ordinary differential equation at its midpoint, a one-stage collocation method.

Work on a fixed finite interval $[a,a+T]$. Write $Y_n=y(x_n)$, $t=x_n+h/2$, and use the paper's normalized local truncation errorAssume $y$ has bounded continuous derivatives through order three on this interval, with $\Vert y^{\prime \prime }\Vert \le M_2$ and $\Vert y^{\prime \prime \prime }\Vert \le M_3$. Assume also that the Jacobian matrix $f_y$ has operator norm at most $L$ on the segments joining $y(t)$ to the endpoint averages. For a scalar equation this last hypothesis is simply $|\partial f/\partial y|\le L$ there.

Two symmetric estimates separate the errors. By Taylor's theorem applied to $y^{\prime }$ about $t$, the linear term integrates to zero on $[-h/2,h/2]$, givingLikewise the two linear terms in the endpoint expansions cancel, and their second-order remainders giveThe derivative bound on $f$ gives a Lipschitz condition in its second variable. Since $y^{\prime }(t)=f(t,y(t))$, the triangle inequality now yieldsThe constant is independent of $n$ and $h$ as long as the steps remain in the fixed interval. This implicit midpoint consistency bound proves numerical consistency of order two. The corresponding unnormalized residual is $h{\tau }_n=O(h^3)$; the two powers conventionally called local truncation error differ because the paper divides by the step size. If local error instead means the error of an exactly solved numerical step started at $Y_n$, the small-step Lipschitz condition gives the bound $h\Vert {\tau }_n\Vert /(1-hL/2)=O(h^3)$ as well.

Let $f$ be continuous in its first variable and uniformly Lipschitz continuous in its second variable, with constant $L\ge 0$. Retain the finite-interval smoothness assumptions and the bound $\Vert {\tau }_n\Vert \le K{h}^2$ from the preceding subpart. We first establish solvability of each implicit step. For a proposed starting state $u$, the maphas Lipschitz constant $hL/2$. If $hL/2<1$, the contraction mapping theorem on the complete state space gives a unique next state $v={\Psi }_{n,h}(u)$.

To compute a useful Lipschitz constant for the method, introduce its increment function byFor two starting states $u,w$,Rearranging gives the implicit midpoint Lipschitz estimateThus, for example, $hL\le 1$ gives an increment Lipschitz constant at most $2L$, uniformly as $h\to 0$.

For an explicit global error estimate, set $e_n=Y_n-y_n$. Subtract the numerical formula from the exact-value identity containing $h{\tau }_n$:The Lipschitz condition and the triangle inequality implyand thereforeIterating this inequality is the elementary discrete Gronwall inequality argument. For $L>0$, summing the geometric series and using $R_h-1=hL/(1-hL/2)$ givesFor $hL\le 1$, $R_h=1+hL/(1-hL/2)\le 1+2hL\le e^{2hL}$. If $nh\le T$, this provesWith the exact initial value, $e_0=0$, the global error is $O(h^2)$ and tends uniformly to zero on every such finite interval. If $L=0$, the recurrence instead gives $\Vert e_n\Vert \le \Vert e_0\Vert +KT{h}^2$, with the same conclusion. Starting values merely converging to the exact initial value also suffice for numerical convergence; second-order accuracy requires their errors to be $O(h^2)$. The result concerns exact implicit solves and convergence on fixed finite intervals, not uniform accuracy over an unbounded time interval.

Applying the implicit midpoint rule to the simple harmonic oscillator system givesFactor the change in the quadratic form $y^2+z^2$ and substitute these equations:Thus every step preserves the first integralFor completeness, the implicit equations are solvable for every real $h$: putting $s=h/2$, their coefficient matrix has determinant $1+s^2>0$ and the update isThis is an orthogonal matrix, as follows from $(1-s^2{)}^2+4s^2=(1+s^2{)}^2$. It therefore preserves the Euclidean norm, explaining geometrically the exact energy conservation. This is a concrete instance of Runge-Kutta conservation of quadratic invariants. Preservation of the circle does not assert that the numerical and exact solutions have the same phase at a fixed nonzero step size.

Source: Oxford Part A, A7 Numerical Analysis, A13759W1, 24 June 2021; parent PDF pages 128–131, paper pages 1–4. These original solutions cover all three questions, all nine lettered parts and all twenty-one answer units, including every optional question. Each question carries 25 marks; the examination counted the best two answers. The local truncation error in Question 3(a) is divided by the step size. In Question 3(c), the leading coefficient of a genuine $k$-step linear multistep method is understood to be nonzero.

There are $n+2$ distinct interpolation nodes but only $n+1$ free coefficients in a polynomial of polynomial degree at most $n$. Thus polynomial interpolation is possible precisely when the data satisfy one compatibility equation. To identify it, form the Lagrange interpolation polynomial through all the data:This is the unique interpolation polynomial of polynomial degree at most $n+1$. Its leading coefficient, which is also the order-$(n+1)$ divided difference, isConsequently the required condition isNecessity follows because any qualifying $p_n$ must equal $q$ by uniqueness of polynomial interpolation; sufficiency follows because $D=0$ removes the only possible term of polynomial degree $n+1$ from $q$.

When this condition holds, a slightly shorter expression uses only the first $n+1$ interpolation nodes:This polynomial agrees with $q$ at $n+1$ distinct points, so the root bound for a polynomial gives $p_n=q$ and hence $p_n(x_{n+1})=f_{n+1}$ as required. Equivalently, the compatibility condition says that this expression, evaluated at the last interpolation node, gives $f_{n+1}$.

Write ${\mathcal{P}}_n$ for the vector space of polynomials of polynomial degree at most $n$, and put $p_n(x)={\sum }_{j=0}^n{a}_j{x}^j$. As the question calls the weighted integral an inner product, its finiteness is implicit; for example, a positive integrable weight $w$ suffices for all the continuous functions involved here.

For every $q\in {\mathcal{P}}_n$, differentiating the squared weighted norm $\Vert f-p_n-tq{\Vert }_w^2$ at its minimum $t=0$ gives the orthogonality equationsIn particular, the requested linear system isIts Gram matrix $G_{kj}={\int }_{-1}^{1}w(x)x^{k+j}dx$ is a symmetric positive-definite matrix: for any nonzero vector $a$,Indeed a nonzero polynomial is nonzero on some open interval, and the weight is positive. Therefore $G$ is an invertible matrix, proving uniqueness. Conversely, if these equations hold, the Pythagorean theorem in an inner-product space gives, for every $q\in {\mathcal{P}}_n$,so the solution really is the unique best approximation.

Let $r=f-p_n$ be the continuous function describing the error. The preceding orthogonality equations sayIf $r$ vanishes identically, any $n+1$ distinct points work. If $r$ already has at least $n+1$ distinct zeros, the assertion also follows. Suppose instead that it has at most $n$ distinct zeros. Among its interior zeros, list those at which its sign changes as $t_1<\cdots <t_m$, where $m\le n$, and define the polynomialwith $q=1$ if $m=0$. Between consecutive zeros of $r$, its sign is constant by the intermediate value theorem. At each $t_j$, both $r$ and $q$ change sign; at every other zero of $r$, neither sign changes. Hence $rq$ has one fixed sign throughout $[-1,1]$ wherever it is nonzero. Replacing $q$ by $-q$ if necessary makes $rq\ge 0$ everywhere.

Since $r$ is not identically zero and $q$ has only finitely many zeros, $rq$ is strictly positive on some open interval. Positivity of $w$ givescontradicting orthogonality because $q\in {\mathcal{P}}_n$. ThusThis is the interpolation property of weighted least-squares polynomial approximation. It supplies interpolation nodes depending on $f$ and $w$, rather than nodes prescribed in advance.

Let $L_{n+1}$ denote the conventional Legendre polynomial, normalized by $L_{n+1}(1)=1$. The $n+1$ quadrature nodes are its distinct polynomial roots in $(-1,1)$. For these nodes form the Lagrange cardinal polynomialsChoose the quadrature weights by integrating these polynomials:An equivalent standard formula isThe reason for selecting these quadrature nodes is the resulting degree of precision of a quadrature rule. For any polynomial $p$ of polynomial degree at most $2n+1$, polynomial division givesThe orthogonality of the Legendre polynomial makes ${\int }_{-1}^1{L}_{n+1}s=0$, while $p(x_i)=r(x_i)$ at each quadrature node. The definition of the quadrature weights integrates $r$ exactly, soIn particular, ${\sum }_i{w}_i=2$. Applying this exactness to ${\ell }_i^2$ also proves $w_i={\int }_{-1}^1{\ell }_i(x{)}^{2}dx>0$.

Truncate the Chebyshev series at the full degree of precision of a quadrature rule:Since $|T_j(x)|\le 1$ on $[-1,1]$, the geometric series estimate givesThe same bound proves uniform convergence of the original Chebyshev series, so its sum and the tail estimate are justified. Gauss-Legendre quadrature is exact on $s_{2n+1}$. Its positive quadrature weights sum to $2$, henceTherefore one explicit choice proving exponential convergence of Gaussian quadrature isUsing exactness through polynomial degree $2n+1$ doubles the exponent compared with truncating at polynomial degree $n$.