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

## Problem 1

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

## Problem 1a

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

## Solution 1a

We wish to show that ${\displaystyle \,\!\|b-Ax(c)\|=\|\nu e_{1}-H_{k}c\|}$

{\displaystyle {\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 ${\displaystyle 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 ${\displaystyle H_{k}\!\,}$ , a ${\displaystyle (k+1)\times k\!\,}$  upper Hessenberg matrix, into QR form.

The cost is on the order of ${\displaystyle 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 ${\displaystyle k\!\,}$  Given's Rotation multiplies to zero out each of the ${\displaystyle k\!\,}$  subdiagonal entries and hence transform ${\displaystyle H_{k}\!\,}$  into upper triangular form ${\displaystyle 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 ${\displaystyle k\!\,}$  Given's Rotations multiplies is

${\displaystyle 4k+4(k-1)+\ldots +4\cdot 1<4k\cdot k=4k^{2}}$

### Back Solving Cost

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

${\displaystyle 1+2+\ldots +k={\frac {k(k+1)}{2}}={\frac {k^{2}}{2}}+{\frac {k}{2}}}$

## Problem 2

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

## Problem 2a

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

## Solution 2a

The composite trapezoidal rule is

${\displaystyle Q_{T,n}(f)={\frac {h}{2}}(f(a)+f(b))+h\sum _{k=1}^{n-1}f(x_{k})\!\,}$

The error is

${\displaystyle I(f)-Q_{T,n}(f)=-{\frac {(b-a)f^{\prime \prime }(\xi )}{12}}h^{2}}$

where ${\displaystyle \xi \in [a,b]\!\,}$ .

### Derivation of Composite Trapezoidal Error

The local error, the error on one interval, is

${\displaystyle -{\frac {h^{3}}{12}}f^{\prime \prime }(\eta )\!\,}$ .

Observe that

${\displaystyle 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

${\displaystyle \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 ${\displaystyle \xi \in [a,b]\!\,}$  such that

${\displaystyle f^{\prime \prime }(\xi )=\sum _{i=1}^{n}{\frac {f^{\prime \prime }(\eta _{i})}{n}}\!\,}$ .

Multiplying both sides of the equation by ${\displaystyle n\!\,}$ ,

${\displaystyle nf^{\prime \prime }(\xi )=\sum _{i=1}^{n}f^{\prime \prime }(\eta _{i})\!\,}$

Using this relationship, we have

{\displaystyle {\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:

 Assume ${\displaystyle f^{n+1}\!\,}$  exists on ${\displaystyle [a,b]\!\,}$  and ${\displaystyle \{x_{0},x_{1},\ldots ,x_{n}|x\in [a,b]\}.\!\,}$  ${\displaystyle P_{n}\!\,}$  interpolates ${\displaystyle f\!\,}$  at ${\displaystyle \{x_{j}\}_{j=0}^{n}\!\,}$ .
Then there is a ${\displaystyle \xi _{x}\!\,}$  (${\displaystyle \xi \!\,}$  is dependent on ${\displaystyle x\!\,}$ ) such that

${\displaystyle f(x)-p_{n}(x)={\frac {(x-x_{0})(x-x_{1})...(x-x_{n})}{(n+1)!}}f^{(n+1)}(\xi _{x})\!\,}$

where ${\displaystyle \xi _{x}\!\,}$  lies in ${\displaystyle (m,M)\!\,}$ ,  ${\displaystyle m=\min\{x_{0},x_{1},\ldots ,x_{n},X\}\!\,}$ ,  ${\displaystyle M=\max\{x_{0},x_{1},\ldots ,x_{n},X\}\!\,}$


Applying the theorem yields,

${\displaystyle f(x)-p_{1}(x)=(x-a)(x-b){\frac {f^{\prime \prime }(\xi _{x})}{2}}}$

Hence,

{\displaystyle {\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 ${\displaystyle a\!\,}$  is the start of the interval, ${\displaystyle x-a\!\,}$  is always positive. Conversely, since ${\displaystyle b\!\,}$  is the end of the interval, ${\displaystyle x-b\!\,}$  is always negative. Hence, ${\displaystyle w(x)<0\!\,}$  is always of one sign. Hence from the mean value theorem of integration, there exists a ${\displaystyle \zeta \in [a,b]\!\,}$  such that

${\displaystyle E(f)={\frac {f^{\prime \prime }(\zeta )}{2}}\int _{a}^{b}(x-a)(x-b)dx\!\,}$

Note that ${\displaystyle f^{\prime \prime }(\zeta )}$  is a constant and does not depend on ${\displaystyle x\!\,}$ .

Integrating ${\displaystyle \int _{a}^{b}(x-a)(x-b)dx\!\,}$ , yields,

{\displaystyle {\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 ${\displaystyle h\!\,}$ ? You may assume here that ${\displaystyle 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 ${\displaystyle n\!\,}$  points in ${\displaystyle [a,b]\!\,}$  is

${\displaystyle E_{n}=I(f)-Q_{T,n}={\frac {-(b-a)h^{2}}{12}}+{\mathcal {O}}(h^{3})\!\,}$

With ${\displaystyle 2n\!\,}$  points, there are twice as many intervals, but the intervals are half as wide. Hence, the error for the composite trapezoidal rule at ${\displaystyle 2n\!\,}$  points in ${\displaystyle [a,b]\!\,}$  is

${\displaystyle 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 ${\displaystyle h^{2}\!\,}$  term by choosing using an appropriate linear combination of ${\displaystyle E_{n}\!\,}$  and ${\displaystyle E_{2n}\!\,}$ . This gives a new error rule with ${\displaystyle h^{3}\!\,}$ error.

{\displaystyle {\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 ${\displaystyle E_{n}\!\,}$  and ${\displaystyle E_{2n}\!\,}$  on the left had side gives

${\displaystyle I(f)-Q_{T,n}-2I(f)-2Q_{T,2n}={\mathcal {O}}(h^{3})\!\,}$

${\displaystyle I(f)=2Q_{T,2n}-Q_{T,n}+{\mathcal {O}}(h^{3})\!\,}$

If we call our new rule ${\displaystyle {\hat {Q}}\!\,}$  we have

${\displaystyle {\hat {Q}}=2Q_{T,2n}-Q_{T,n}\!\,}$  whose error is on the order of ${\displaystyle {\mathcal {O}}(h^{3})\!\,}$

## Problem 3

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

## Problem 3a

 If ${\displaystyle A_{k}={\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}},}$  specify the orthogonal matrix ${\displaystyle 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 ${\displaystyle {\tilde {A}}_{k}\!\,}$  be the ${\displaystyle \sigma _{k}\!\,}$ -shifted matrix of ${\displaystyle A_{k}\!\,}$  i.e.

${\displaystyle {\tilde {A}}_{k}=A_{k}-\sigma _{k}I={\begin{bmatrix}a_{11}-\sigma _{k}&a_{12}\\a_{21}&a_{22}-\sigma _{k}\end{bmatrix}},}$

${\displaystyle Q_{k}^{T}\!\,}$  is an orthogonal 2x2 Given's rotations matrix. ${\displaystyle Q_{k}^{T}\!\,}$ 's entries are chosen such that when ${\displaystyle Q_{k}^{T}\!\,}$  is pre-multiplied against the 2x2 matrix ${\displaystyle {\tilde {A}}_{k}\!\,}$ , ${\displaystyle Q_{k}^{T}\!\,}$  will zero out the (2,1) entry of ${\displaystyle {\tilde {A}}_{k}\!\,}$  and scale ${\displaystyle {\tilde {A}}_{k}\!\,}$ 's remaining three entries i.e.

${\displaystyle Q_{k}^{T}{\tilde {A}}_{k}={\begin{bmatrix}*&*\\0&*\end{bmatrix}}=R_{k}}$

where ${\displaystyle *\!\,}$  denotes calculable scalar values we are not interested in and ${\displaystyle R_{k}\!\,}$  is our desired upper triangular matrix sought by the QR algorithm.

Since ${\displaystyle Q_{k}^{T}\!\,}$  is orthogonal, the above equation implies

{\displaystyle {\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 ${\displaystyle Q_{k}^{T}\!\,}$  is given by

${\displaystyle Q_{k}^{T}={\begin{bmatrix}c&s\\-s&c\end{bmatrix}}}$

where

{\displaystyle {\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 ${\displaystyle Q_{k}^{T}\!\,}$  yields

${\displaystyle Q_{k}={\begin{bmatrix}c&-s\\s&c\end{bmatrix}}}$

## Problem 3b

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

## Solution 3b

From hypothesis and part (a),

{\displaystyle {\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 ${\displaystyle A_{k+1}(2,1)\!\,}$  be the (2,1) entry of ${\displaystyle A_{k}\!\,}$ . Using the above equation, we can find ${\displaystyle A_{k+1}(2,1)\!\,}$  by finding the inner product of the second row of ${\displaystyle R_{k}\!\,}$  and first column ${\displaystyle Q_{k}\!\,}$  and adding the (2,1) entry of ${\displaystyle \sigma _{k}I\!\,}$  i.e.

{\displaystyle {\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 ${\displaystyle r_{22}\!\,}$  so we need to calculate ${\displaystyle R_{k}\!\,}$ .

From hypothesis and the orthogonality of ${\displaystyle Q_{k}\!\,}$ , we have

{\displaystyle {\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 ${\displaystyle R_{k}\!\,}$ , we can find its (2,2) entry ${\displaystyle r_{2}2\!\,}$  by using inner products.

${\displaystyle r_{22}=-s\cdot a_{12}+c\cdot (a_{22}-\sigma _{k})\!\,}$

Now that we have ${\displaystyle r_{22}\!\,}$  we can calculate ${\displaystyle A_{k+1}(2,1)\!\,}$  by using appropriate substitutions

{\displaystyle {\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 ${\displaystyle \delta \!\,}$  is a small number.

Let our shift ${\displaystyle \sigma _{k}=a_{22}\!\,}$ . Then the above equation yields,

{\displaystyle {\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,

{\displaystyle {\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 ${\displaystyle |a_{12}|\leq (a_{11}-a_{22})^{2}\!\,}$

Hence we have shown that ${\displaystyle A_{k+1}(2,1)\!\,}$  is of order ${\displaystyle \delta ^{2}\!\,}$ .

If ${\displaystyle \delta \!\,}$  is small, the QR convergence rate will be quadratic.