Preliminary Examination, Numerical Analysis, Fall 2004

Sensitivity Analysis In this problem, we consider a real, non-singular \(n\times n\) matrix \(A\) and vectors \(f\), \(b\), \(x\in \mathbb {R}^n\).

a) (6%) First we consider the conditioning of matrix-vector multiplication \(f=Ax\), that is, find the smallest factor \(\kappa >0\) such that, when \(x\) is perturbed by any small \(\delta x\), the perturbation for \(f\) satisfies \[ \frac {\lVert \delta f\rVert }{\lVert f\rVert } \le \kappa \frac {\lVert \delta x\rVert }{\lVert x\rVert }. \] Here we assume no perturbation in matrix \(A\).

b) (8%) For the linear system \(Ax=b\), a numerical solution \(\hat {x}\) actually solves a slightly different problem due to roundoff errors: \[ (A+\delta A)\hat {x}=b+\delta b. \] If we denote \(\delta x=\hat {x}-x\), derive an upper bound for the relative error \[ E=\frac {\lVert \delta x\rVert }{\lVert x\rVert } \] in terms of the condition number of \(A\), \(\kappa (A)\), and the norms of \(A\), \(\delta A\), \(b\) and \(\delta b\), where \(\lVert \cdot \rVert \) is a suitable vector norm.

c) (6%) In practice, \(x\) is not available so \(\hat {x}\) is often used in the above bound to generate an approximate bound for \(\delta x\). Comment on this practice as what the advantages and disadvantages are, point out cases where this approximation may be inappropriate.

LS Problems and SVD Let \(A\) be a real \(m\times n\) matrix \((m>n)\) of full rank and let \(b\in \mathbb {R}^m\), the least square problem is to find \(x\in \mathbb {R}^n\) such that \[ \lVert Ax-b\rVert _2 \] is minimized.

a) (13%) Describe the following three approaches to solve the problem: 1) normal equations, 2) QR decomposition, and 3) singular value decomposition (SVD). Discuss briefly the advantages and disadvantages of each approach. In particular, give a geometric interpretation of the normal equation approach, and explain why it is often to be avoided.

b) (7%) In addition to minimizing the above norm, we impose the constraint: \(\lVert x\rVert _2\le \beta \), for some constant \(\beta >0\). This is a constrained LS problem. Use SVD to formulate and solve this problem, notice that there are different cases where a Lagrange multiplier may be needed.

Eigenvalue Problems

a) (8%) QR iteration can be viewed as an extension of the power method. We assume \(n\times n\) real symmetric, nonsingular matrix \(A\), with distinct eigenvalues. Here is the algorithm of orthogonal iteration for two linearly independent vectors \(z_1\) and \(z_2\), listed in an \(n\times 2\) matrix \(Z_0=[z_1,z_2]\):

\(i=0\)

repeat \[ \begin {aligned} Y_{i+1} &= AZ_i,\\ \text {QR Factor }Y_{i+1} &= Z_{i+1}R_{i+1},\\ i &= i+1 \end {aligned} \] until convergence

Notice that \(Y_i,Z_i\) are \(n\times 2\) matrices, and \(R_i\) is a \(2\times 2\) upper triangular matrix. Our assumptions allow us to write \(A=S\Lambda S^{-1}\), where \(\Lambda \) is a diagonal matrix containing all eigenvalues \(\lvert \lambda _1\rvert >\lvert \lambda _2\rvert >\cdots >\lvert \lambda _n\rvert \), \(S=[s_1,s_2,\ldots ,s_n]\) contains all the eigenvectors.

Show that \[ \operatorname {span}\{Z_{i+1}\} = \operatorname {span}\{Y_{i+1}\} = \operatorname {span}\{AZ_i\}, \] so \[ \operatorname {span}\{Z_i\} = \operatorname {span}\{A^iZ_0\} = \operatorname {span}\{S\Lambda ^iS^{-1}Z_0\}. \] Analyze it by writing \[ S\Lambda ^iS^{-1}Z_0=\lambda _2^i S D_i(S^{-1}Z_0). \] Derive the diagonal \(D_i\) and the \(n\times 2\) matrix \(Z'_0=S^{-1}Z_0\), find the limiting span of \(Z_i\) as \(i\to \infty \). How do we obtain the approximating eigenvalues \(\lambda _1,\lambda _2\) and approximating eigenvectors \(s_1\) and \(s_2\), once iterations are terminated?

b) (6%) Assume that we can continue this extension to \(n\) linearly independent vectors, contained in an \(n\times n\) nonsingular \(Z_0\); in particular when \(Z_0=I\), we will have the method of QR iteration. Describe the QR algorithm (without shift) and argue for the rate of convergence based on the analysis in part a).

c) (6%) To speed up convergence, we should use shifting. Describe the QR iteration with a shift, and make a suggestion as what shift should be used. How fast is the convergence?

Numerical ODEs

a) (8%) Consider the general system of ODEs for \(y(x)\in \mathbb {R}^n\), \(x\in \mathbb {R}\), and \[ f=[f_1(x,y),\ldots ,f_n(x,y)]^T: \qquad y'=f(x,y), \] where we assume \(f\) is continuously differentiable in both \(x\) and \(y\). To study local stability, we compare two solutions \(y_1(x)\) and \(y_2(x)\) where \(y_1(x_0)\approx y_2(x_0)\). Derive the linearized system for \(e(x)=y_1(x)-y_2(x)\) near \(x=x_0\) to show that \[ e'\approx Ae \] and find this square matrix \(A\). What is the most important property of \(A\) we need to find that will affect the choice of numerical schemes and step sizes?

b) (12%) Now consider the initial value problem \[ y''+2\mu y'+(1+\mu ^2)y=0, \] \[ y(0)=0,\qquad y'(0)=1, \] where \(\mu >0\) is a given constant. If we use a simple explicit Euler’s method, determine the maximum step size \(h\) that is allowed for the solution to be stable.

Iterative Methods for Linear Systems

a) (6%) Consider the linear system \(Ax=b\) and a splitting of \(A:A=M-K\), with \(M\) being nonsingular, an iterative method is to iterate \[ x_{m+1}=Rx_m+c, \] where \(R=M^{-1}K\) and \(c=M^{-1}b\). Show that the method converges for any \(x_0\) if \(\lVert R\rVert <1\), where \(\lVert \cdot \rVert \) is any operator norm.

b) (14%) To solve the boundary value problem \[ -u''(x)=f(x),\qquad \text {for }0<x<1, \] with \(u'(0)=0\) and \(u(1)=0\), we use second-order central differences for the discretization of the equation at mesh points \(x_j=j/N\): \[ -u_{j-1}+2u_j-u_{j+1}=f_j=\frac {f(x_j)}{N^2}, \qquad j=0,1,\ldots ,N-1 \] and the central difference approximation for \(u'\) at \(x_0=0\): \(u_1-u_{-1}=0\), to eliminate the ghost point.

Show that the finite difference equations and the boundary conditions lead to the following linear system \(Au=b\), where \[ \setlength {\arraycolsep }{5pt} A= \begin {bmatrix} 1 & -1 & 0 & \cdots & 0 & 0\\ -1 & 2 & -1 & \cdots & 0 & 0\\ 0 & -1 & 2 & \ddots & 0 & 0\\ \vdots & \vdots & \ddots & \ddots & -1 & 0\\ 0 & 0 & 0 & -1 & 2 & -1\\ 0 & 0 & 0 & 0 & -1 & 2 \end {bmatrix}, \quad u= \begin {pmatrix} u_0\\ u_1\\ \vdots \\ u_{N-2}\\ u_{N-1} \end {pmatrix}, \quad b=\frac {1}{N^2} \begin {pmatrix} \frac {1}{2}f_0\\ f_1\\ \vdots \\ f_{N-2}\\ f_{N-1} \end {pmatrix}. \] Notice that equations are scaled so that \(A\) is symmetric. Determine if \(A\) has eigenvectors \[ v^{(k)}= \begin {bmatrix} 1\\ \cos \frac {1}{2}k\pi x_1\\ \vdots \\ \cos \frac {1}{2}k\pi x_j\\ \vdots \\ \cos \frac {1}{2}k\pi x_{N-2}\\ \cos \frac {1}{2}k\pi x_{N-1} \end {bmatrix}, \qquad k=1,3,5,\ldots ,2N-3,2N-1, \] and use them to find the eigenvalues of the iteration matrix \(R\) for the Jacobi method. Analyze the convergence properties of the Jacobi method applied to this problem, and express the convergence speed as a function of \(N\). How does the number of iterations required to reduce the initial error by a factor of \(\delta \) depend on \(N\)? In reality, should you use an iterative method to solve this matrix problem? If not, what other method is preferred?

Wave Equations Consider the standard wave equation in one space dimension: \[ u_{tt}=c^2u_{xx},\qquad 0<x<1,\quad t>0, \] with initial conditions \[ u(x,0)=f(x),\qquad u_t(x,0)=g(x),\qquad 0<x<1 \] for given smooth functions \(f\) and \(g\), and boundary condition \[ u(0,t)=u(1,t)=0,\qquad t>0. \] The leap-frog scheme takes advantage of the fact that the equation can be reformulated as a system of first order equations \[ u_t+cv_x=0, \qquad v_t+cu_x=0. \] Denote \(x_j=j\Delta x\), \(U_j^n\approx u(x_j,n\Delta t)\), and \(V_{j+1/2}^{n+1/2}\approx v(x_{j+1/2},(n+\frac {1}{2})\Delta t)\); use central difference approximations for both derivatives to derive the leap-frog scheme for the system of equations. How do you implement the initial and boundary conditions for both \(u\) and \(v\)? Derive the truncation error for this approximation and determine the order of truncation error. For stability, consider solutions in the form \[ \begin {bmatrix} U^n\\ V^{n+1/2} \end {bmatrix} = \lambda ^n e^{ikx} \begin {bmatrix} \widetilde {U}\\ \widetilde {V} \end {bmatrix}, \] for some constants \(\widetilde {U}\) and \(\widetilde {V}\). Substitute this solution into the finite difference equations, find out what condition \(\lambda \) has to satisfy in order for \(\widetilde {U}\) and \(\widetilde {V}\) to have nontrivial solutions. For those \(\lambda \) values, compute \(\lvert \lambda \rvert \) to see if the stability condition \(\lvert \lambda \rvert \le 1\) is satisfied.

Convection-Diffusion Equation For the convection-diffusion equation \[ u_t+au_x=bu_{xx}, \] where \(a>0\) and \(b>0\), a natural Crank-Nicolson scheme is \[ \begin {aligned} \frac {U_j^{n+1}-U_j^n}{\Delta t} &+\frac {a}{2}\left ( \frac {U_{j+1}^n-U_{j-1}^n}{2\Delta x} + \frac {U_{j+1}^{n+1}-U_{j-1}^{n+1}}{2\Delta x} \right )\\ &= \frac {b}{2}\left ( \frac {U_{j+1}^n-2U_j^n+U_{j-1}^n}{(\Delta x)^2} + \frac {U_{j+1}^{n+1}-2U_j^{n+1}+U_{j-1}^{n+1}}{(\Delta x)^2} \right ). \end {aligned} \] Denote \(\nu =b\Delta t/(\Delta x)^2\) and \(\mu =a\Delta t/(2\Delta x)\); show that the scheme is stable in the max norm if \(\mu \le \nu \le 1\). Also use Fourier analysis to derive the stability condition in the \(\ell _2\) norm and show that the scheme is always stable in the \(\ell _2\) norm. Suppose \(u\) stands for some chemical concentration, which norm is more appropriate? Explain why this scheme is so restrictive that it is not often used.