Preliminary Examination, Numerical Analysis, Fall 2002

Sensitivity Analysis For the linear system \(Ax=b\) where \(A\) is a square, real non-singular matrix, a numerical solution \(\hat {x}\) actually solves a slightly different problem due to roundoff errors: \[ (A+\delta A)\hat {x}=b+\delta b. \]

(a) If we denote \(\delta x=\hat {x}-x\), derive an upper bound for the relative error \[ E_1=\frac {\lVert \delta x\rVert }{\lVert \hat {x}\rVert } \] in terms of the condition number of \(A\), \(\kappa (A)\), and the norms of \(A,\hat {x},\delta A,\delta b\), where \(\lVert \cdot \rVert \) is a suitable vector norm.

(b) We can define another relative error \[ E_2=\frac {\lVert \delta x\rVert }{\lVert x\rVert }, \] based on the exact solution. Derive an upper bound for \(E_2\) in terms of \(\kappa (A)\) and the norms of \(\delta A,A,\delta b\) and \(b\).

(c) Strictly speaking, the relative error is given by \(E_2\), even though in many applications we simply use \(E_1\) since it is more readily available. Point out cases where this approximation may not be appropriate. Explain the significance of the condition number \(\kappa (A)\), and express the relative distance of \(A\) to its nearest singular matrix in terms of \(\kappa (A)\).

Linear Least Squares and Singular Value Decomposition 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.

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. Characterize the difficulties in a rank-deficient least square problem and make a suggestion about which method you should use.

Eigenvalue Problems

(a) For the general eigenvalue problem of a nonsingular matrix \(A\), the simplest method is the power method, which converges to the largest eigenvalue in magnitude, assuming it is simple. The algorithm is the following:

\(i=0\)

repeat \[ \begin {aligned} y_{i+1} &= Ax_i,\\ x_{i+1} &= y_{i+1}/\lVert y_{i+1}\rVert _2,\\ \tilde {\lambda }_{i+1} &= x_{i+1}^T A x_{i+1},\\ i &= i+1 \end {aligned} \] until convergence

The inverse iteration algorithm takes advantage of the fact the if \(\sigma \) (the shift) is closest to one of the simple eigenvalues \(\lambda \), then \(1/(\lambda -\sigma )\) will be the largest eigenvalue in magnitude of the matrix \((A-\sigma I)^{-1}\). Use this fact to construct the inverse iteration algorithm. Note that we will not need to construct the explicit inverse of \(A-\sigma I\). How fast does it converge to the desired eigenvalue? When do you expect the convergence to be slow? How do we find the corresponding eigenvector?

(b) Now consider the symmetric eigenvalue problem where \(A\) is symmetric. What is the Rayleigh quotient? The Rayleigh quotient method is the inverse iteration algorithm described above applied to a symmetric matrix, when the shift is chosen to be the Rayleigh quotient \(\rho _i\) of \(A\) and the updated vector \(x_i\). Explain the rationale of choosing \(\rho _i\) to be the shift and determine the speed of convergence of \(\rho _i\) to \(\lambda \).

Numerical ODEs

(a) Consider the general system of ODEs for \(y(x)\in \mathbb {R}^n\), \(x\in \mathbb {R}\) and \(f=[f_1,\ldots ,f_n]^T\): \[ y'=f(x,y), \] where we assume \(f\) to be Lipschitz continuous. To study the numerical stability, we linearize \(f\) near \(x=x_0,y=y_0\) to obtain: \[ y'\approx f(x_0,y_0)+f_x(x_0,y_0)(x-x_0)+\nabla _y f(x_0,y_0)(y-y_0). \] If we denote \(A=\nabla _y f(x_0,y_0)\), then we claim that the local stability of the above system can be analyzed by studying the following constant coefficient system \[ y'=Ay. \] Justify this claim by analyzing the simple scalar equation \[ y'=a+bx+cy, \] where \(a,b,c\) are constants, in the explicit Euler’s method, to show that the stability is determined by how we solve \(y'=cy\). What is the most important property of \(A\) we need to find that will affect the choice of numerical schemes and step sizes?

(b) 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) 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) 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 at the mesh points \(x_j=j/N\): \[ -u_{j-1}+2u_j-u_{j+1}=f_j=\frac {f(x_j)}{N^2}, \qquad j=1,2,\ldots ,N \] and the central difference for \(u'\) at \(x_N\) to eliminate the ghost point \(x_{N+1}\): \[ u_{N+1}=u_{N-1}. \] Use all the finite difference equations and the boundary conditions to obtain the linear system \(Au=b\), where \[ \setlength {\arraycolsep }{5pt} A= \begin {bmatrix} 2 & -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 & -2 & 2 \end {bmatrix}, \quad u= \begin {pmatrix} u_1\\ u_2\\ \vdots \\ u_{N-1}\\ u_N \end {pmatrix}, \quad b=\frac {1}{N^2} \begin {pmatrix} f_1\\ f_2\\ \vdots \\ f_{N-1}\\ f_N \end {pmatrix}. \]

Verify that \(A\) has eigenvectors \[ v^{(k)}= \begin {bmatrix} \sin \frac {k\pi }{2}x_1\\ \sin \frac {k\pi }{2}x_2\\ \cdot \\ \cdot \\ \sin \frac {k\pi }{2}x_j\\ \cdot \\ \cdot \\ \sin \frac {k\pi }{2}x_{N-1}\\ \sin \frac {k\pi }{2}x_N \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\)?

Wave Equations Consider the first-order wave equation in one space dimension: \[ u_t+au_x=0,\qquad -\infty <x<\infty ,\qquad t>0, \] where \(a\) is a constant, with the initial condition \[ u(x,0)=u_0(x),\qquad -\infty <x<\infty \] for a given function \(u_0\). We use the usual discretization \[ U_j^n\approx u(x_j,t^n), \] where \(x_j=j\Delta x\) and \(t^n=n\Delta t\). For a finite difference scheme to solve this problem, what is the CFL condition? Give an intuitive explanation for the condition. Derive the first-order upwind scheme for this equation and its truncation error. For numerical stability, derive the conditions for stability in the maximum norm (maximum principle) and the \(\ell _2\) norm (Fourier analysis).

Heat Equation Derive the Crank-Nicolson scheme for the following problem \[ u_t=u_{xx},\qquad t>0,\qquad 0<x<1, \] with the boundary conditions \[ u(0,t)=u(1,t)=0,\qquad t\ge 0, \] and the initial data \[ u(x,0)=u_0(x),\qquad 0\le x\le 1, \] for some given function \(u_0\). Use the usual discretization with uniform mesh points \(x_j=jh\) where \(h=1/N\) and \(j=0,1,2,\ldots ,N\). Show that the method is second order by deriving the truncation error. What other conditions do we need to guarantee the convergence of the numerical solution? Derive the stability condition in \(\ell _2\) norm by using Fourier analysis. For any value of \(\nu =\Delta t/(\Delta x)^2\), do we always have the maximum principle satisfied? If we have some initial data such that \(0\le u_0(x)\le 1\), derive the condition that should be imposed on \(\nu \) so that \(0\le u_j^n\le 1\) also holds for all \(j\) and \(n>0\).