Numerical Methods Qualification Exam Problems and Solutions (University of Maryland)/August 2005

Problem 1

Derive the one-, two-, and three-point Gaussian quadrature formulas such that $\underset{-1}{\overset{1}{\int }}f(x)x^{4}dx\approx \sum _{j=1}^{n}f(x_j)w_j$ Give Bounds on the error of these formulas.

Solution 1

Find Orthogonal Polynomials w.r.t. Weighting Function

For this problem, we need the first three orthogonal polynomials $P_1,P_2,\text{ and }{P}_3$ with respect to a weighted inner product

$⟨f,g⟩=\underset{-1}{\overset{1}{\int }}f(x)g(x)w(x)dx$

where $w(x)=x^4$ in our case. This can be done using Gram-Schmidt Orthogonalization on the basis $\{1,x,x^2,x^3\}$

Zero Order Polynomial

Let $P_0(x)=1$

First Order Polynomial

$\begin{array}{cc}{P}_1(x)& =x-\frac{⟨1,x⟩}{⟨1,1⟩}\cdot 1\\ & =x-\underset{\frac{0}{2/5}}{\underbrace{\frac{\underset{-1}{\overset{1}{\int }}{x}^{5}dx}{\underset{-1}{\overset{1}{\int }}{x}^{4}dx}}}\\ & =x\end{array}$

Second Order Polynomial

$\begin{array}{cc}{P}_2(x)& =x^2-\frac{⟨1,x^2⟩}{⟨1,1⟩}\cdot 1-\frac{⟨x,x^2⟩}{⟨x,x⟩}\cdot x\\ & =x^2-\frac{\underset{-1}{\overset{1}{\int }}{x}^{6}dx}{\underset{-1}{\overset{1}{\int }}{x}^{4}dx}-x\cdot \underset{\frac{0}{2/7}}{\underbrace{\frac{\underset{-1}{\overset{1}{\int }}{x}^{7}dx}{\underset{-1}{\overset{1}{\int }}{x}^{6}dx}}}\\ & =x^2-\frac{5}{7}\end{array}$

Third Order Polynomial

$\begin{array}{cc}{P}_3(x)& =x^3-\frac{⟨1,x^3⟩}{⟨1,1⟩}\cdot 1-\frac{⟨x,x^3⟩}{⟨x,x⟩}\cdot x-\frac{⟨x^2-\frac{5}{7},x^3⟩}{⟨x^2-\frac{5}{7},x^2-\frac{5}{7}⟩}\cdot (x^2-\frac{5}{7})\\ & =x^3-0-\frac{7}{2}x\cdot \frac{\underset{-1}{\overset{1}{\int }}{x}^{8}dx}{\underset{-1}{\overset{1}{\int }}{x}^{6}dx}-0\\ & =x^3-\frac{7}{9}x\end{array}$

Find Roots of Polynomials

The roots of these polynomials will be the interpolation nodes used in Gauss Quadrature.

$\begin{array}{cc}\text{Polynomial}& \text{Roots}\\ P_1(x)=x& 0\\ \\ P_2(x)=x^2-\frac{5}{7}& \pm \sqrt{\frac{5}{7}}\\ \\ P_3(x)=x^3-\frac{7}{9}x& 0,\pm \sqrt{\frac{7}{9}}\\ \end{array}$

Formula to Compute Coefficients

In Gaussian quadrature we have

$w_j=\underset{-1}{\overset{1}{\int }}{l}_i(x)x^{4}dx$, where

$l_i(x)=\prod _{j=1,j\ne i}^n\frac{x-x_j}{x_i-x_j}$

1 point formula

$l_1(x)=1$

$w_1=\underset{-1}{\overset{1}{\int }}{x}^{4}dx=\frac{2}{5}$

$\underset{-1}{\overset{1}{\int }}f(x)x^{4}dx\approx \frac{2}{5}f(x_1)$

where $x_1$ is the root of $P_1(x)$.

2 point formula

$l_1(x)=\frac{x-x_1}{x_1-x_2}$

$l_2(x)=\frac{x-x_1}{x_2-x_1}$

$w_1=\underset{-1}{\overset{1}{\int }}\frac{x-x_1}{x_1-x_2}{x}^{4}dx=\frac{2x_2}{5(x_2-x_1)}$

$w_2=\underset{-1}{\overset{1}{\int }}\frac{x-x_1}{x_2-x_1}{x}^{4}dx=\frac{2x_1}{5(x_1-x_2)}$

$\underset{-1}{\overset{1}{\int }}f(x)x^{4}dx\approx w_{1}f(x_1)+w_{2}f(x_2)$

where $x_1\text{ and }{x}_2$ are the roots of $P_2(x)$.

3 point formula

$l_1(x)=\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}$

$l_2(x)=\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}$

$l_3(x)=\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}$

$w_1=\underset{-1}{\overset{1}{\int }}\frac{x-x_2}{x_1-x_2}\frac{x-x_3}{x_1-x_3}{x}^{4}dx=\frac{1}{(x_1-x_2)(x_1-x_3)}(\frac{2}{7}+\frac{2x_2{x}_3}{5})$

$w_2=\underset{-1}{\overset{1}{\int }}\frac{x-x_1}{x_2-x_1}\frac{x-x_3}{x_2-x_3}{x}^{4}dx=\frac{1}{(x_2-x_1)(x_2-x_3)}(\frac{2}{7}+\frac{2x_1{x}_3}{5})$

$w_3=\underset{-1}{\overset{1}{\int }}\frac{x-x_1}{x_3-x_1}\frac{x-x_2}{x_3-x_2}{x}^{4}dx=\frac{1}{(x_3-x_1)(x_3-x_2)}(\frac{2}{7}+\frac{2x_1{x}_2}{5})$

$\underset{-1}{\overset{1}{\int }}f(x)x^{4}dx\approx w_{1}f(x_1)+w_{2}f(x_2)+w_{3}f(x_3)$

where $x_1\text{, }{x}_2\text{ and }{x}_3$ are the roots of $P_3(x)$.

Derive Error Bound

We know that this quadrature is exact for all polynomials of degree at most $2n-1$.

We choose a polynomial $p$ of degree at most $2n-1$ that Hermite interpolates i.e.

$p(x_i)=f(x_i)p^{\prime }(x_i)=f^{\prime }(x_i)(0\le i\le n-1)$

The error for this interpolation is

$f(x)-p(x)=\frac{f^{(2n)}(\zeta (x))}{(2n)!}\prod _{i=0}^{n-1}(x-x_i{)}^2$

Compute the error of the quadrature as follows:

$\begin{array}{cc}E& =\underset{-1}{\overset{1}{\int }}f(x)x^{4}dx-\sum _{i=0}^{n-1}{w}_{i}f(x_i)\\ & =\underset{-1}{\overset{1}{\int }}f(x)x^{4}dx-\sum _{i=0}^{n-1}{w}_{i}p(x_i)\\ & =\underset{-1}{\overset{1}{\int }}f(x)x^{4}dx-\underset{-1}{\overset{1}{\int }}p(x)x^{4}dx\\ & =\underset{-1}{\overset{1}{\int }}(f(x)-p(x))x^{4}dx\\ & =\underset{-1}{\overset{1}{\int }}\frac{f^{(2n)}(\zeta (x))}{(2n)!}\prod _{i=0}^{n-1}(x-x_i{)}^2\cdot x^{4}dx\\ & =\frac{f^{(2n)}(\xi )}{(2n)!}\underset{-1}{\overset{1}{\int }}{x}^4\cdot \underset{(p_n(x){)}^2}{\underbrace{\prod _{i=0}^{n-1}(x-x_i{)}^2}}dx\end{array}$

where the last line follows from the mean value theorem of integrals.

Notice that $p_n(x)$ is simply the polynomials orthogonal with respect to weight function since $x_i$ are the roots of the polynomials.

Hence, the error for 1 point gaussian quadrature is

$\begin{array}{cc}E& =\frac{f^{(2)}(\xi )}{(2)!}\underset{-1}{\overset{1}{\int }}{x}^4\cdot x^{2}dx\\ & =\frac{f^{(2)}(\xi )}{7}\end{array}$

The error for 2 point quadrature:

$\begin{array}{cc}E& =\frac{f^{(4)}(\xi )}{(4)!}\underset{-1}{\overset{1}{\int }}{x}^4\cdot (x^2-\frac{5}{7}{)}^{2}dx\end{array}$

The error for 3 point quadrature:

$\begin{array}{cc}E& =\frac{f^{(6)}(\xi )}{(6)!}\underset{-1}{\overset{1}{\int }}{x}^4\cdot (x^3-\frac{7}{9}x{)}^{2}dx\end{array}$

Problem 2

We wish to solve $Ax=b$ iteratively where $A=\left[\begin{array}{ccc}1& 1& -1\\ 1& 2& 1\\ 1& 1& 1\end{array}\right]$ Show that for this $A$ the Jacobi method and the Gauss-Seidel method both converge. Explain why for this $A$ one of these methods is better than the other.

Solution 2

Jacobi Method

Decompose $A$ into its diagonal, lower, and upper parts i.e.

$A=D+(L+U)$

Derive Jacobi iteration by solving for x as follows:

$\begin{array}{cc}Ax& =b\\ (D+(L+U))x& =b\\ Dx+(L+U)x& =b\\ Dx& =b-(L+U)x\\ x& =D^{-1}(b-(L+U)x)\end{array}$

Gauss-Seidel Method

Decompose $A$ into its diagonal, lower, and upper parts i.e.

$A=(D+L)+U$

Derive Gauss-Sediel iteration by solving for x as follows:

$\begin{array}{cc}Ax& =b\\ ((D+L)+U)x& =b\\ (D+L)x+Ux& =b\\ (D+L)x& =b-Ux\\ x& =(D+L{)}^{-1}(b-Ux)\end{array}$

Jacobi Convergence

Convergence occurs for the Jacobi iteration if the spectral radius of $D^{-1}(L+U)$ is less than 1 i.e.

$\rho (D^{-1}(L+U))<1$

The eigenvalues of the matrix

$D^{-1}(L+U)=\left[\begin{array}{ccc}0& 1& -1\\ \frac{1}{2}& 0& \frac{1}{2}\\ 1& 1& 0\end{array}\right]$

are $\lambda =0,0,0$ i.e. the eigenvalue $0$ has order 3.

Therefore, the spectral radius is $\rho =0<1$.

Gauss-Seidel Convergence

Convergence occurs for the Gauss-Seidel iteration if the spectral radius of $(D+L{)}^{-1}U$ is less than 1 i.e.

$\rho ((D+L{)}^{-1}U)<1$

The eigenvalues of the matrix

$(D+L{)}^{-1}U=\left[\begin{array}{ccc}0& 1& -1\\ 0& -\frac{1}{2}& 1\\ 0& -\frac{1}{2}& 0\end{array}\right]$

are $\lambda =0,-\frac{1}{4}+\frac{\sqrt{7}i}{4},-\frac{1}{4}-\frac{\sqrt{7}i}{4}$

Therefore, the spectral radius is $\rho =\sqrt{\frac{1}{16}+\frac{7}{16}}=\frac{\sqrt{2}}{2}\approx .7<1$.

Comparison of Methods

In general, the Gauss-Seidel iteration converges faster than the Jacobi iteration since Gauss-Seidel uses the new components of $x_i$ as they become available, but in this case ${\rho }_{Jacobi}<{\rho }_{Gauss-Seidel}$, so the Jacobi iteration is faster.

Problem 3

Consider the three-term polynomial recurrence $p_{k+1}(x)=(x-{\mu }_k)p_k(x)-{\nu }_k^2{p}_{k-1}(x)\text{ for }k=1,2,\dots ,$ initialized by $p_0(x)=1\text{ and }{p}_1(x)=x-{\mu }_0$, where each ${\mu }_k\text{ and }{\nu }_k$ is real and each ${\nu }_k$ is nonzero.

Problem 3a

Prove that each $p_k(x)$ is a monic polynomial of degree $k$, and for every $n=0,1,\dots$, one has $\mathrm{span}\{p_0(x),p_1(x),\dots ,p_n(x)\}=\mathrm{span}\{1,x,x^2,\dots ,x^n\}$

Solution 3a

We prove the claim by by induction.

Base Case

$p_0(x)=1$ is monic with degree zero, and $\mathrm{span}\{p_0(x)\}=\mathrm{span}\left\{1\right\}$.

$p_1(x)=x-{\mu }_0$ is monic with degree one, and $\mathrm{span}\{p_0(x),p_1(x)\}=\mathrm{span}\{1,x\}$.

$p_2(x)=(x-{\mu }_0)(x-{\mu }_1)-{\nu }_1^2$ is monic with degree 2, and $\mathrm{span}\{p_0(x),p_1(x),p_2(x)\}=\mathrm{span}\{1,x,x^2\}$.

Induction Step

Assume $p_{k-1}(x)$ is monic with degree $k-1$, and $\mathrm{span}\{p_0(x),p_1(x),\dots ,p_{k-1}(x)\}=\mathrm{span}\{1,x,\dots ,x^{k-1}\}$.

Also, assume $p_k(x)$ is monic with degree $k$, and $\mathrm{span}\{p_0(x),p_1(x),\dots ,p_k(x)\}=\mathrm{span}\{1,x,\dots ,x^k\}$.

Then from the hypothesis, $p_{k+1}(x)=(x-{\mu }_k)p_k(x)-{\nu }_k^2{p}_{k-1}(x)$ is monic with degree $k+1$.

Also,

$\mathrm{span}\{p_0(x),p_1(x),\dots ,p_{k+1}(x)\}=\mathrm{span}\{1,x,\dots ,x^{k+1}\}$

since $p_{k+1}(x)$ is just $x^{n+1}$ plus a linear combination of lower order terms already in $\mathrm{span}\{1,x,\dots ,x^k\}$.

Problem 3b

Show that for every $k=0,1,\dots$ the polynomial $p_k(x)$ has $k$ simple real roots that interlace with the roots of $p_{k-1}(x)$.

Solution 3b

We prove the claim by induction.

Base Case

Let's consider the case $k=2$. We know that

$p_2(x)=(x-{\mu }_0)(x-{\mu }_1)-{\nu }_1^2.$

The quadratic formula shows that $p_2(x)$ has two simple real roots.

Let $r_{21}$ and $r_{22}$ be the zeros of $p_2(x)$. Then, we have (because of sign changes) that

$p_2({\mu }_0)=-{\nu }_1^2$ and $\underset{x\to \mathrm{\infty }}{\lim }{p}_2(x)>0$ $\Rightarrow$ there exists $r_{22}\in ({\mu }_0,\mathrm{\infty })$ such that $p_2(r_{22})=0$.

$p_2({\mu }_0)=-{\nu }_1^2$ and $\underset{x\to -\mathrm{\infty }}{\lim }{p}_2(x)>0$ $\Rightarrow$ there exists $r_{21}\in (-\mathrm{\infty },{\mu }_0)$ such that $p_2(r_{21})=0$.

In conclusion, $r_{21}<\underset{r_{11}}{\underbrace{{\mu }_0}}<r_{22}.$.

Induction Step

Let $r_{k-1,1},\dots ,r_{k-1,k-1}$ and $r_{k,1},\dots ,r_{k,k}$ be the simple real roots of $p_{k-1}(x)$ and $p_k(x)$, respectively, such that the roots of $p_k$ are interlaced with the roots of $p_{k-1}$, i.e.,

$r_{k,1}<r_{k-1,1}<r_{k,2}<r_{k-1,2}<\dots <r_{k-1,k-1}<r_{k,k}.$

Then, we have that

$p_{k+1}(r_{k,1})=(r_{k,1}-{\mu }_k)\underset{0}{\underbrace{p_k(r_{k,1})}}-{\nu }_k^2{p}_{k-1}(r_{k,1})=-{\nu }_k^2{p}_{k-1}(r_{k,1}),$

$p_{k+1}(r_{k,2})=(r_{k,2}-{\mu }_k)\underset{0}{\underbrace{p_k(r_{k,2})}}-{\nu }_k^2{p}_{k-1}(r_{k,2})=-{\nu }_k^2{p}_{k-1}(r_{k,2}).$

From the induction hypothesis $p_{k-1}(r_{k,1})$ and $p_{k-1}(r_{k,2})$ have different signs since $r_{k,1}<r_{k-1,1}<r_{k,2}$. Then, there exists $r_{k+1,1}\in (r_{k,1},r_{k,2})$ such that $p_{k+1}(r_{k+1,1})=0$.

Proceeding in the same manner for all the intervals $(r_{k,j},r_{k,j+1})$, we obtain that there exists $r_{k+1,j}\in (r_{k,j},r_{k,j+1})$ such that $p_{k+1}(r_{k+1,j})=0$ for $j=1,\dots ,k-1.$

We now consider the smallest and largest roots of $p_{k+1}$ i.e. $r_{k+1,0},r_{k+1,k+1}$

Since for $k=0,1,2,\dots$, $p_k$ is a monic polynomial,

$\underset{x\to \mathrm{\infty }}{\lim }{p}_k(x)=\mathrm{\infty }$

Hence for any $x>r_k$ (any $x$ larger than the largest root of $p_k$)

$p_k(x)>0$

Therefore

$\underset{x\to \mathrm{\infty }}{\lim }{p}_{k+1}(x)>0$ and $p_{k+1}(r_{k,k})=-{\nu }_k^2{p}_{k-1}(r_{k,k})<0,$

implies there exists $r_{k+1,k}\in (r_{k,k},\mathrm{\infty })$ such that $p_{k+1}(r_{k+1,k})=0$.

If $k+1$ is even, then by similar reasoning

$\underset{x\to -\mathrm{\infty }}{\lim }{p}_{k+1}(x)>0$ and $p_{k+1}(r_{k,1})=-{\nu }_k^2{p}_{k-1}(r_{k,1})<0,$ $\Rightarrow$ there exists $r_{k+1,0}\in (-\mathrm{\infty },r_{k,1})$ such that $p_{k+1}(r_{k+1,0})=0$.

If $k+1$ is odd,

$\underset{x\to -\mathrm{\infty }}{\lim }{p}_{k+1}(x)<0$ and $p_{k+1}(r_{k,1})=-{\nu }_k^2{p}_{k-1}(r_{k,1})>0,$ $\Rightarrow$ there exists $r_{k+1,0}\in (-\mathrm{\infty },r_{k,1})$ such that $p_{k+1}(r_{k+1,0})=0$.

Problem 3c

Show that for every $n=0,1,\dots$ the roots of $p_{k+1}(x)$ are the eigenvalues of the symmetric tridiagonal matrix $T=\left(\begin{array}{ccccc}{\mu }_0& {\nu }_1& 0& \cdots & 0\\ {\nu }_1& {\mu }_1& {\nu }_2& \ddots & ⋮\\ 0& {\nu }_2& {\nu }_2& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & {\nu }_n\\ 0& \cdots & 0& {\nu }_n& {\mu }_n\end{array}\right)$

Solution 3c

Let $T_{n+1}=\left(\begin{array}{ccccc}{\mu }_0& {\nu }_1& 0& \cdots & 0\\ {\nu }_1& {\mu }_1& {\nu }_2& \ddots & ⋮\\ 0& {\nu }_2& {\nu }_2& \ddots & 0\\ ⋮& \ddots & \ddots & \ddots & {\nu }_n\\ 0& \cdots & 0& {\nu }_n& {\mu }_n\end{array}\right)$

Then $\det (x{I}_{n+1}-T_{n+1})$ is a monic polynomial in $x$ of degree $n+1$.

Expanding this determinant along the last row, we have

$\det (x{I}_{n+1}-T_{n+1})=(x-{\mu }_n)\det (x{I}_n-T_n)-{\nu }_n^2\det (x{I}_{n-1}-T_{n-1})(1)$

where $\det (x{I}_n-T_n)$ and $\det (x{I}_{n-1}-T_{n-1})$ are monic polynomials of degree $n$ and $n-1$, respectively.

Notice that $\det (x{I}_1-T_1)=x-{\mu }_0\equiv p_1(x)$ and if we let $p_0(x)\equiv 1$, then (1) is equivalent to the three-term recurrence stated in the problem.

Thus, $p_{n+1}(x)\equiv \det (x{I}_{n+1}-T_{n+1})=0$ shows that the roots of $p_{n+1}(x)$ are the eigenvalues of $T_{n+1}$.