Preliminary Examination, Numerical Analysis, August 2005

Ordinary Differential Equations Consider the initial value problem for a system of two ODEs \[ y'=Ay \] where \(A\) has eigenvalues \(-10^4\) and \(1\). Which of the following schemes would be the best choice for solving this problem? Justify your answer in terms of stability, accuracy, and efficiency. \[ y^{n+1}=y^n+hAy^n, \tag {1} \] \[ y^{n+1}=y^n+hAy^{n+1}, \tag {2} \] \[ y^{n+1}=y^n+\frac {h}{2}(Ay^n+Ay^{n+1}). \tag {3} \]

Wave Equation Consider the Lax-Friedrichs scheme below for the PDE \(v_t+cv_x=0\), where you may assume that \(c>0\): \[ u_j^{n+1}=\frac {1}{2}(u_{j-1}^n+u_{j+1}^n) -\frac {ck}{2h}(u_{j+1}^n-u_{j-1}^n), \] where \(h\) is the spacestep and \(k\) is the timestep. Note that the scheme can also be written \[ u_j^{n+1}=u_j^n-\frac {ck}{2h}(u_{j+1}^n-u_{j-1}^n) +\frac {1}{2}(u_{j-1}^n-2u_j^n+u_{j+1}^n). \]

a) Analyze the stability of this scheme in terms of \(\alpha =\frac {ck}{h}\).

b) Recall that a modified equation for a scheme is the PDE satisfied, to a certain order, by a smooth solution \(u(x,t)\) of the difference scheme. To zeroth order in \(k\) and \(h\), the modified equation for the Lax-Friedrichs scheme is \(u_t+cu_x=0\). Derive the modified equation through terms of first order in \(k\) and \(h\). You may assume that \(\alpha \) is constant. Does this modified equation predict the same stability condition as your analysis in part a)? Explain your answer. What does the modified equation suggest about the numerical solution that will be obtained for the step-function initial data \(u(x,0)=1\) for \(x\le 0\) and \(u(x,0)=0\) for \(x>0\)?

Heat Equation Consider the variable coefficient diffusion equation \[ v_t=(\beta v_x)_x,\qquad 0<x<1,\quad t>0 \] with Dirichlet boundary conditions \[ v(0,t)=0.0,\qquad v(1,t)=0.0 \] and initial data \(v(x,0)=f(x)\). Assume that \(\beta (x)\ge \beta _0>0\), and that \(\beta (x)\) is smooth. For any fixed value of \(\theta \in [0,1]\), consider the \(\theta \)-scheme for this problem: \[ \begin {aligned} \frac {u_j^{n+1}-u_j^n}{k} =\frac {1}{h^2}\{& (1-\theta )(\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})\\ &+\theta (\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} \] with appropriate treatment of the boundary conditions in the relevant equations. Here, \(h\) is the spacestep, \((J+1)h=1\), \(k\) is the timestep, and \(\beta _{j+1/2}=\beta (x_{j+1/2})\), where \(x_{j+1/2}=(j+1/2)h\). Let \(M\) be the symmetric \(J\times J\) matrix defined by \[ \setlength {\arraycolsep }{3pt} M= \begin {pmatrix} -(\beta _{1/2}+\beta _{3/2}) & \beta _{3/2} & \cdots & \cdots & \cdots & \cdots \\ \beta _{3/2} & -(\beta _{3/2}+\beta _{5/2}) & \beta _{5/2} & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \beta _{J-3/2} & -(\beta _{J-3/2}+\beta _{J-1/2}) & \beta _{J-1/2}\\ \cdots & \cdots & \cdots & \cdots & \beta _{J-1/2} & -(\beta _{J-1/2}+\beta _{J+1/2}) \end {pmatrix}. \] Analyze the stability of this scheme in terms of \(\theta \), \(\frac {k}{h^2}\), and the eigenvalues of \(M\). Are there values of \(\theta \) for which the scheme is stable for all choices of \(\frac {k}{h^2}>0\)? (Note that answering this question should not require lots of writing.)

Sensitivity Consider a \(6\times 6\) 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.

Iterative Methods Consider the fixed-point problem \(x=Tx+c\), where \(T\) is a real \(N\times N\) matrix and \(x\) and \(c\) are real \(N\)-vectors.

a) Show that the fixed point iteration \[ x^{(k+1)}=Tx^{(k)}+c \] converges for any choice of initial vector \(x^{(0)}\) if and only if the spectral radius of \(T\), \(\rho (T)\), satisfies \(\rho (T)<1\).

b) Consider the linear system \(Ax=b\) where \(A\) is a strictly diagonally dominant matrix, that is, \(a_{ii}>\sum _{j\ne i}\lvert a_{ij}\rvert \) for all \(i\), and strict inequality holds for at least one value of \(i\). Show that the Jacobi iterative scheme converges for this problem.

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)=\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)=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)\).

Projectors Recall that a square matrix \(P\) that satisfies \(P^2=P\) is called a projector and that the subspaces \(\operatorname {range}(P)\) and \(\operatorname {null}(P)\) are complementary subspaces. Recall also that a projector \(P\) is an orthogonal projector if these two subspaces are perpendicular to one another, or equivalently, if \(P\) is Hermitian.

a) Let \(P\) be a nonzero projector. Show that \(\lVert P\rVert _2\ge 1\), with equality if and only if \(P\) is an orthogonal projector.

b) Given an \(m\times n\) matrix \(A\) \((m\ge n)\) with rank \(r\). Describe a stable algorithm for constructing an orthogonal projector onto \(\operatorname {range}(A)\). Why do you think your algorithm is stable?

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 \(\lVert x\rVert _2=\lVert y\rVert _2=1\). Show that the condition number for \(\lambda \) to changes in \(A\) is given by \(s(\lambda )=1/\lvert y^*x\rvert \). Give an example of a matrix \(A\) which has a badly conditioned eigenvalue, and another matrix \(A\) which has only well-conditioned eigenvalues.

b) Consider the problem of computing an eigenvalue/eigenvector pair of a real symmetric \(m\times m\) matrix \(A\) with eigenvalues of distinct magnitudes. One algorithm for trying to do this is Rayleigh Quotient Iteration:

Guess \(v^{(0)}\) with \(\lVert v^{(0)}\rVert _2=1\).

Set \(\lambda ^{(0)}=(v^{(0)})^TAv^{(0)}\).

For \(k=1,2,3,\ldots \): \[ \text {Solve }(A-\lambda ^{(k-1)}I)w=v^{(k-1)}\text { for }w, \] \[ v^{(k)}=\frac {w}{\lVert w\rVert _2}, \] \[ \lambda ^{(k)}=(v^{(k)})^TAv^{(k)}. \] Until convergence.

Suppose that the Rayleigh Quotient Iteration iterates \(\lambda ^{(k)}\) converge to a simple eigenvalue \(\lambda \) of \(A\). Analyze the speed of convergence of \(\lambda ^{(k)}\) to \(\lambda \).