# Numerical Methods Qualification Exam Problems and Solutions (University of Maryland)/January 2008
## Problem 1

* Consider the system $Ax=b\,\!$. The GMRES method starts with a point $\,\!x_{0}$ and normalizes the residual $\,\!r_{0}=b-Ax_{0}$ so that $\,\!v_{1}={\frac {r_{0}}{\nu }}$ has 2-norm one. It then constructs orthonormal Krylov bases $\,\!V_{k}=(v_{1}\,v_{2}\,\cdots v_{m})$ satisfying
* $\,\!AV_{k}=V_{k+1}H_{k}$
* where $\,\!H_{k}$ is a $\,\!(k+1)\times k$ upper Hessenberg matrix. One then looks for an approximation to $\displaystyle x$ of the form
* $\,\!x(c)=x_{0}+V_{k}c$
* choosing $\,\!c_{k}$ so that $\,\!\|r(c)\|=\|b-Ax(c)\|$ is minimized, where $\,\!\|\cdot \|$ is the usual Euclidean norm.

## Problem 1a

* Show that $\,\!c_{k}$ minimizes $\|\nu e_{1}-H_{k}c\|$.

## Solution 1a

We wish to show that $\,\!\|b-Ax(c)\|=\|\nu e_{1}-H_{k}c\|$
${\begin{aligned}\|b-Ax(c)\|&=\|b-A(x_{0}+V_{k}c)\|\\&=\|b-Ax_{0}-AV_{k}c\|\\&=\|r_{0}-AV_{k}c\|\\&=\|r_{0}-V_{k+1}H_{k}c\|\\&=\|\nu v_{1}-V_{k+1}H_{k}c\|\\&=\|V_{k+1}\underbrace {(\nu e_{1}-H_{k}c)} _{h_{c}}\|\\&=(V_{k+1}h_{c},V_{k+1}h_{c})^{\frac {1}{2}}\\&=((V_{k+1}h_{c})^{T}V_{k+1}h_{c})^{\frac {1}{2}}\\&=(h_{c}^{T}V_{k+1}^{T}V_{k+1}h_{c})^{\frac {1}{2}}\\&=(h_{c}^{T}h_{c})^{\frac {1}{2}}\\&=\|h_{c}\|\\&=\|\nu e_{1}-H_{k}c\|\end{aligned}}$

## Problem 1b

* Suppose we choose to solve the least squares problem in part a for $c_{k}\!\,$ by the method of orthogonal triangularization (QR). What is the order of the floating point operation for this method. Give reasons.

## Solution 1b

We would like to transform $H_{k}\!\,$, a $(k+1)\times k\!\,$ upper Hessenberg matrix, into QR form.
The cost is on the order of $4k^{2}+{\frac {1}{2}}k^{2}=4.5k^{2}$ from the cost of Given's Rotations and backsolving.

### Given's Rotations Cost

We need $k\!\,$ Given's Rotation multiplies to zero out each of the $k\!\,$ subdiagonal entries and hence transform $H_{k}\!\,$ into upper triangular form $R\!\,$,
Each successive Given's Rotations multiply on an upper Hessenberg matrix requires four fewer multiplies because each previous subdiagonal entry has been zeroed out by a Given's Rotation multiply.
Hence the cost of $k\!\,$ Given's Rotations multiplies is
    $4k+4(k-1)+\ldots +4\cdot 1<4k\cdot k=4k^{2}$

### Back Solving Cost

$R\!\,$ is a $(k+1)\times k\!\,$ upper triangular matrix with last row zero. Hence, we need to backsolve $k\!\,$ upper triangular rows.
    $1+2+\ldots +k={\frac {k(k+1)}{2}}={\frac {k^{2}}{2}}+{\frac {k}{2}}$

## Problem 2

* We wish to approximate the integral $I(f)\equiv \int _{a}^{b}f(x)dx$

## Problem 2a

* State the composite trapezoidal rule $Q_{T,n}\!\,$ for approximating $I(f)\!\,$ on a uniform partitioning of width $h={\frac {b-a}{n}}\!\,$, and give a formula for the error $I(f)-Q_{T,n}\!\,$ that is in a form suitable for extrapolation.

## Solution 2a

The composite trapezoidal rule is
    $Q_{T,n}(f)={\frac {h}{2}}(f(a)+f(b))+h\sum _{k=1}^{n-1}f(x_{k})\!\,$
The error is
    $I(f)-Q_{T,n}(f)=-{\frac {(b-a)f^{\prime \prime }(\xi )}{12}}h^{2}$
where $\xi \in [a,b]\!\,$.

### Derivation of Composite Trapezoidal Error

The local error, the error on one interval, is
    $-{\frac {h^{3}}{12}}f^{\prime \prime }(\eta )\!\,$.
Observe that
    $n\min _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\leq \sum _{i=1}^{n}f^{\prime \prime }(\eta _{i})\leq n\max _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\!\,$
which implies
    $\min _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\leq \sum _{i=1}^{n}{\frac {f^{\prime \prime }(\eta _{i})}{n}}\leq \max _{\eta \in [a,b]}f^{\prime \prime }(\eta _{i})\!\,$
Hence, the Intermediate Value Theorem implies there exists a $\xi \in [a,b]\!\,$ such that
    $f^{\prime \prime }(\xi )=\sum _{i=1}^{n}{\frac {f^{\prime \prime }(\eta _{i})}{n}}\!\,$.
Multiplying both sides of the equation by $n\!\,$,
    $nf^{\prime \prime }(\xi )=\sum _{i=1}^{n}f^{\prime \prime }(\eta _{i})\!\,$
Using this relationship, we have
${\begin{aligned}I(f)-Q_{T,n}(f)&=\sum _{i=1}^{n}I_{i}(f)-Q_{i}(f)\\&=\sum _{i=1}^{n}\left[-{\frac {h^{3}}{12}}f^{\prime \prime }(\eta _{i})\right]\\&=-{\frac {h^{3}}{12}}f^{\prime \prime }(\xi )n\\&=-{\frac {h^{3}}{12}}f^{\prime \prime }(\xi )({\frac {b-a}{h}})\\&=-{\frac {h^{2}}{12}}f^{\prime \prime }(\xi )(b-a)\end{aligned}}$

### Derivation of Local Error

The error in polynomial interpolation can be found by using the following theorem:
Applying the theorem yields,
    $f(x)-p_{1}(x)=(x-a)(x-b){\frac {f^{\prime \prime }(\xi _{x})}{2}}$
Hence,
${\begin{aligned}E(f)&=\int _{a}^{b}f(x)-p_{1}(x)dx\\&=\int _{a}^{b}\underbrace {\underbrace {(x-a)} _{>0}\underbrace {(x-b)} _{<0}} _{w(x)}{\frac {f^{\prime \prime }(\xi _{x})}{2}}dx\end{aligned}}$
Since $a\!\,$ is the start of the interval, $x-a\!\,$ is always positive. Conversely, since $b\!\,$ is the end of the interval, $x-b\!\,$ is always negative. Hence, $w(x)<0\!\,$ is always of one sign. Hence from the mean value theorem of integration, there exists a $\zeta \in [a,b]\!\,$ such that
    $E(f)={\frac {f^{\prime \prime }(\zeta )}{2}}\int _{a}^{b}(x-a)(x-b)dx\!\,$
Note that $f^{\prime \prime }(\zeta )$ is a constant and does not depend on $x\!\,$.
Integrating $\int _{a}^{b}(x-a)(x-b)dx\!\,$, yields,
${\begin{aligned}E(f)&={\frac {f^{\prime \prime }(\zeta )}{2}}\left(-{\frac {(b-a)}{6}}\right)\\&=-{\frac {f^{\prime \prime }(\zeta )(b-a)}{12}}\!\,\end{aligned}}$

## Problem 2b

* Use the error formula to derive a new quadrature rule obtained by performing one step of extrapolation on the composite trapezoidal rule. What is this rule, and how does its error depend on $h\!\,$? You may assume here that $f\!\,$ is as smooth as you need it to be.

## Solution 2b

STOER AND BUESCH pg 162
INTRODUCTION TO APPLIED NUMERICAL ANALYSIS by RICHARD HAMMING pg 178
The error for the composite trapezoidal rule at $n\!\,$ points in $[a,b]\!\,$ is
$E_{n}=I(f)-Q_{T,n}={\frac {-(b-a)h^{2}}{12}}+{\mathcal {O}}(h^{3})\!\,$
With $2n\!\,$ points, there are twice as many intervals, but the intervals are half as wide. Hence, the error for the composite trapezoidal rule at $2n\!\,$ points in $[a,b]\!\,$ is
$E_{2n}=I(f)-Q_{T,2n}=2\cdot {\frac {-(b-a)({\frac {h}{2}})^{2}}{12}}+{\mathcal {O}}(h^{3})={\frac {-(b-a)h^{2}}{24}}+{\mathcal {O}}(h^{3})\!\,$
We can eliminate the $h^{2}\!\,$ term by choosing using an appropriate linear combination of $E_{n}\!\,$ and $E_{2n}\!\,$. This gives a new error rule with $h^{3}\!\,$error.
${\begin{aligned}E_{n}-2E_{2n}&={\frac {-(b-a)h^{2}}{12}}-2\cdot {\frac {-(b-a)h^{2}}{24}}+{\mathcal {O}}(h^{3})\\&={\mathcal {O}}(h^{3})\end{aligned}}$
Substituting our equations for $E_{n}\!\,$ and $E_{2n}\!\,$ on the left had side gives
$I(f)-Q_{T,n}-2I(f)-2Q_{T,2n}={\mathcal {O}}(h^{3})\!\,$
$I(f)=2Q_{T,2n}-Q_{T,n}+{\mathcal {O}}(h^{3})\!\,$
If we call our new rule ${\hat {Q}}\!\,$ we have
${\hat {Q}}=2Q_{T,2n}-Q_{T,n}\!\,$ whose error is on the order of ${\mathcal {O}}(h^{3})\!\,$

## Problem 3

* Consider the shifted QR iteration for computing a 2 x 2 matrix $A\!\,$: starting with $A_{0}=A\!\,$, compute
* $A_{k}-\sigma _{k}I=Q_{k}R_{k}\quad \quad \quad A_{k+1}=R_{k}Q_{k}+\sigma _{k}I$
* where $\sigma _{k}\!\,$ is a scalar shift.

## Problem 3a

* If
* $A_{k}={\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}},$
* specify the orthogonal matrix $Q_{k}\!\,$ used to perform this step when Givens rotations are used. The matrix should be described in terms of the entries of the shifted matrix.

## Solution 3a

Let ${\tilde {A}}_{k}\!\,$ be the $\sigma _{k}\!\,$-shifted matrix of $A_{k}\!\,$ i.e.
    ${\tilde {A}}_{k}=A_{k}-\sigma _{k}I={\begin{bmatrix}a_{11}-\sigma _{k}&a_{12}\\a_{21}&a_{22}-\sigma _{k}\end{bmatrix}},$
$Q_{k}^{T}\!\,$ is an orthogonal 2x2 Given's rotations matrix. $Q_{k}^{T}\!\,$'s entries are chosen such that when $Q_{k}^{T}\!\,$ is pre-multiplied against the 2x2 matrix ${\tilde {A}}_{k}\!\,$, $Q_{k}^{T}\!\,$ will zero out the (2,1) entry of ${\tilde {A}}_{k}\!\,$ and scale ${\tilde {A}}_{k}\!\,$'s remaining three entries i.e.
    $Q_{k}^{T}{\tilde {A}}_{k}={\begin{bmatrix}*&*\\0&*\end{bmatrix}}=R_{k}$
where $*\!\,$ denotes calculable scalar values we are not interested in and $R_{k}\!\,$ is our desired upper triangular matrix sought by the QR algorithm.
Since $Q_{k}^{T}\!\,$ is orthogonal, the above equation implies
${\begin{aligned}Q_{k}^{T}{\tilde {A}}_{k}&=R_{k}\\Q_{k}Q_{k}^{T}{\tilde {A}}_{k}&=Q_{k}R_{k}\\{\tilde {A}}_{k}&=Q_{k}R_{k}\end{aligned}}$
The Given's rotation $Q_{k}^{T}\!\,$ is given by
    $Q_{k}^{T}={\begin{bmatrix}c&s\\-s&c\end{bmatrix}}$
where
${\begin{aligned}c&={\frac {a_{11}-\sigma _{k}}{r}}\\s&={\frac {a_{21}}{r}}\\r&={\sqrt {(a_{11}-\sigma _{k})^{2}+a_{21}^{2}}}\end{aligned}}$
Taking the transpose of $Q_{k}^{T}\!\,$ yields
    $Q_{k}={\begin{bmatrix}c&-s\\s&c\end{bmatrix}}$

## Problem 3b

* Suppose $a_{21}=\delta \!\,$, a small number, and $|a_{12}|\leq (a_{11}-a_{22})^{2}\!\,$. Demonstrate that with an appropriate shift $\sigma _{k}\!\,$, the (2,1)-entry of $A_{k+1}\!\,$ is of magnitude $O(\delta ^{2})\!\,$. What does this suggest about the convergence rate of the QR iteration?

## Solution 3b

From hypothesis and part (a),
${\begin{aligned}A_{k+1}&=R_{k}Q_{k}+\sigma _{k}I\\&={\begin{bmatrix}*&*\\0&r_{22}\end{bmatrix}}{\begin{bmatrix}c&-s\\s&c\end{bmatrix}}+{\begin{bmatrix}\sigma _{k}&0\\0&\sigma _{k}\end{bmatrix}}\end{aligned}}$
Let $A_{k+1}(2,1)\!\,$ be the (2,1) entry of $A_{k}\!\,$. Using the above equation, we can find $A_{k+1}(2,1)\!\,$ by finding the inner product of the second row of $R_{k}\!\,$ and first column $Q_{k}\!\,$ and adding the (2,1) entry of $\sigma _{k}I\!\,$ i.e.
${\begin{aligned}A_{k+1}(2,1)&=(0\cdot c+r_{22}\cdot s)+0\\&=r_{22}s\end{aligned}}$
We need to find the value of $r_{22}\!\,$ so we need to calculate $R_{k}\!\,$.
From hypothesis and the orthogonality of $Q_{k}\!\,$, we have
${\begin{aligned}A_{k}-\sigma _{k}I&=Q_{k}R_{k}\\Q_{k}^{T}(A_{k}-\sigma _{k}I)&=Q_{k}^{T}Q_{k}R_{k}\\Q_{k}^{T}(A_{k}-\sigma _{k}I)&=R_{k}\\R_{k}&=Q_{k}^{T}(A_{k}-\sigma _{k}I)\\&={\begin{bmatrix}c&s\\-s&c\end{bmatrix}}{\begin{bmatrix}a_{11}-\sigma _{k}&a_{12}\\\delta &a_{22}-\sigma _{k}\end{bmatrix}}\end{aligned}}$
From $R_{k}\!\,$, we can find its (2,2) entry $r_{2}2\!\,$ by using inner products.
    $r_{22}=-s\cdot a_{12}+c\cdot (a_{22}-\sigma _{k})\!\,$
Now that we have $r_{22}\!\,$ we can calculate $A_{k+1}(2,1)\!\,$ by using appropriate substitutions
${\begin{aligned}A_{k+1}(2,1)&=r_{22}\cdot s\\&=(-sa_{12}+c(a_{22}-\sigma _{k}))s\\&=-s^{2}a_{12}+cs(a_{22}-\sigma _{k})\\&={\frac {-\delta ^{2}a_{12}}{(a_{11}-\sigma _{k})^{2}+\delta ^{2}}}+{\frac {(a_{11}-\sigma _{k})(a_{22}-\sigma _{k})\delta }{(a_{11}-\sigma _{k})^{2}+\delta ^{2}}}\\&={\frac {-\delta ^{2}a_{12}+\delta (a_{11}-\sigma _{k})(a_{22}-\sigma _{k})}{(a_{11}-\sigma _{k})^{2}+\delta ^{2}}}\\&\approx {\frac {-\delta ^{2}a_{12}+\delta (a_{11}-\sigma _{k})(a_{22}-\sigma _{k})}{(a_{11}-\sigma _{k})^{2}}}\end{aligned}}$
since $\delta \!\,$ is a small number.
Let our shift $\sigma _{k}=a_{22}\!\,$. Then the above equation yields,
${\begin{aligned}A_{k+1}(2,1)&\approx {\frac {-\delta ^{2}a_{12}+\delta (a_{11}-\sigma _{k})(a_{22}-\sigma _{k})}{(a_{11}-\sigma _{k})^{2}}}\\&\approx {\frac {-\delta ^{2}a_{12}+\delta (a_{11}-a_{22})(a_{22}-a_{22})}{(a_{11}-a_{22})^{2}}}\\&={\frac {-\delta ^{2}a_{12}}{(a_{11}-a_{22})^{2}}}\end{aligned}}$
Hence,
${\begin{aligned}|A_{k+1}(2,1)|&\approx \left|{\frac {\delta ^{2}a_{12}}{(a_{11}-a_{22})^{2}}}\right|\\&\leq |\delta ^{2}|\end{aligned}}$
since $|a_{12}|\leq (a_{11}-a_{22})^{2}\!\,$
Hence we have shown that $A_{k+1}(2,1)\!\,$ is of order $\delta ^{2}\!\,$.
If $\delta \!\,$ is small, the QR convergence rate will be quadratic.