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 $\int _{-1}^{1}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},{\mbox{ and }}P_{3}\!\,$ with respect to a weighted inner product

$\left\langle f,g\right\rangle =\int _{-1}^{1}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 $\left\{1,x,x^{2},x^{3}\right\}\!\,$

Zero Order Polynomial

Let $P_{0}(x)=1\!\,$

First Order Polynomial

${\begin{aligned}P_{1}(x)&=x-{\frac {\left\langle 1,x\right\rangle }{\left\langle 1,1\right\rangle }}\cdot 1\\&=x-\underbrace {\frac {\int _{-1}^{1}x^{5}dx}{\int _{-1}^{1}x^{4}dx}} _{\frac {0}{2/5}}\\&=x\end{aligned}}$

Second Order Polynomial

${\begin{aligned}P_{2}(x)&=x^{2}-{\frac {\left\langle 1,x^{2}\right\rangle }{\left\langle 1,1\right\rangle }}\cdot 1-{\frac {\left\langle x,x^{2}\right\rangle }{\left\langle x,x\right\rangle }}\cdot x\\&=x^{2}-{\frac {\int _{-1}^{1}x^{6}dx}{\int _{-1}^{1}x^{4}dx}}-x\cdot \underbrace {\frac {\int _{-1}^{1}x^{7}dx}{\int _{-1}^{1}x^{6}dx}} _{\frac {0}{2/7}}\\&=x^{2}-{\frac {5}{7}}\end{aligned}}$

Third Order Polynomial

${\begin{aligned}P_{3}(x)&=x^{3}-{\frac {\left\langle 1,x^{3}\right\rangle }{\left\langle 1,1\right\rangle }}\cdot 1-{\frac {\left\langle x,x^{3}\right\rangle }{\left\langle x,x\right\rangle }}\cdot x-{\frac {\left\langle x^{2}-{\frac {5}{7}},x^{3}\right\rangle }{\left\langle x^{2}-{\frac {5}{7}},x^{2}-{\frac {5}{7}}\right\rangle }}\cdot (x^{2}-{\frac {5}{7}})\\&=x^{3}-0-{\frac {7}{2}}x\cdot {\frac {\int _{-1}^{1}x^{8}dx}{\int _{-1}^{1}x^{6}dx}}-0\\&=x^{3}-{\frac {7}{9}}x\end{aligned}}$

Find Roots of Polynomials

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

${\begin{array}{|c|c|}{\mbox{Polynomial}}&{\mbox{Roots}}\\\hline 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}}}\\\\\hline \end{array}}$

Formula to Compute Coefficients

In Gaussian quadrature we have

$w_{j}=\int _{-1}^{1}l_{i}(x)x^{4}dx\!\,$, where

$l_{i}(x)=\prod _{j=1,j\neq i}^{n}{\frac {x-x_{j}}{x_{i}-x_{j}}}\!\,$

1 point formula

$l_{1}(x)=1\!\,$

$w_{1}=\int _{-1}^{1}x^{4}dx={\frac {2}{5}}\!\,$

$\int _{-1}^{1}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}=\int _{-1}^{1}{\frac {x-x_{1}}{x_{1}-x_{2}}}x^{4}dx={\frac {2x_{2}}{5(x_{2}-x_{1})}}\!\,$

$w_{2}=\int _{-1}^{1}{\frac {x-x_{1}}{x_{2}-x_{1}}}x^{4}dx={\frac {2x_{1}}{5(x_{1}-x_{2})}}\!\,$

$\int _{-1}^{1}f(x)x^{4}dx\approx w_{1}f(x_{1})+w_{2}f(x_{2})\!\,$

where $x_{1}{\mbox{ 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}=\int _{-1}^{1}{\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})}}\left({\frac {2}{7}}+{\frac {2x_{2}x_{3}}{5}}\right)\!\,$

$w_{2}=\int _{-1}^{1}{\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})}}\left({\frac {2}{7}}+{\frac {2x_{1}x_{3}}{5}}\right)\!\,$

$w_{3}=\int _{-1}^{1}{\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})}}\left({\frac {2}{7}}+{\frac {2x_{1}x_{2}}{5}}\right)\!\,$

$\int _{-1}^{1}f(x)x^{4}dx\approx w_{1}f(x_{1})+w_{2}f(x_{2})+w_{3}f(x_{3})\!\,$

where $x_{1}{\mbox{, }}x_{2}{\mbox{ 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})\quad p'(x_{i})=f'(x_{i})\quad (0\leq i\leq 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{aligned}E&=\int _{-1}^{1}f(x)x^{4}dx-\sum _{i=0}^{n-1}w_{i}f(x_{i})\\&=\int _{-1}^{1}f(x)x^{4}dx-\sum _{i=0}^{n-1}w_{i}p(x_{i})\\&=\int _{-1}^{1}f(x)x^{4}dx-\int _{-1}^{1}p(x)x^{4}dx\\&=\int _{-1}^{1}(f(x)-p(x))x^{4}dx\\&=\int _{-1}^{1}{\frac {f^{(2n)}(\zeta (x))}{(2n)!}}\prod _{i=0}^{n-1}(x-x_{i})^{2}\cdot x^{4}dx\\&={\frac {f^{(2n)}(\xi )}{(2n)!}}\int _{-1}^{1}x^{4}\cdot \underbrace {\prod _{i=0}^{n-1}(x-x_{i})^{2}} _{(p_{n}(x))^{2}}dx\end{aligned}}\!\,$

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{aligned}E&={\frac {f^{(2)}(\xi )}{(2)!}}\int _{-1}^{1}x^{4}\cdot x^{2}dx\\&={\frac {f^{(2)}(\xi )}{7}}\end{aligned}}\!\,$

The error for 2 point quadrature:

${\begin{aligned}E&={\frac {f^{(4)}(\xi )}{(4)!}}\int _{-1}^{1}x^{4}\cdot (x^{2}-{\frac {5}{7}})^{2}dx\end{aligned}}\!\,$

The error for 3 point quadrature:

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

Problem 2

We wish to solve $Ax=b\!\,$ iteratively where $A={\begin{bmatrix}1&1&-1\\1&2&1\\1&1&1\end{bmatrix}}\!\,$ 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{aligned}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{aligned}}\!\,$

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{aligned}Ax&=b\\((D+L)+U)x&=b\\(D+L)x+Ux&=b\\(D+L)x&=b-Ux\\x&=(D+L)^{-1}(b-Ux)\end{aligned}}\!\,$

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)={\begin{bmatrix}0&1&-1\\{\frac {1}{2}}&0&{\frac {1}{2}}\\1&1&0\end{bmatrix}}\!\,$

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={\begin{bmatrix}0&1&-1\\0&-{\frac {1}{2}}&1\\0&-{\frac {1}{2}}&0\end{bmatrix}}\!\,$

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)\quad {\mbox{ for }}k=1,2,\dots ,\!\,$ initialized by $p_{0}(x)=1{\mbox{ and }}p_{1}(x)=x-\mu _{0}\!\,$, where each $\mu _{k}{\mbox{ 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 $\operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{n}(x)\right\}=\operatorname {span} \left\{1,x,x^{2},\dots ,x^{n}\right\}\!\,$

Solution 3a

We prove the claim by by induction.

Base Case

$p_{0}(x)=1\!\,$ is monic with degree zero, and $\operatorname {span} \left\{p_{0}(x)\right\}=\operatorname {span} \left\{1\right\}\!\,$.

$p_{1}(x)=x-\mu _{0}\!\,$ is monic with degree one, and $\operatorname {span} \left\{p_{0}(x),p_{1}(x)\right\}=\operatorname {span} \left\{1,x\right\}\!\,$.

$p_{2}(x)=(x-\mu _{0})(x-\mu _{1})-\nu _{1}^{2}\!\,$ is monic with degree 2, and $\operatorname {span} \left\{p_{0}(x),p_{1}(x),p_{2}(x)\right\}=\operatorname {span} \left\{1,x,x^{2}\right\}\!\,$.

Induction Step

Assume $p_{k-1}(x)\!\,$ is monic with degree $k-1\!\,$, and $\operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{k-1}(x)\right\}=\operatorname {span} \left\{1,x,\dots ,x^{k-1}\right\}\!\,$.

Also, assume $p_{k}(x)\!\,$ is monic with degree $k\!\,$, and $\operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{k}(x)\right\}=\operatorname {span} \left\{1,x,\dots ,x^{k}\right\}\!\,$.

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,

$\operatorname {span} \left\{p_{0}(x),p_{1}(x),\dots ,p_{k+1}(x)\right\}=\operatorname {span} \left\{1,x,\dots ,x^{k+1}\right\}\!\,$

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

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 $\lim _{x\rightarrow \infty }p_{2}(x)>0\!\,$ $\Rightarrow \!\,$ there exists $r_{22}\in (\mu _{0},\infty )\!\,$ such that $p_{2}(r_{22})=0\!\,$.

$p_{2}(\mu _{0})=-\nu _{1}^{2}\!\,$ and $\lim _{x\rightarrow -\infty }p_{2}(x)>0\!\,$ $\Rightarrow \!\,$ there exists $r_{21}\in (-\infty ,\mu _{0})\!\,$ such that $p_{2}(r_{21})=0\!\,$.

In conclusion, $r_{21}<\underbrace {\mu _{0}} _{r_{11}}<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})\underbrace {p_{k}(r_{k,1})} _{0}-\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})\underbrace {p_{k}(r_{k,2})} _{0}-\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,\ldots \!\,$, $p_{k}\!\,$ is a monic polynomial,

$\lim _{x\rightarrow \infty }p_{k}(x)=\infty \!\,$

Hence for any $x>r_{k}\!\,$ (any $x\!\,$ larger than the largest root of $p_{k}\!\,$)

$p_{k}(x)>0\!\,$

Therefore

$\lim _{x\rightarrow \infty }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},\infty )\!\,$ such that $p_{k+1}(r_{k+1,k})=0\!\,$.

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

$\lim _{x\rightarrow -\infty }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 (-\infty ,r_{k,1})\!\,$ such that $p_{k+1}(r_{k+1,0})=0\!\,$.

If $k+1\!\,$ is odd,

$\lim _{x\rightarrow -\infty }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 (-\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={\begin{pmatrix}\mu _{0}&\nu _{1}&0&\cdots &0\\\nu _{1}&\mu _{1}&\nu _{2}&\ddots &\vdots \\0&\nu _{2}&\nu _{2}&\ddots &0\\\vdots &\ddots &\ddots &\ddots &\nu _{n}\\0&\cdots &0&\nu _{n}&\mu _{n}\end{pmatrix}}$

Solution 3c

Let $T_{n+1}={\begin{pmatrix}\mu _{0}&\nu _{1}&0&\cdots &0\\\nu _{1}&\mu _{1}&\nu _{2}&\ddots &\vdots \\0&\nu _{2}&\nu _{2}&\ddots &0\\\vdots &\ddots &\ddots &\ddots &\nu _{n}\\0&\cdots &0&\nu _{n}&\mu _{n}\end{pmatrix}}$

Then $\operatorname {det} (xI_{n+1}-T_{n+1})\!\,$ is a monic polynomial in $x\!\,$ of degree $n+1\!\,$.

Expanding this determinant along the last row, we have

$\operatorname {det} (xI_{n+1}-T_{n+1})=(x-\mu _{n})\operatorname {det} (xI_{n}-T_{n})-\nu _{n}^{2}\operatorname {det} (xI_{n-1}-T_{n-1})\qquad (1)\!\,$

where $\operatorname {det} (xI_{n}-T_{n})\!\,$ and $\operatorname {det} (xI_{n-1}-T_{n-1})\!\,$ are monic polynomials of degree $n\!\,$ and $n-1\!\,$, respectively.

Notice that $\operatorname {det} (xI_{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 \operatorname {det} (xI_{n+1}-T_{n+1})=0\!\,$ shows that the roots of $p_{n+1}(x)\!\,$ are the eigenvalues of $T_{n+1}\!\,$.