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

## Problem 1

 Let ${\displaystyle x_{1},x_{2},\ldots ,x_{n}}$  be distinct real numbers, and consider the matrix ${\displaystyle V=\left({\begin{array}{ccccc}1&x_{1}&x_{1}^{2}&\cdots &x_{1}^{n-1}\\1&x_{2}&x_{2}^{2}&\cdots &x_{2}^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&x_{n}&x_{n}^{2}&\cdots &x_{n}^{n-1}\\\end{array}}\right)}$ Prove that ${\displaystyle V^{-1}}$  can be expressed as ${\displaystyle V^{-1}=\left({\begin{array}{c}\\\\\alpha ^{(1)},\alpha ^{(2)},\alpha ^{(3)},\cdots ,\alpha ^{(n)}\\\\\\\end{array}}\right)}$ where the column vector ${\displaystyle \alpha ^{(j)}=(\alpha _{1}^{(j)},\alpha _{2}^{(j)},\cdots ,\alpha _{n}^{(j)})^{T}}$  generates a polynomial ${\displaystyle p_{j}(x)=\sum _{i=1}^{n}\alpha _{i}^{(j)}x^{i-1}}$ satisfying ${\displaystyle p_{j}(x)={\frac {\prod _{i=1,i\neq j}^{n}(x-x_{i})}{\prod _{i=1,i\neq j}^{n}(x_{j}-x_{i})}}}$ You may cite anything you know about polynomial interpolation to justify that ${\displaystyle V}$  is nonsingular.

## Solution 1

We want to show that

${\displaystyle \,\!VV^{-1}=I}$

or equivalently the ${\displaystyle \,\!i}$ th, ${\displaystyle \,\!j}$ th entry of ${\displaystyle \,\!VV^{-1}}$  is ${\displaystyle \,\!1}$  for ${\displaystyle \,\!i=j}$  and ${\displaystyle \,\!0}$  for ${\displaystyle \,\!i\neq j}$  i.e.

${\displaystyle \{VV^{-1}\}_{ij}=\delta _{ij}=\left\{{\begin{array}{cc}1&{\mbox{ if }}i=j\\0&{\mbox{ if }}i\neq j\end{array}}\right.}$

First notice that

${\displaystyle \,\!\{VV^{-1}\}_{ij}={\begin{bmatrix}1+x_{i}^{1}+x_{i}^{2}+\cdots +x_{i}^{n-1}\end{bmatrix}}{\begin{bmatrix}\alpha _{1}^{(j)}\\\alpha _{2}^{(j)}\\\vdots \\\alpha _{n}^{(j)}\end{bmatrix}}=\sum _{k=1}^{n}\alpha _{k}^{(j)}x_{i}^{k-1}}$

Also notice that

${\displaystyle \,\!p_{j}(x_{i})=\sum _{k=1}^{n}\alpha _{k}^{(j)}x_{i}^{k-1}=\delta _{ij}}$

Hence,

${\displaystyle \,\!\{VV^{-1}\}_{ij}=p_{j}(x_{i})=\delta _{ij}}$

## Problem 2

 Let ${\displaystyle A\,\!}$  be a real matrix of order ${\displaystyle n\,\!}$  with eigenvalues ${\displaystyle \lambda _{1},\lambda _{2},\ldots ,\lambda _{n}\,\!}$  and ${\displaystyle n\,\!}$  linearly independent eigenvectors ${\displaystyle v_{1},v_{2},\ldots ,v_{n}\,\!}$ . Suppose you want to find an eigenvector ${\displaystyle v_{j}\,\!}$  corresponding to the eigenvalue ${\displaystyle \lambda _{j}\,\!}$  and you are given ${\displaystyle \sigma \,\!}$  such that ${\displaystyle |\lambda _{j}-\sigma |<|\lambda _{i}-\sigma |{\mbox{ for all }}i\neq j\,\!}$ Specify a numerical algorithm for finding ${\displaystyle v_{j}\,\!}$  and give a convergence proof for this algorithm demonstrating that it is convergent under appropriate circumstances

## Solution 2

### Shifted Inverse Power Method

Let

{\displaystyle {\begin{aligned}{\tilde {A}}&=(A-\sigma I)^{-1}\\\end{aligned}}}

Then,

${\displaystyle {\tilde {\lambda _{i}}}={\frac {1}{|\lambda _{i}-\sigma |}}}$

which implies

${\displaystyle {\tilde {\lambda _{j}}}>{\tilde {\lambda }}_{i}{\mbox{ for all }}i\neq j}$

Since shifting the eigenvalues and inverting the matrix does not affect the eigenvectors,

${\displaystyle {\tilde {A}}v_{i}={\tilde {\lambda _{i}}}v_{i}\quad \quad i=1,2,\ldots ,n}$

Assume ${\displaystyle \|v_{i}\|=1}$  for all ${\displaystyle \!\,i}$ . Generate ${\displaystyle w_{0},w_{1},w_{2},\ldots }$  to find ${\displaystyle \,\!v_{j}}$ . Start with arbitrary ${\displaystyle \,\!w_{0}}$  such that ${\displaystyle \,\!\|w_{0}\|=1}$ .

For ${\displaystyle \,\!k=0,1,2,\ldots }$
${\displaystyle \,\!{\hat {w}}_{k+1}={\tilde {A}}w_{k}}$
${\displaystyle \,\!w_{k+1}={\frac {{\hat {w}}_{k+1}}{\|{\hat {w}}_{k+1}\|}}}$
${\displaystyle \lambda _{j}^{(k+1)}={\frac {({\hat {w}}_{k+1},w_{k})}{(w_{k},w_{k})}}={\frac {({\tilde {A}}w_{k},w_{k})}{(w_{k},w_{k})}}}$  (Rayleigh quotient)
End


### Convergence of Power Method

If ${\displaystyle \,\!|\lambda _{1}|>|\lambda _{i}|}$ , for all ${\displaystyle \!\,i\neq 1}$ , then ${\displaystyle \,\!A^{k}w_{0}}$  will be dominated by ${\displaystyle \,\!v_{1}}$ .

Since ${\displaystyle v_{1},v_{2},\ldots ,v_{n}}$  are linearly independent, they form a basis of ${\displaystyle R^{n}\,\!}$ . Hence,

${\displaystyle w_{0}=\alpha _{1}v_{1}+\alpha _{2}v_{2}+\ldots +\alpha _{n}v_{n}\,\!}$

From the definition of eigenvectors,

{\displaystyle {\begin{aligned}Aw_{0}&=\alpha _{1}\lambda _{1}v_{1}+\alpha _{2}\lambda _{2}v_{2}+\ldots +\alpha _{n}\lambda _{n}v_{n}\\A^{2}w_{0}&=\alpha _{1}\lambda _{1}^{2}v_{1}+\alpha _{2}\lambda _{2}^{2}v_{2}+\ldots +\alpha _{n}\lambda _{n}^{2}v_{n}\\&\vdots \\A^{k}w_{0}&=\alpha _{1}\lambda _{1}^{k}v_{1}+\alpha _{2}\lambda _{k}^{2}v_{2}+\ldots +\alpha _{n}\lambda _{n}^{k}v_{n}\\\end{aligned}}}

To find a general form of ${\displaystyle w_{k}\,\!}$ , the approximate eigenvector at the kth step, examine a few steps of the algorithm:

{\displaystyle {\begin{aligned}{\hat {w}}_{1}&=Aw_{0}\\w_{1}&={\frac {{\hat {w}}_{1}}{\|{\hat {w}}_{1}\|}}={\frac {Aw_{0}}{\|Aw_{0}\|}}\\{\hat {w}}_{2}&=Aw_{1}={\frac {A^{2}w_{0}}{\|Aw_{0}\|}}\\w_{2}&={\frac {{\hat {w}}_{2}}{\|{\hat {w}}_{2}\|}}={\frac {\frac {A^{2}w_{0}}{\|Aw_{0}\|}}{\frac {\|A^{2}w_{0}\|}{\|Aw_{0}\|}}}={\frac {A^{2}w_{0}}{\|A^{2}w_{0}\|}}\\{\hat {w}}_{3}&={\frac {A^{3}w_{0}}{\|A^{2}w_{0}\|}}\\w_{3}&={\frac {\frac {A^{3}w_{0}}{\|A^{2}w_{0}\|}}{\frac {\|A^{3}w_{0}\|}{\|A^{2}w_{0}\|}}}={\frac {A^{3}w_{0}}{\|A^{3}w_{0}\|}}\\\end{aligned}}}

From induction,

${\displaystyle w_{k}={\frac {A^{k}w_{0}}{\|A^{k}w_{0}\|}}}$

Hence,

{\displaystyle {\begin{aligned}w_{k}&={\frac {A^{k}w_{0}}{\|A^{k}w_{0}\|}}\\&={\frac {\alpha _{1}\lambda _{1}^{k}v_{1}+\alpha _{2}\lambda _{2}^{k}v_{2}+\ldots +\alpha _{n}\lambda _{n}^{k}v_{n}}{\|A^{k}w_{0}\|}}\\&={\frac {\lambda _{1}^{k}(\alpha _{1}v_{1}+\alpha _{2}({\frac {\lambda _{2}}{\lambda _{1}}})^{k}v_{2}+\ldots +\alpha _{n}({\frac {\lambda _{n}}{\lambda _{1}}})^{k}v_{n})}{\|A^{k}w_{0}\|}}\end{aligned}}}

Comparing a weighted ${\displaystyle w_{k}\!\,}$  and ${\displaystyle v_{1}\!\,}$ ,

{\displaystyle {\begin{aligned}\left\|{\frac {\|A^{k}w_{0}\|}{\lambda _{1}^{k}}}w_{k}-\alpha _{1}v_{1}\right\|&=\left\|\alpha _{2}\left({\frac {\lambda _{2}}{\lambda _{1}}}\right)^{k}v_{2}+\ldots +\alpha _{n}\left({\frac {\lambda _{n}}{\lambda _{1}}}\right)^{k}v_{n}\right\|\\&\leq \alpha _{2}\left|{\frac {\lambda _{2}}{\lambda _{1}}}\right|^{k}+\ldots +\alpha _{n}\left|{\frac {\lambda _{n}}{\lambda _{1}}}\right|^{k}\end{aligned}}}

since ${\displaystyle \|v_{i}\|=1}$  by assumption.

The above expression goes to ${\displaystyle \,\!0}$  as ${\displaystyle k\rightarrow \infty }$  since ${\displaystyle \,\!|\lambda _{1}|>|\lambda _{i}|}$ , for all ${\displaystyle \!\,i\neq 1}$ . Hence as ${\displaystyle \!\,k}$  grows, ${\displaystyle \!\,w_{k}}$  is parallel to ${\displaystyle \!\,v_{1}}$ . Because ${\displaystyle \|w_{k}\|=1}$ , it must be that ${\displaystyle w_{k}\rightarrow \pm v_{1}}$ .

## Problem 3

 Suppose A is an upper triangular, nonsingular matrix. Show that both Jacobi iterations and Gauss-Seidel iterations converge in finitely many steps when used to solve ${\displaystyle A\mathbf {x} =\mathbf {b} }$ .

## Solution 3

### Derivation of Iterations

Let ${\displaystyle A=D+L+U\!\,}$  where ${\displaystyle D\!\,}$  is a diagonal matrix, ${\displaystyle L\!,}$  is a lower triangular matrix with a zero diagonal and ${\displaystyle U\!\,}$  is an upper triangular matrix also with a zero diagonal.

The Jacobi iteration can be found by substituting into ${\displaystyle Ax=b\!\,}$ , grouping ${\displaystyle L+U\!\,}$ , and solving for ${\displaystyle x\!\,}$  i.e.

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

Since ${\displaystyle L=0\!\,}$  by hypothesis, the iteration is

${\displaystyle x^{(i+1)}=D^{-1}(b-Ux^{(i)})\!\,}$

Similarly, the Gauss-Seidel iteration can be found by substituting into ${\displaystyle Ax=b\!\,}$ , grouping ${\displaystyle D+L\!\,}$ , and solving for ${\displaystyle x\!\,}$  i.e.

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

Since ${\displaystyle L=0\!\,}$  by hypothesis, the iteration has identical from as Jacopi:

${\displaystyle x^{(i+1)}=D^{-1}(b-Ux^{(i)})\!\,}$

### Convergence in Finite Number of Steps

Jacobi and Gauss-Seidel are iterative methods that split the matrix ${\displaystyle A\in \mathbb {R} ^{n\times n}\!\,}$  into ${\displaystyle D\!\,}$ , ${\displaystyle U\!\,}$ , and ${\displaystyle L\!\,}$ : Diagonal, Upper (everything above the diagonal), and Lower (everything below the diagonal) triangular matrices, respectively. Their iterations go as such

${\displaystyle x^{(i+1)}=(D+L)^{-1}(b-Ux^{(i)})\!\,}$  (Gauss-Seidel)
${\displaystyle x^{(i+1)}=D^{-1}(b-(U+L)x^{(i)})\!\,}$  (Jacobi)

In our case ${\displaystyle A\!\,}$  is upper triangular, so ${\displaystyle L\!\,}$  is the zero matrix. As a result, the Gauss-Seidel and Jacobi methods take on the following identical form

${\displaystyle x^{(i+1)}=D^{-1}(b-Ux^{(i)})\!\,}$

Additionally, ${\displaystyle x\!\,}$  can be written

${\displaystyle x=D^{-1}(b-Ux)\!\,}$

Subtracting ${\displaystyle x\!\,}$  from ${\displaystyle x^{(i+1)}\!\,}$  we get

{\displaystyle {\begin{aligned}e^{(i+1)}&=x^{(i+1)}-x\\&=D^{-1}(b-Ux^{(i)})-D^{-1}(b-Ux)\\&=D^{-1}U(-x^{(i)}+x)\\&=D^{-1}Ue^{(i)}\\&=D^{-1}U(D^{-1}Ue^{(i-1)})=(D^{-1}U)^{2}e^{(i-1)}\\&\vdots \\&=(D^{-1}U)^{i+1}e^{(0)}\end{aligned}}}

In our problem, ${\displaystyle D^{-1}\!\,}$  is diagonal and ${\displaystyle U\!\,}$  is upper triangular with zeros along the diagonal. Notice that the product ${\displaystyle D^{-1}U\!\,}$  will also be upper triangular with zeros along the diagonal.

Let ${\displaystyle R=D^{-1}U\!\,}$ , ${\displaystyle R\in \mathbb {R} ^{n\times n}\!\,}$

${\displaystyle R={\begin{pmatrix}0&*&*&\cdots &*\\0&0&*&\cdots &*\\\vdots &\vdots &\ddots &\ddots &\vdots \\0&\cdots &\cdots &0&*\\0&\cdots &\cdots &0&0\end{pmatrix}}}$

Also, let ${\displaystyle {\tilde {R}}\in \mathbb {R} ^{n\times n}\!\,}$  be the related matrix

{\displaystyle {\begin{aligned}{\tilde {R}}&=\left({\begin{array}{ccc|cccc}0&\cdots &0&*&*&\cdots &*\\0&\cdots &0&0&*&\cdots &*\\\vdots &\vdots &\vdots &\vdots &\ddots &\ddots &\vdots \\0&\cdots &0&0&\cdots &0&*\\\hline 0&\cdots &0&0&\cdots &0&0\\\vdots &\vdots &\vdots &\vdots &\cdots &\vdots &\vdots \\0&\cdots &0&0&\cdots &0&0\\\end{array}}\right)\\&=\left({\begin{array}{c|c}\mathbf {a} &T\\\hline \mathbf {0} &\mathbf {a} ^{T}\\\end{array}}\right)\end{aligned}}}

Where ${\displaystyle \mathbf {a} \!\,}$  is ${\displaystyle k\times (n-k)\!\,}$ , ${\displaystyle T\!\,}$  is ${\displaystyle k\times k\!\,}$ , and ${\displaystyle \mathbf {0} \!\,}$  is ${\displaystyle (n-k)\times (n-k)\!\,}$ .

Finally, the product ${\displaystyle R{\tilde {R}}\!\,}$  (call it (1)) is

{\displaystyle {\begin{aligned}R{\tilde {R}}&={\begin{pmatrix}0&*&*&\cdots &*\\0&0&*&\cdots &*\\\vdots &\vdots &\ddots &\ddots &\vdots \\0&\cdots &\cdots &0&*\\0&\cdots &\cdots &0&0\end{pmatrix}}\left({\begin{array}{ccc|cccc}0&\cdots &0&*&*&\cdots &*\\0&\cdots &0&0&*&\cdots &*\\\vdots &\vdots &\vdots &\vdots &\ddots &\ddots &\vdots \\0&\cdots &0&0&\cdots &0&*\\\hline 0&\cdots &0&0&\cdots &0&0\\\vdots &\vdots &\vdots &\vdots &\cdots &\vdots &\vdots \\0&\cdots &0&0&\cdots &0&0\\\end{array}}\right)\\&=\left({\begin{array}{ccc|ccccc}0&\cdots &0&0&*&\cdots &\cdots &*\\0&\cdots &0&0&0&*&\cdots &*\\\vdots &\vdots &\vdots &\vdots &\cdots &\cdots &\ddots &\vdots \\0&\cdots &0&0&\cdots &\cdots &0&*\\\hline 0&\cdots &0&0&\cdots &\cdots &0&0\\\vdots &\vdots &\vdots &\vdots &\cdots &\cdots &\vdots &\vdots \\0&\cdots &0&0&\cdots &\cdots &0&0\\\end{array}}\right)\\&=\left({\begin{array}{c|c}\mathbf {a} &{\tilde {T}}\\\hline \mathbf {0} &\mathbf {a} ^{T}\\\end{array}}\right)\end{aligned}}}

Here ${\displaystyle {\tilde {T}}\!\,}$  is almost identical in structure to ${\displaystyle T\!\,}$ , except that its diagonal elements are zeros.

At this point the convergence in ${\displaystyle n\!\,}$  steps (the size of the starting matrix) should be apparent since ${\displaystyle R\!\,}$  is just ${\displaystyle {\tilde {R}}\!\,}$  where ${\displaystyle k=0\!\,}$  and each time ${\displaystyle R\!\,}$  is multiplied by ${\displaystyle {\tilde {R}}\!\,}$ , the k-th super-diagonal is zeroed out (where k=0 is the diagonal itself). After ${\displaystyle n-1\!\,}$  applications of ${\displaystyle R\!\,}$ , the result will be the zero matrix of size ${\displaystyle n\!\,}$ .

In brief, ${\displaystyle R^{n}\!\,}$  is the zero matrix of size ${\displaystyle n\!\,}$ .

Therefore ${\displaystyle e^{(n)}=(D^{-1}U)^{n}e^{(0)}\equiv 0\!\,}$ , i.e. the Jacobi and Gauss-Seidel Methods used to solve ${\displaystyle A\mathbf {x} =b\!\,}$  converge in ${\displaystyle n\!\,}$  steps when ${\displaystyle A\in \mathbb {R} ^{n\times n}\!\,}$  is upper triangular.

### Examples of (1)

${\displaystyle n=3\!\,}$

${\displaystyle R={\begin{pmatrix}0&*&*\\0&0&*\\0&0&0\\\end{pmatrix}}}$

${\displaystyle R^{2}={\begin{pmatrix}0&*&*\\0&0&*\\0&0&0\\\end{pmatrix}}{\begin{pmatrix}0&*&*\\0&0&*\\0&0&0\\\end{pmatrix}}={\begin{pmatrix}0&0&*\\0&0&0\\0&0&0\\\end{pmatrix}}}$

${\displaystyle R^{3}={\begin{pmatrix}0&*&*\\0&0&*\\0&0&0\\\end{pmatrix}}{\begin{pmatrix}0&0&*\\0&0&0\\0&0&0\\\end{pmatrix}}={\begin{pmatrix}0&0&0\\0&0&0\\0&0&0\\\end{pmatrix}}}$

${\displaystyle n=4\!\,}$

${\displaystyle R={\begin{pmatrix}0&*&*&*\\0&0&*&*\\0&0&0&*\\0&0&0&0\\\end{pmatrix}}}$

${\displaystyle R^{2}={\begin{pmatrix}0&*&*&*\\0&0&*&*\\0&0&0&*\\0&0&0&0\\\end{pmatrix}}{\begin{pmatrix}0&*&*&*\\0&0&*&*\\0&0&0&*\\0&0&0&0\\\end{pmatrix}}={\begin{pmatrix}0&0&*&*\\0&0&0&*\\0&0&0&0\\0&0&0&0\\\end{pmatrix}}}$

${\displaystyle R^{3}={\begin{pmatrix}0&*&*&*\\0&0&*&*\\0&0&0&*\\0&0&0&0\\\end{pmatrix}}{\begin{pmatrix}0&0&*&*\\0&0&0&*\\0&0&0&0\\0&0&0&0\\\end{pmatrix}}={\begin{pmatrix}0&0&0&*\\0&0&0&0\\0&0&0&0\\0&0&0&0\\\end{pmatrix}}}$

${\displaystyle R^{4}={\begin{pmatrix}0&*&*&*\\0&0&*&*\\0&0&0&*\\0&0&0&0\\\end{pmatrix}}{\begin{pmatrix}0&0&0&*\\0&0&0&0\\0&0&0&0\\0&0&0&0\\\end{pmatrix}}={\begin{pmatrix}0&0&0&0\\0&0&0&0\\0&0&0&0\\0&0&0&0\\\end{pmatrix}}}$