Preliminary Examination, Numerical Analysis, January 2010

**Least Squares Problems** **a)** For a full rank real \(m\times n\) matrix \(A\), show that \(X=A^\dagger\), the pseudoinverse of \(A\), minimizes \(\|AX-I\|_F\) over all \(n\times m\) matrices \(X\). What is the value of the minimum? (Hint: Relate the problem to a set of least-squares problems). **b)** For a real full rank \(m\times n\) matrix \(A\) and vector \(b\in\mathbb{R}^m\), explain how to solve the least-squares problem of finding \(x\in\mathbb{R}^n\) that minimizes \(\|Ax-b\|_2\) using **i)** the normal equations, and **b)** a QR factorization of the matrix \(A\). What are the advantages and disadvantages of each of these methods?

**Interpolation and Integration:** **a)** Consider equally spaced points \(x_j=a+jh\), \(j=0,\ldots,n\) on the interval \([a,b]\), where \(nh=b-a\). Let \(f(x)\) be a smooth function defined on \([a,b]\). Show that there is a unique polynomial \(p(x)\) of degree \(n+1\) which interpolates \(f\) at all of the points \(x_j\). Derive the formula for the interpolation error at an arbitrary point \(x\) in the interval \([a,b]\): \[ f(x)-p(x)\equiv E(x) = \frac{1}{(n+1)!}(x-x_0)(x-x_1)\cdots(x-x_n)f^{n+1}(\eta). \] for some \(\eta\in[a,b]\). **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)\).

**Sensitivity:** Consider a \(6\times6\) symmetric positive definite matrix \(A\) with singular values \(\sigma_1=1000\), \(\sigma_2=500\), \(\sigma_3=300\), \(\sigma_4=20\), \(\sigma_5=1\), \(\sigma_6=0.01\). **a)** Suppose you use a Cholesky factorization package on a computer with a machine epsilon \(10^{-14}\) to solve the system \(Ax=b\) for some nonzero vector \(b\). How many digits of accuracy do you expect in the computed solution? Justify your answer in terms of condition and stability. You may assume that the entries of \(A\) and \(b\) are exactly represented in the computer's floating-point system. **b)** Suppose that instead you use an iterative method to find an approximate solution to \(Ax=b\) and you stop iterating and accept iterate \(x^{(k)}\) when the residual \(r^{(k)}=Ax^{(k)}-b\) has 2-norm less than \(10^{-9}\). Give an estimate of the maximum size of the relative error in the final iterate? 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)\). A scheme for this problem is: \[ \frac{u_j^{n+1}-u_j^n}{k} = \frac{\beta}{h^2} \left\{u_{j-1}^{n+1}-2u_j^{n+1}+u_{j+1}^{n+1}\right\}. \] Analyze the 2-norm stability of this scheme. For which values of \(k>0\) and \(h>0\) is the scheme stable? (Note that this there are no boundary conditions here.) **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})\). A scheme for this problem is: \[ \frac{u_j^{n+1}-u_j^n}{k} = \frac{1}{h^2} \left\{ \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} \right\}. \] 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!

{Numerical Methods for ODEs:} Consider the Linear Multistep Method \[ y_{n+2}-\frac{4}{3}y_{n+1}+\frac{1}{3}y_n=\frac{2}{3}k f_{n+2} \] for solving an initial value problem \(y'=f(y,x)\), \(y(0)=\eta\). You may assume that \(f\) is Lipschitz continuous with respect to \(y\) uniformly for all \(x\). **a)** Analyze the consistency, stability, accuracy, and convergence properties of this method. **b)** Sketch a graph of the solution to the following initial value problem. \[ y'=-10^8[y-\cos(x)]-\sin(x),\qquad y(0)=2. \] Would it be more reasonable to use this method or Euler's method for this problem? What would you consider in choosing a timestep \(k\) for each of the methods? Justify your answer. **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^T Ax}{x^T x}. \] The gradient of \(r(x)\) is \[ \nabla r(x)=\frac{2}{x^T x}(Ax-r(x)x). \] If \(x\) is an eigenvector of \(A\) then \(r(x)\) is the corresponding eigenvalue and \(\nabla r(x)=0\).