Preliminary Examination, Numerical Analysis, Fall 1997

Sensitivity of Solution of Linear System Suppose that \(A\) is a real non-singular \(n\times n\) matrix. Consider the family of linear systems \(Ax(\epsilon )=b+\epsilon c\), where \(b\) and \(c\) are given nonzero vectors in \(\mathbb {R}^n\), and \(\epsilon \) is a real parameter. Let \(x=x(0)\) be the solution for \(\epsilon =0\). Let \(\lVert \cdot \rVert _2\) denote the Euclidean vector norm. Show that \[ \frac {\lVert x(\epsilon )-x\rVert _2}{\lVert x\rVert _2} \le \lVert A\rVert _2\lVert A^{-1}\rVert _2 \frac {\lVert \epsilon c\rVert _2}{\lVert b\rVert _2}. \] The expression \(\lVert A\rVert _2\lVert A^{-1}\rVert _2\) is the condition number of \(A\). Interpret the above inequality and explain its significance?

Singular Value Decomposition Suppose that \(A\) is a real \(m\times n\) matrix with \(m\ge n\). What is the singular value decomposition (SVD) of \(A\)? Does every matrix \(A\) of this type have an SVD? If your answer is yes, prove the assertion. If your answer is no, give an example of a matrix \(A\) which has no SVD. For a matrix \(A\) which does have an SVD, how would having the SVD give easy access to information about the rank of \(A\)? If \(A\) were \(n\times n\) and nonsingular, how would knowing \(A\)’s SVD help you determine the sensitivity of the solution of a linear system of the following form \[ Ax=b? \]

Eigenvalue Problems 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 \(\lVert x\rVert _2=\lVert y\rVert _2=1\). a) Show that the condition number for \(\lambda \) to changes in \(A\) is given by \(s(\lambda )=1/\lvert y^*x\rvert \). b) Give an example of a matrix \(A\) which has a badly conditioned eigenvalue, and another matrix \(A\) which has only well-conditioned eigenvalues. c) Sketch a reasonably efficient algorithm for finding the largest eigenvalue of a symmetric matrix \(A\) with distinct eigenvalues. How fast do you expect this algorithm to converge?

Linear Least Squares The Linear Least Squares problem for an \(m\times n\) real matrix \(A\) and \(b\in \mathbb {R}^m\) is the problem: \[ \text {Find }x\in \mathbb {R}^n\text { such that }\lVert Ax-b\rVert _2\text { is minimized.} \]

a) Suppose that you have data \(\{(t_j,y_j)\}\), \(j=1,2,\ldots ,m\) that you wish to approximate by an expansion \[ p(t)=\sum _{k=1}^n x_k\phi _k(t). \] Here, the functions \(\phi _k(t)\) are given functions. Which norm on the difference between the approximation function \(p\) and the data gives rise to a linear least squares problem for the unknown expansion coefficients \(x_k\)? What is the matrix \(A\) in this case, and what is the vector \(b\)?

b) Suppose that \(A\) is a real \(m\times n\) matrix of full rank and let \(b\in \mathbb {R}^m\). What are the ‘normal equations’ for the Least Squares problem? How can they be used to solve the Least Squares problem? What is the QR factorization of \(A\) and how can it be used to solve the Least Squares problem? Compare and contrast these approaches for numerically solving the Least Square problem.

Iterative Methods for Linear Systems Consider the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0\le x\le 1 \] with \(u(0)=u(1)=0\), and the following discretization of it: \[ -U_{j-1}+2U_j-U_{j+1}=F_j \] for \(j=1,2,\ldots ,N-1\) where \(Nh=1\) and \(F_j=h^2f(jh)\).

Analyze the convergence properties of both the Jacobi and Gauss Seidel iterative methods for this problem. In particular, for each, express the speed of convergence as a function of the discretization stepsize \(h\). How does the number of iterations required to reduce the initial error by a factor \(\delta \) depend on \(h\) for each method? In practice, would you use one of these methods to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Newton’s Method Suppose \(f(x)\) is a real-valued function defined for all values of the real variable \(x\), and that \(f'\) is Lipschitz continuous. State and prove a local convergence theorem for Newton’s method for finding a simple zero of \(f\). How does this result change if the root is a double root? Is Newton’s method globally convergent? If you answered yes, prove your claim. If you answered no, sketch an example which illustrates your claim, and outline a method that has the Newton method’s local convergence properties but which is globally convergent.

Broyden’s Method Consider the nonlinear system: \[ F(x)=0\qquad \text {where }F:\mathbb {R}^n\to \mathbb {R}^n \] and the iteration \[ x^{(k+1)}=x^{(k)}-(A^{(k)})^{-1}F(x^{(k)}) \] where \(x^{(0)}\in \mathbb {R}^n\) is given.

Newton’s Method is defined by the choice \(A^{(k)}=F'(x)\) where \(F'(x)\) is the Jacobean matrix of \(F\). An alternative choice for \(A^{(k+1)}\) is made in Broyden’s Method and is motivated as follows:

Let \(A^{(0)}\) be a given and suitably chosen matrix. The scalar Secant Method suggests imposing on \(A^{(k+1)}\) the condition \[ A^{(k+1)}(x^{(k+1)}-x^{(k)})=F(x^{(k+1)})-F(x^{(k)}). \] This however does not uniquely determine \(A^{k+1}\) as it does not give any information on how \(A^{k+1}\) acts on vectors perpendicular to \((x^{(k+1)}-x^{(k)})\). Thus require also that: \[ A^{(k+1)}t=A^{(k)}t \qquad \text {for all }t\in \mathbb {R}^n\text { such that } t^T(x^{(k+1)}-x^{(k)})=0 \] Use these conditions to derive an explicit expression for \(A^{(k+1)}\) in terms of \(A^{(k)}\), \(x^{(k)}\), and \(x^{(k+1)}\).

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}hf_{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\) for all \(x\).

a) Analyze the consistency, stability, accuracy, and convergence properties of this method.

b) Would it be more reasonable to use this method or Euler’s method for the initial value problem: \[ y'=-10^8[y-\cos (x)]-\sin (x),\qquad y(0)=1? \] Justify your answer.