Preliminary Examination, Numerical Analysis, August 2007

**Matrix Factorizations** **a)** Prove any two of the following statements: **(i) Schur Decomposition:** Any matrix \(A\in\mathbb{C}^{m\times m}\) can be factored as \(A=Q^*TQ\), where \(Q\) is unitary and \(T\) is upper triangular. **(ii) Singular Value Decomposition:** Any matrix \(A\in\mathbb{C}^{m\times n}\) can be factored as \(A=U\Sigma V^*\), where \(U\in\mathbb{C}^{m\times m}\) and \(V\in\mathbb{C}^{n\times n}\) are unitary and \(\Sigma\in\mathbb{R}^{m\times n}\) is a rectangular matrix whose only nonzero entries are non-negative entries on its diagonal. **(iii) QR Factorization:** Any full-rank matrix \(A\in\mathbb{R}^{m\times n}\) for \(m\geq n\) can be factored \(A=QR\) where \(Q\in\mathbb{R}^{m\times m}\) is orthogonal and \(R\in\mathbb{R}^{m\times n}\) is upper triangular with positive diagonal entries. **b)** Discuss situations in which each of these factorizations is useful in numerical analysis and explain why the factorizations are useful in those situations.

**Interpolation and Integration** **a)** Consider equally spaced points \(x_j=a+jh\), \(j=0,\ldots,J+1\) on the interval \([a,b]\), where \((J+1)h=b-a\). Let \(f(x)\) be a function defined on \([a,b]\). Consider the problem of finding a cubic spline approximation \(s(x)\) to \(f(x)\) that interpolates \(f\) at the points \(x_j\), is twice continuously differentiable, and satisfies \(s''(a)=s''(b)=0\). Does this problem always have a solution? If your answer is yes, derive formulas by which to determine the spline. If your answer is no, explain your reasoning. **b)** Let \(I_n(f)\) denote the result of using the composite Trapezoidal rule to approximate \[ I(f)\equiv \int_a^b f(x)\,dx \] using \(n\) equally sized subintervals of length \(h=(b-a)/n\). It can be shown that the integration error \(E_n(f)\equiv I(f)-I_n(f)\) satisfies \[ E_n(f)=d_2h^2+d_4h^4+d_6h^6+\cdots \] where \(d_2,d_4,d_6,\ldots\) are numbers that depend only on the values of \(f\) and its derivatives at \(a\) and \(b\). Suppose you have a black-box program that, given \(f\), \(a\), \(b\), and \(n\), calculates \(I_n(f)\). Show how to use this program to obtain an \(O(h^4)\) approximation and an \(O(h^6)\) approximation to \(I(f)\).

**Eigenvalue Problems** **a)** Let \(A\) be a real \(n\times n\) matrix with simple eigenvalue \(\lambda\), right eigenvector \(x\) (i.e., \(Ax=\lambda x\)), and left eigenvector \(y\) (i.e., \(y^*A=\lambda y^*\)). Assume that \(\|x\|_2=\|y\|_2=1\). Consider the perturbed matrix \(A+\delta A\), where \(\|\delta A\|_2/\|A\|_2=\epsilon\ll1\). Show that the matrix \(A+\delta A\) has an eigenvalue \(\mu\) such that \(|\lambda-\mu|\) is approximately \(\epsilon/|y^*x|\). (Note that this implies that \(|\lambda-\mu|\leq(2/|y^*x|)\epsilon\) for \(\epsilon\) sufficiently small.) **b)** Consider a matrix \(A\) with double eigenvalue \(\lambda\) and a perturbation \(\delta A\) of \(A\) as above. Must there be a constant \(C\) such that \(A+\delta A\) has an eigenvalue \(\mu\) for which \(|\lambda-\mu|\leq C\epsilon\) for all sufficiently small \(\epsilon\)? If yes, prove it. If no, construct and analyze an example that shows such a \(C\) need not exist. **c)** Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues that satisfy \(|\lambda_1|>|\lambda_2|>\cdots>|\lambda_m|\). One algorithm for trying to do this is Power Iteration with Rayleigh Quotient approximation of the eigenvalue: Guess \(v^{(0)}\) with \(\|v^{(0)}\|_2=1\). Set \(\lambda^{(0)}=(v^{(0)})^TAv^{(0)}\). For \(k=1,2,3,\ldots\) \[ \begin{aligned} &\text{Let } w=Av^{(k-1)},\\ &v^{(k)}=w/\|w\|_2,\\ &\lambda^{(k)}=(v^{(k)})^TAv^{(k)}. \end{aligned} \] End Analyze the convergence of the approximate eigenvector \(v^{(k)}\) and the approximate eigenvalue \(\lambda^{(k)}\). To what do these quantities converge, and how fast do they converge?

**Iterative Methods** Consider the fixed-point iteration \[ u^{(k+1)}=Tu^{(k)}+c \] for finding a solution of the problem \[ u=Tu+c \] where \(T\) is an \(m\times m\) real matrix and \(c\) is a real \(m\)-vector. It is a fact that the fixed-point iteration converges for any choice of starting point \(u^{(0)}\) if and only if the spectral radius of \(T\), \(\rho(T)\), satisfies \(\rho(T)<1\). **a)** Consider a linear system \(Ax=b\) where \(A\) is strictly diagonally dominant (that is, \(a_{ii}>\sum_{j\ne i}|a_{ij}|\) for all \(i\)). Show that the Gauss-Seidel iterative scheme converges for this problem. Describe a situation in numerical analysis where a strictly diagonally-dominant matrix system arises. **b)** For convergence of fixed-point iteration from any starting point, is it sufficient that \(\|T\|<1\) for some operator norm \(\|\cdot\|\)? Justify your answer.

**Heat Equation Stability** **a)** Consider the initial value problem for the constant-coefficient diffusion equation \[ v_t=\beta v_{xx},\qquad t>0 \] with initial data \(v(x,0)=f(x)\). The Crank-Nicolson scheme for this problem is: \[ \frac{u_j^{n+1}-u_j^n}{k} = \frac{\beta}{2h^2} \left\{ u_{j-1}^{n+1}-2u_j^{n+1}+u_{j+1}^{n+1} +u_{j-1}^n-2u_j^n+u_{j+1}^n \right\}. \] Analyze the 2-norm stability of this scheme and show that the scheme is stable for any choice of \(k>0\) and \(h>0\). **b)** Consider the variable coefficient diffusion equation \[ v_t=(\beta v_x)_x,\qquad 00 \] with Dirichlet boundary conditions \[ v(0,t)=0,\qquad v(1,t)=0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta(x)\geq\beta_0>0\), and that \(\beta(x)\) is smooth. Let \(\beta_{j+1/2}=\beta(x_{j+1/2})\). The Crank-Nicolson scheme for this problem is: \[ \begin{aligned} \frac{u_j^{n+1}-u_j^n}{k} = \frac{1}{2h^2}\{& \beta_{j-1/2}u_{j-1}^{n+1} -(\beta_{j-1/2}+\beta_{j+1/2})u_j^{n+1} +\beta_{j+1/2}u_{j+1}^{n+1}\\ &+\beta_{j-1/2}u_{j-1}^n -(\beta_{j-1/2}+\beta_{j+1/2})u_j^n +\beta_{j+1/2}u_{j+1}^n \}. \end{aligned} \] Analyze the 2-norm stability of this scheme for solving this initial boundary value problem. Do not neglect the fact that there are boundary conditions!

**Elliptic Problems** Consider the standard five-point difference approximation (centered difference for both the gradient and divergence operators) for the variable coefficient Poisson equation \[ -\nabla\cdot(a\nabla v)=f \] with Dirichlet boundary conditions, in a two-dimensional rectangular region. We assume that \(a(x,y)\geq a_0>0\). The approximate solution \(\{u_{i,j}\}\) satisfies a linear system \(Au=b\). **1.** State and prove the maximum principle for the numerical solution \(u_{i,j}\). **2.** Derive the matrix \(A\) in the one-dimensional case and show that it is symmetric and positive definite. **3.** For the one-dimensional and constant-coefficient case, show that the global error \(e_j=v(x_j)-u_j\) satisfies \(\|e\|_2=O(h^2)\) as the space step \(h\to0\). **4.** Discuss the advantages and disadvantages of trying to solve the system for the two-dimensional problem using **(i)** the SOR (Successive Over Relaxation) method and **(ii)** the (preconditioned) Conjugate Gradient method. **Fact 1:** A real symmetric \(n\times n\) matrix \(A\) can be diagonalized by an orthogonal similarity transformation, and \(A\)'s eigenvalues are real. **Fact 2:** The \((N-1)\times(N-1)\) matrix \(M\) defined by \[ M= \begin{bmatrix} -2 & 1 & 0 & \cdots & 0 & 0\\ 1 & -2 & 1 & \ddots & \vdots & \vdots\\ 0 & 1 & -2 & \ddots & 0 & 0\\ \vdots & \ddots & \ddots & \ddots & 1 & 0\\ 0 & \cdots & 0 & 1 & -2 & 1\\ 0 & \cdots & 0 & 0 & 1 & -2 \end{bmatrix} \] has eigenvalues \[ \mu_l=-4\sin^2\left(\frac{\pi l}{2N}\right), \qquad l=1,2,\ldots,N-1. \] **Fact 3:** The \((N+1)\times(N+1)\) matrix: \[ \begin{bmatrix} -1 & 1 & 0 & \cdots & 0 & 0\\ 1 & -2 & 1 & \ddots & \vdots & \vdots\\ 0 & 1 & -2 & \ddots & 0 & 0\\ \vdots & \ddots & \ddots & \ddots & 1 & 0\\ 0 & \cdots & 0 & 1 & -2 & 1\\ 0 & \cdots & 0 & 0 & 1 & -1 \end{bmatrix} \] has eigenvalues \[ \mu_l=-4\sin^2\left(\frac{\pi l}{2(N+1)}\right), \qquad l=0,1,\ldots,N. \] **Fact 4:** For a real \(n\times n\) matrix \(A\), the Rayleigh quotient of a vector \(x\in\mathbb{R}^n\) is the scalar \[ r(x)=\frac{x^TAx}{x^Tx}. \] The gradient of \(r(x)\) is \[ \nabla r(x)=\frac{2}{x^Tx}(Ax-r(x)x). \] If \(x\) is an eigenvector of \(A\) then \(r(x)\) is the corresponding eigenvalue and \(\nabla r(x)=0\).