Past exam of the numerical analysis course of the University of Oxford / 2021 (A7, A13759W1)

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$.

Set $W=\mathrm{diag}(w_0,\dots ,w_n)$. For column indices $j,k\in \{1,\dots ,n+1\}$, the product $P_{j-1}{P}_{k-1}$ has polynomial degree at most $2n$, so exactness of Gauss-Legendre quadrature and orthonormality giveThus $Q=W^{1/2}A$ is a square matrix with $Q^{T}Q=I$, and therefore is an orthogonal matrix: $Q{Q}^T=I$ as well. It follows thatEquivalently, the rows of the generalized Vandermonde matrix are orthogonal andThe normalization matters: here $P_{n+1}=\sqrt{(2n+3)/2}{L}_{n+1}$, where $L_{n+1}$ is the conventional Legendre polynomial used above. Thus a formula containing only the nodes and the question's orthonormal polynomials isThis identity is the discrete orthogonality of Gaussian quadrature.

Throughout Question 2(a) and 2(b), write $U=[u_1,\dots ,u_n]$, $V=[v_1,\dots ,v_n]$, and $\Sigma =\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_n)$. Then the supplied full singular value decomposition can also be written $A=U\Sigma V^T$, with $U^{T}U=V^{T}V=I_n$ and $[U,U_{\perp }]$ an $m\times m$ orthogonal matrix.

All $n$ singular values are positive, so $\Sigma$ and $V$ are invertible matrices and $A$ has matrix rank $n$. The image of a linear map is thereforewith those columns furnishing an orthonormal basis. Also $Ax=0$ implies $\Sigma V^{T}x=0$ and hence $x=0$. ThusFinally $A^{T}z=V\Sigma U^{T}z=0$ exactly when $U^{T}z=0$. The remaining columns of the full orthogonal matrix span this orthogonal complement:If $m=n$, this last null space also equals $\{0\}$ and has the empty basis.

For $m=n$, $U$ and $V$ are square orthogonal matrices. Multiply the linear system successively by $U^T$, ${\Sigma }^{-1}$ and $V$ to obtainWith the singular value decomposition already supplied, compute $c=U^{T}b$, divide $c_j$ by ${\sigma }_j$, then form $x=Vc$. The two dense matrix-vector products each cost $O(n^2)$ arithmetic operations in Big O notation, and the diagonal scaling costs $O(n)$. ThusIf the singular value decomposition first has to be computed, its usual dense cost is $O(n^3)$; that is a separate setup cost and can be reused for further right-hand sides.

The Euclidean norm is invariant under an orthogonal matrix, so the full singular value decomposition givesThe second term is independent of $x$; the first attains its unique minimum zero because $\Sigma$ is an invertible matrix. The linear least squares solution is thereforeHere $A^{+}=V{\Sigma }^{-1}{U}^T$ is the Moore-Penrose inverse, and the fitted vector is the orthogonal projection $A{x}_{\ast }=U{U}^{T}b$ onto the image of a linear map defined by $A$.

With the singular value decomposition supplied, forming $U^{T}b$ costs $O(mn)$ arithmetic operations, diagonal scaling costs $O(n)$, and multiplication by $V$ costs $O(n^2)$. Since $m\ge n$,One need not form $U_{\perp }$, $A^{+}$ or $A^{T}A$ to calculate $x_{\ast }$. Computing the needed thin singular value decomposition first has standard dense cost $O(m{n}^2)$; explicitly forming an entire $m\times m$ orthogonal completion can require additional work when $m$ is much larger than $n$.

Retain the first $r$ terms in the singular value decomposition:The Eckart-Young-Mirsky theorem says that this truncated singular value decomposition is a best low-rank approximation in the spectral norm. The remaining singular values of the difference are ${\sigma }_{r+1},\dots ,{\sigma }_n$, together with zeros, soFor $r\ge 1$ the displayed matrix has matrix rank exactly $r$, since every retained singular value is positive. For $r=0$, the only matrix of matrix rank zero is the zero matrix, and the same formula gives error ${\sigma }_1$.

Take the identity matrixBoth $B_1$ and $B_2$ have matrix rank one, andThey are distinct and both attain the minimum spectral norm error from the preceding part. Hence a best low-rank approximation of rank one need not be unique. More generally, every orthogonal projection $u{u}^T$ with $\Vert u{\Vert }_2=1$ gives another such minimizer.

No. Distinct singular values do not force uniqueness in the spectral norm. For an explicit example takeThe singular values of $A$ are $2>1>0$. Every $B_t$ has matrix rank one, and its error in the spectral norm isThus there are infinitely many best rank-one approximations even with a strict gap between the singular values.

The same mechanism works for every $1\le r<n$ under the question's positive singular values. Starting with $A_r$, defineWhen the singular values are distinct, ${\sigma }_1+t\ge {\sigma }_1-{\sigma }_{r+1}>0$, so every $B_t$ has matrix rank $r$, whileThe spectral norm measures only the largest remaining singular value, leaving room to change retained components without increasing that maximum. This is the nonuniqueness of spectral-norm low-rank approximation. The exceptional case $r=0$ is unique for the elementary reason given above.

We use the spectral theorem to construct a sufficiently accurate approximate eigenpair. Complete $x$ to an orthonormal eigenbasis $x,x_2,\dots ,x_n$ of the symmetric matrix $B$, with corresponding eigenvalues $\lambda ,{\lambda }_2,\dots ,{\lambda }_n$. Simplicity of $\lambda$ means that every ${\lambda }_j-\lambda$ below is nonzero. PutThe orthogonal decomposition of $Ex$ shows that $x^{T}v=0$ andWith $u_{\epsilon }=x+\epsilon v$ and ${\mu }_{\epsilon }=\lambda +\epsilon a$, direct expansion cancels every first-order term:Also $\Vert u_{\epsilon }{\Vert }_2^2=1+{\epsilon }^2\Vert v{\Vert }_2^2\ge 1$. Thus its eigenpair residual has Euclidean norm $O({\epsilon }^2)$ in Big O notation.

For completeness, a real symmetric matrix $M$ always has an eigenvalue within $\Vert (M-\mu I)u{\Vert }_2/\Vert u{\Vert }_2$ of any real $\mu$ and nonzero $u$. Indeed, writing $u={\sum }_j{c}_j{q}_j$ in an orthonormal eigenbasis with eigenvalues ${\theta }_j$ givesApply this symmetric eigenvalue residual bound to $M=B+\epsilon E$. There is an eigenvalue ${\lambda }_{\epsilon }$ such thatwhich provesThis establishes the first-order perturbation of a simple symmetric eigenvalue without presupposing differentiability of the perturbed eigenvector. If $n=1$, the sum defining $v$ is empty and the displayed eigenvalue formula is exact.

Use $x_n=nh$, $h=1/N$, and $y_n$ for the numerical approximation in the scalar initial value problem $y^{\prime }=f(x,y)$, $y(0)=y_0$, on $[0,1]$.

Here $y^{\prime }=x/(1+y^2)$ and $y(0)=1$. Since $x\ge 0$, the solution is nondecreasing, so $y(x)\ge 1$. In fact $0\le y^{\prime }\le 1/2$ on $[0,1]$, which also bounds the solution and ensures it extends across this entire interval. The smooth right-hand side makes $y$ twice continuously differentiable, with second derivativeFor $y\ge 1$ and $0\le x\le 1$,so $y^{\prime \prime }\ge 0$. The subtracted term in the first expression is nonnegative, giving the useful boundThe integral form of the Taylor remainder now givesTherefore the normalized local truncation error satisfies the stronger estimateThis proves the requested bound for $n=1,\dots ,N-1$, and also for $n=0$, which will be needed to start the global error estimate. For $N=1$ the stated range is empty, but the estimate at $n=0$ still applies.

For $0\le x\le 1$ and every real $z$,where $2|z|\le 1+z^2$ was used. The mean value theorem therefore proves Lipschitz continuity in $y$ with constant $1$.

Let $e_n=y(x_n)-y_n$. The explicit Euler method and the exact-solution identity from part (i) giveUsing the requested local truncation error bound $|{\tau }_n|\le h/2$, including at $n=0$, yieldsIteration and the finite geometric series giveSince $1+h\le e^h$ and $Nh=1$,Keeping the stronger bound $|{\tau }_n|\le h/4$ from part (i) improves the same proof to $(e-1)h/4$. The requested inequality follows already from the weaker bound.

For $f(x,y)=-x{y}^2$, the explicit Euler method evaluates the right-hand side at the known point $(x_n,y_n)$. ThereforeThis update requires no root finding.

The Backward Euler method uses the new point:Writing $a=h(x_n+h)>0$ gives the quadratic equation $a{y}_{n+1}^2+y_{n+1}-y_n=0$. Starting from $y_n>0$, its two polynomial roots have opposite signs. The positive root, which continues the solution as $h\to 0$, isThe other root is negative and does not approach $y_n$ as $h\to 0$. Rationalizing the numerator gives an equivalent form that avoids cancellation error when the step size is small:The defining quadratic equation gives $0<y_{n+1}<y_n$ for this root, so positivity is maintained inductively from $y_0=1$. The branch choice is part of specifying the Backward Euler method here: the algebraic update alone has two solutions.

With $f(x,y)=-x{y}^2-y^5$, write $t=x_{n+1}=x_n+h$ and $z=y_{n+1}$. The Backward Euler method requiresor the monic quintic polynomial equationConstruct its $5\times 5$ companion matrixIts characteristic polynomial is $\det (zI-C_n)=p(z)$. Thus the polynomial roots needed for this update are exactly the eigenvalues of $C_n$.

At each step, assemble $C_n$, compute its eigenvalues, for example with the QR algorithm, and choose its unique positive real eigenvalue as $y_{n+1}$. To justify this selection, letIf $y_n>0$, then $F(0)=-y_n<0$, $F(y_n)=h{y}_n^5+ht{y}_n^2>0$, andThe intermediate value theorem and strict monotonicity prove there is exactly one positive polynomial root, located in $(0,y_n)$. Starting at $y_0=1$ therefore makes every update well defined, positive and decreasing. In numerical arithmetic, select the nearly real eigenvalue in this interval and check the polynomial residual $|F(z)|$; a final Newton root-finding iteration can refine it. This reduces each update to the requested eigenvalue problem.

No. A necessary condition for consistency of a numerical method is reproduction of the constant solution $y(x)=1$ of the ordinary differential equation $y^{\prime }=0$. In a linear multistep method this givesBut the reverse triangle inequality gives, under the strict coefficient inequality,Thus the necessary order conditions for a linear multistep method already fail at order zero, regardless of the values of the ${\beta }_j$.

As appropriate for a genuine $k$-step linear multistep method, take ${\alpha }_k\ne 0$. Suppose the characteristic polynomial has a polynomial root $\xi$ with modulus $|\xi |>1$. From $\rho (\xi )=0$ and the triangle inequality,Dividing by the positive number $|{\alpha }_k||\xi {|}^{k-1}$ gives $|\xi |\le 1$, a contradiction. HenceThe nonzero leading coefficient matters: if ${\alpha }_k=0$, the assumed inequality forces all the ${\alpha }_j$ to vanish, so $\rho$ is the zero polynomial rather than a characteristic polynomial with $k$ roots. This degenerate case does not define the method in the question.

Yes: for $k=2$ the stated inequality guarantees zero-stability. Part (ii) puts every polynomial root in the closed unit disc. The remaining root condition for a multistep method is that every root on the unit circle has multiplicity of a root one.

If this failed for a quadratic polynomial, it would have a repeated root $\xi$ with $|\xi |=1$, andConsequently $|{\alpha }_1|=2|{\alpha }_2|$ and $|{\alpha }_0|=|{\alpha }_2|$, contradictingbecause ${\alpha }_2\ne 0$. Thus all roots on the unit circle are simple andThis conclusion does not assert consistency of a numerical method: part (i) rules out consistency of a numerical method whenever the coefficient inequality is strict. Zero-stability and consistency of a numerical method are separate requirements.