Iterative Methods for Linear Systems

Fixed-point iteration for linear systems.

1. For \(T \in \mathbb {R}^{m \times m}\) and any matrix norm, show that \(\|T^k\| \to 0\) as \(k \to \infty \) if and only if \(\rho (T) < 1\) where \(\rho (T)\) is the spectral radius of \(T\). Do not assume that the matrix \(T\) is diagonalizable.
2. Suppose \(A \in \mathbb {R}^{m \times m}\) is nonsingular, that \(b \in \mathbb {R}^m\), and that \(Ax = b\). Suppose that \(x = T x + c\) is a fixed point problem equivalent to this linear system. What are the implications of part (a) for the convergence of the fixed-point iteration \(x^{(k+1)} = T x^{(k)} + c\) to the solution \(x\) of the linear system?

3. Suppose \(A\) is a symmetric, weakly diagonally dominant, and irreducible matrix. Define the Jacobi iterative method for solving \(Ax = b\), and show that this method converges to the solution of the linear system starting with any initial vector \(x^{(0)}\). (Hint: You might find it helpful to use the fact that a matrix which is both weakly diagonally dominant and irreducible is nonsingular.)

   Recall that: An \(m \times m\) matrix \(A\) is irreducible iff for all \(i,j \in \{1,2,\ldots ,m\}\), \(i \ne j\), either \(a_{ij} \ne 0\) or there are \(i_1,i_2,\ldots ,i_s\), each between 1 and \(m\), such that \(a_{i,i_1}a_{i_1,i_2}a_{i_2,i_3}\cdots a_{i_s,j} \ne 0\). A symmetric \(m \times m\) matrix \(A\) is weakly diagonally dominant if \(|a_{ii}| \ge \sum _{j\ne i}|a_{ij}|\) for all \(i\), and \(|a_{ii}| > \sum _{j\ne i}|a_{ij}|\) for at least one \(i\).

Fixed-point iteration for linear systems.

1. For \(T \in \mathbb {R}^{m \times m}\) and any matrix norm, show that \(\|T^k\| \to 0\) as \(k \to \infty \) if and only if \(\rho (T) < 1\) where \(\rho (T)\) is the spectral radius of \(T\). Do not assume that the matrix \(T\) is diagonalizable.
2. Suppose \(A \in \mathbb {R}^{m \times m}\) is nonsingular, that \(b \in \mathbb {R}^m\), and that \(Ax = b\). Suppose that \(x = T x + c\) is a fixed point problem equivalent to this linear system. What are the implications of part (a) for the convergence of the fixed-point iteration \(x^{(k+1)} = T x^{(k)} + c\) to the solution \(x\) of the linear system?

3. Suppose \(A\) is a symmetric, weakly diagonally dominant, and irreducible matrix, define the Jacobi iterative method for solving \(Ax = b\), and show that this method converges to the solution of the linear system starting with any initial vector \(x^{(0)}\). (Hint: You might find it helpful to use the fact that a matrix which is both weakly diagonally dominant and irreducible is nonsingular.)

   Recall that: An \(m \times m\) matrix \(A\) is irreducible iff for all \(i,j \in \{1,2,\ldots ,m\}\), \(i \ne j\), either \(a_{ij} \ne 0\) or there are \(i_1,i_2,\ldots ,i_s\), each between 1 and \(m\), such that \(a_{i,i_1}a_{i_1,i_2}a_{i_2,i_3}\cdots a_{i_s,j} \ne 0\). A symmetric \(m \times m\) matrix \(A\) is weakly diagonally dominant if \(|a_{ii}| \ge \sum _{j\ne i}|a_{ij}|\) for all \(i\), and \(|a_{ii}| > \sum _{j\ne i}|a_{ij}|\) for at least one \(i\).

(Iterative Solvers).

1. State Gauss-Seidel method for the solution of the linear system \(Ax = b\), where \(A \in \mathbb {R}^{m \times m}\).
2. Show that if \(A\) is a strictly diagonally dominant matrix by rows, the Gauss-Seidel method is convergent for any \(x^{(0)}\).

(Moore-Penrose pseudo-inverse)

Let \(A \in \mathbb {C}^{m \times n}\) have rank \(r\) and reduced SVD \[ A = \tilde {U}\tilde {\Sigma }\tilde {V}^* . \]

The Moore-Penrose pseudoinverse of \(A\) is defined as

a. In the induced matrix \(\ell ^2\) norm \(\|\cdot \|\), is the operation \(A \mapsto A^+\) well-conditioned? That is, given an arbitrary but fixed \(A\) and a perturbation matrix \(B\), is \(\|(A + B)^+ - A^+\|/\|A^+\|\) controllable by \(\|B\|/\|A\|\)? Prove it, or give a counterexample.

b. If \(m \ge n\) and \(A\) has full rank \(n\) and \(b \in \mathbb {C}^m\) is a given vector, prove that \(x = A^+b\) is identical to the least-squares solution \(z\) of the over-determined linear system \(Az = b\).

Let \(X\) be a real Banach space, \(f:X\to \mathbb {R}\) be a bounded linear functional and \(T:X\to X\) a bounded linear operator. In the following \(y\in X\) is fixed and \(\alpha \in \mathbb {R}\).

1. Prove that there exists a constant \(C>0\) such that if \(\lvert \alpha \rvert <C\), then the non-linear equation \[ x+\alpha f(x)Tx=y, \tag {1} \] has a unique solution \(x\) in the ball \(B=\{x\in X\mid \lVert x-y\rVert \le 1\}\).
2. Suggest an iterative procedure for approximating the unique solution \(x\) to equation (1).

CG Method

a) Formulate the Conjugate Gradient (CG) algorithm for the solution of the linear system \(Ax=b\), where \(A\in \mathbb {R}^{n\times n}\) is symmetric positive definite matrix.

b) Show that for any initial guess \(x_0\in \mathbb {R}^n\), Conjugate Gradient algorithm converges to the solution \(x\) of the linear system \(Ax=b\) in at most \(n\) steps.

Overdetermined Linear System

Solve the following overdetermined linear system in the sense of the least squares: \[ z+y=1,\qquad z-y=2,\qquad 4z+2y=4.8. \]

Let \(X\) be a Banach space, \(f:X\to \mathbb {R}\) be a bounded linear functional, \(y\in X\) fixed and \(\alpha \in \mathbb {R}\).

1. Prove that there exists a constant \(C>0\) such that if \(\lvert \alpha \rvert <C\), then the non-linear equation \[ x+\alpha f(x)x=y, \tag {2} \] has a unique solution \(x\) in the ball \(B=\{x\in X\mid \lVert x-y\rVert \le 1\}\).
2. Suggest an iterative procedure for approximating the unique solution \(x\) to equation (2).

Iterative Methods for Linear Systems Consider the square linear system \(Ax=b\) and the Gauss-Jacobi Method \[ x_i^{[k+1]}= \frac {1}{a_{ii}}\left ( b_i-\sum _{j=1}^{i-1}a_{ij}x_j^{[k]} -\sum _{j=i+1}^{n}a_{ij}x_j^{[k]} \right ) \] where \(i=1,\ldots ,n\) and \(k=0,1,\ldots \). Show that this method converges for all starting vectors \(x^{[0]}\) if \(A\) is diagonally dominant, i.e., \[ \sum _{j=1}^{k}\lvert a_{ij}\rvert + \sum _{j=k+1}^{n}\lvert a_{ij}\rvert < \lvert a_{ii}\rvert \qquad \text {for all } i=1,2,\ldots ,n. \]

Consider the (square) linear system \(Ax=b\). Suppose you have obtained a numerical approximation \(\hat {x}\). Define the error and residual, respectively, by \[ e=x-\hat {x}\qquad \text {and}\qquad r=b-A\hat {x}. \] Let \(\|\cdot \|\) denote a vector norm, and the induced matrix norm. Show that \[ \frac {1}{\|A\|\|A^{-1}\|}\frac {\|r\|}{\|b\|} \leq \frac {\|e\|}{\|x\|} \leq \|A\|\|A^{-1}\|\frac {\|r\|}{\|b\|} \] and discuss the significance of these inequalities.

Consider the variational principle: \[ I(u)=\int _0^1 \left (u'(x)\right )^2+k^2\left (u(x)\right )^2 dx=\min \qquad \text {subject to}\qquad u(0)=u(1)=1. \] Suppose you want to use the Ritz Method to approximate the solution of this problem by a polynomial of degree \(n\). Derive a linear system for the coefficients of that polynomial.

Let \(T:L^1[0,1]\to L^2[0,1]\) be a bounded linear operator. Let \(a\in L^\infty [0,1]\) and \(g\in L^1[0,1]\).

1. Prove that there exists a constant \(C>0\) such that if \(\lVert a\rVert _{L^\infty }<C\) then the equation \[ f+a(Tf)^2=g \tag {1} \] has a unique solution \(f\) in the set \[ B=\{u\in L^1[0,1]:\lVert u-g\rVert _{L^1}\le 1\}. \] Note: the constant \(C\) may depend on \(T\) and/or \(g\).
2. Suggest an iterative procedure for approximating the unique solution \(f\) to equation (1).

Iterative Methods for Linear Systems

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

Analyze the convergence properties of the Jacobi iterative method for this problem. In particular, 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\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Elliptic Problems:

Consider the standard five-point difference approximation (centered difference for both the gradient and divergence operators) for the variable coefficient Poisson equation \[ -\nabla \cdot (a\nabla v)=f \] with Dirichlet boundary conditions, in a two-dimensional rectangular region. We assume that \(a(x,y)\geq a_0>0\). The approximate solution \(\{u_{i,j}\}\) satisfies a linear system \(Au=b\).

1. State and prove the maximum principle for the numerical solution \(u_{i,j}\).

2. Derive the matrix \(A\) in the one-dimensional case and show that it is symmetric and positive definite.

3. For the one-dimensional and constant-coefficient case, show that the global error \(e_j=v(x_j)-u_j\) satisfies \(\|e\|_2=O(h^2)\) as the space step \(h\to 0\).

4. Discuss the advantages and disadvantages of trying to solve the system for the two-dimensional problem using (i) the SOR (Successive Over Relaxation) method and (ii) the (preconditioned) Conjugate Gradient method.

Iterative Methods for Linear Systems

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

Analyze the convergence properties of the Jacobi iterative method for this problem. In particular, 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\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Elliptic Problems:

Consider the standard five-point difference approximation (centered difference for both the gradient and divergence operators) for the variable coefficient Poisson equation \[ -\nabla \cdot (a\nabla v)=f \] with Dirichlet boundary conditions, in a two-dimensional rectangular region. We assume that \(a(x,y)\geq a_0>0\). The approximate solution \(\{u_{i,j}\}\) satisfies a linear system \(Au=b\).

1. State and prove the maximum principle for the numerical solution \(u_{i,j}\).

2. Derive the matrix \(A\) in the one-dimensional case and show that it is symmetric and positive definite.

3. For the one-dimensional and constant-coefficient case, show that the global error \(e_j=v(x_j)-u_j\) satisfies \(\|e\|_2=O(h^2)\) as the space step \(h\to 0\).

4. Discuss the advantages and disadvantages of trying to solve the system for the two-dimensional problem using (i) the SOR (Successive Over Relaxation) method and (ii) the (preconditioned) Conjugate Gradient method.

Error Analysis Consider the linear system \[ Ax=b \tag {1} \] where \(A\) is an invertible matrix. Let \(\hat {x}\) be an approximation of the solution of (1) and let \[ e=x-\hat {x}\qquad \text {and}\qquad r=b-A\hat {x}. \] Let \(\|\cdot \|\) denote any vector norm or the corresponding induced matrix norm. Show that \[ \frac {\|e\|}{\|x\|} \leq \|A\|\|A^{-1}\|\frac {\|r\|}{\|b\|}. \] Comment on the significance of the condition number \(\|A\|\|A^{-1}\|\) and give a lower bound for it in terms of the eigenvalues of \(A\).

Elliptic Problems:

Consider the standard five-point difference approximation (centered difference for both the gradient and divergence operators) for the variable coefficient Poisson equation \[ -\nabla \cdot (a\nabla v)=f \] with Dirichlet boundary conditions, in a two-dimensional rectangular region. We assume that \(a(x,y)\geq a_0>0\). The approximate solution \(\{u_{i,j}\}\) satisfies a linear system \(Au=b\).

1. State and prove the maximum principle for the numerical solution \(u_{i,j}\).

2. Derive the matrix \(A\) in the one-dimensional case and show that it is symmetric and positive definite.

3. For the one-dimensional and constant-coefficient case, show that the global error \(e_j=v(x_j)-u_j\) satisfies \(\|e\|_2=O(h^2)\) as the space step \(h\to 0\).

4. Discuss the advantages and disadvantages of trying to solve the system for the two-dimensional problem using (i) the SOR (Successive Over Relaxation) method and (ii) the (preconditioned) Conjugate Gradient method.

Sensitivity Analysis

(a) Let \(X\in \mathbb {R}^{n\times n}\) such that \(\lVert X\rVert <1\) for some induced norm. Show that \(I-X\) is nonsingular and \[ \lVert (I-X)^{-1}\rVert \le \frac {1}{1-\lVert X\rVert }. \]

(b) Assume that \(A\in \mathbb {R}^{n\times n}\) is nonsingular, and the linear system \[ Ax=b \] is perturbed to \[ (A+\delta A)(x+\delta x)=b+\delta b, \] and the perturbations are small enough so that \[ \lVert \delta A\rVert \le \epsilon \lVert A\rVert ,\qquad \lVert \delta b\rVert \le \epsilon \lVert b\rVert , \] for some small \(\epsilon >0\), and \(\epsilon \kappa (A)=r<1\), where \(\kappa (A)\) is the condition number of \(A\). Show that \(A+\delta A\) is nonsingular.

(c) Let \(\tilde {x}=x+\delta x\), show that \[ \lVert \tilde {x}\rVert \le \frac {1+r}{1-r}\lVert x\rVert ; \] and \[ \frac {\lVert \delta x\rVert }{\lVert x\rVert } \le \frac {2\epsilon }{1-r}\kappa (A). \]

Elliptic Equations Consider the standard five-point difference approximation (centered difference for both the gradient and divergence operators) for the variable coefficient Poisson’s equation \[ -\nabla \cdot (a\nabla u)=f \] with Dirichlet boundary conditions, in a two-dimensional rectangular region. We assume that \(a(\mathbf {x})\ge \delta >0\). The approximate solution \(\{u_{i,j}\}\) satisfies a linear system \(Au=b\).

(a) State and prove the maximum principle for the numerical solution \(u_{i,j}\).

(b) Use the maximum principle to show that the coefficient matrix \(A\) is nonsingular.

(c) Derive the matrix \(A\) in the one-dimensional case and show that it is symmetric and positive definite.

(d) With this linear system, do we get convergence with Jacobi iterations? What about Gauss-Seidel iterations? Give reasons for your answers.

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.

Sensitivity Analysis In this problem, we consider a real, non-singular \(n\times n\) matrix \(A\) and vectors \(b,x\in \mathbb {R}^n\). We are concerned with the error in solving the linear system \(Ax=b\). It turns out that the numerical solution \(\hat {x}\) is only an approximate solution to this equation. However, it satisfies a perturbed equation \(\hat {A}\hat {x}=b\) exactly, where \(\hat {A}=A+\delta A\). Let the error be defined as \(\delta x=\hat {x}-x\); we are interested in bounding \(\lVert \delta x\rVert \) in terms of \(\lVert \delta A\rVert \lVert A^{-1}\rVert \), or \(\lVert \delta A\rVert \lVert \hat {A}^{-1}\rVert \), and \(\lVert x\rVert \), or \(\lVert \hat {x}\rVert \), where \(\lVert \cdot \rVert \) is some suitable norm, and we proceed as follows:

a) First prove the identity \[ A^{-1}-\hat {A}^{-1}=A^{-1}\cdot \delta A\cdot \hat {A}^{-1}, \] and hence deduce that \[ \lVert A^{-1}-\hat {A}^{-1}\rVert \le \lVert A^{-1}\rVert \cdot \lVert \delta A\rVert \cdot \lVert \hat {A}^{-1}\rVert . \]

b) Assuming that \(\hat {\delta }=\lVert \delta A\rVert \lVert \hat {A}^{-1}\rVert <1\), show that \[ \lVert A^{-1}\rVert \le \frac {1}{1-\hat {\delta }}\lVert \hat {A}^{-1}\rVert ,\qquad \lVert A^{-1}-\hat {A}^{-1}\rVert \le \frac {\hat {\delta }}{1-\hat {\delta }}\lVert \hat {A}^{-1}\rVert . \]

c) By comparing the equations \(Ax=b\) and \((A+\delta A)(x+\delta x)=b\), show that \[ \lVert \delta x\rVert \le \frac {\hat {\delta }}{1-\hat {\delta }}\lVert \hat {x}\rVert . \] Similarly, show that if \(\delta =\lVert \delta A\rVert \lVert A^{-1}\rVert <1\), \[ \lVert \delta x\rVert \le \frac {\delta }{1-\delta }\lVert x\rVert . \]

Singular Value Decomposition The singular value decomposition (SVD) of a real \(m\times n\) matrix \(A\) (you may assume \(m\ge n\)) is defined by \[ A=U\Sigma V^T \] where \(U\) is \(m\times m\) and orthogonal. \(\Sigma \) is \(m\times n\) and diagonal. \(V\) is \(n\times n\) and orthogonal.

a) Prove that every matrix \(A\) has an SVD.

b) Use the SVD of \(A\) to prove the famous result from linear algebra that the null space of \(A\) is the orthogonal complement of the range of \(A^T\) (the adjoint of \(A\)).

c) For a nonsingular \(m\times m\) matrix \(A\), the solution to the system \(Ax=b\) can be expressed in terms of \(A\)’s SVD. Write \(b=\sum _{j=1}^m c_j u_j\) where \(u_j\) is the \(j\)th column of \(U\). Then the solution \(x\) is \[ x=\sum _{j=1}^m \frac {c_j}{\sigma _j}v_j \] where \(v_j\) is the \(j\)th column of \(V\) and \(\sigma _j\) is the \(j\)th singular value. Verify this formula. Consider the special case \(b=u_1\). What perturbation \(\delta b\) in \(b\), with \(\lVert \delta b\rVert _2=1\), causes the largest change \(\delta x\) in the solution \(x\)?

Iterative Methods for Linear Systems Consider the system \(Ax=b\), in which \(A\) is a symmetric positive definite real \(m\times m\) matrix, and denote by \(x_*\) its true solution. The conjugate gradient (CG) algorithm for finding the solution to the system is:

Set \(x_0=0\), \(r_0=b\), and \(p_0=r_0\).

For \(n=1,2,3,\ldots \), \[ \alpha _n=\frac {r_{n-1}^T r_{n-1}}{p_{n-1}^T A p_{n-1}}, \] \[ x_n=x_{n-1}+\alpha _n p_{n-1}, \] \[ r_n=r_{n-1}-\alpha _n A p_{n-1}, \] \[ \beta _n=\frac {r_n^T r_n}{r_{n-1}^T r_{n-1}}, \] \[ p_n=r_n+\beta _n p_{n-1}. \] End.

The vectors generated by this algorithm have the property that, provided \(r_{n-1}\ne 0\), \[ \begin {aligned} \mathcal {K}_n &= \operatorname {span}\{b,Ab,A^2b,\ldots ,A^{n-1}b\} = \operatorname {span}\{x_1,x_2,\ldots ,x_n\}\\ &= \operatorname {span}\{r_0,r_1,\ldots ,r_{n-1}\} = \operatorname {span}\{p_0,p_1,\ldots ,p_{n-1}\} \end {aligned} \] and \(r_n^T r_j=0\) and \(p_n^T A p_j=0\) for all \(j<n\).

Let \(e_n=\lVert x_n-x_*\rVert _A\), where the \(A\)-norm of a vector \(x\) is defined by \(\lVert x\rVert _A=(x^T Ax)^{1/2}\).

a) Prove the following: Suppose the CG algorithm is applied to the symmetric positive definite system \(Ax=b\) and that it has not yet converged (\(r_{n-1}\ne 0\)). Then \(e_n<\lVert y-x_*\rVert _A\) for all \(y\in \mathcal {K}_n\), \(y\ne x_n\); \(\lVert e_n\rVert _A\le \lVert e_{n-1}\rVert _A\); and \(e_N\to 0\) for some \(N\le m\).

b) In practice this result is not very useful. Why is this, and what determines the speed of convergence of \(x_n\) to \(x_*\)? What is ‘preconditioning’ and how is it related to the speed of convergence of the CG algorithm?

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.

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?

Consider a linear chain of \(2N\) atoms with the same spacing \(a\) and force constant \(\beta \) but of two different masses \(m,M\), \(M>m\), placed alternately. Label the light atoms by even integers \(2n\), \(n=0,\ldots ,N-1\) and the heavy atoms by odd integers \(2n-1\), \(n=1,\ldots ,N\). Denoting their displacements from equilibrium by the variables \(U_{2n}\) and \(V_{2n-1}\), Newton’s law of motion gives \[ \begin {aligned} m\ddot {U}_{2n} &= \beta [V_{2n-1}+V_{2n+1}-2U_{2n}],\\ M\ddot {V}_{2n-1} &= \beta [U_{2n}+U_{2n-2}-2V_{2n-1}]. \end {aligned} \] Determine the modes of vibration of the system and sketch the resulting pair of dispersion curves. [Assume periodic boundary conditions \(U_0=U_{2N}\) and \(V_1=V_{2N+1}\).]

Sensitivity Analysis Consider the (square, non-singular) linear system \[ Ax=b. \tag {8} \] Suppose you compute (by any method) a numerical solution \(\hat {x}\) and you are concerned about the relative error \[ E=\frac {\lVert x-\hat {x}\rVert }{\lVert x\rVert } \tag {9} \] where \(\lVert \cdot \rVert \) denotes a suitable vector norm. Derive upper and lower bounds on \(E\) in terms of the residual \(r=b-A\hat {x}\) and the condition number of \(A\) with respect to \(\lVert \cdot \rVert \).

Newton’s Method Describe how (the two-variable version of) Newton’s Method can be used to solve the nonlinear system \[ \begin {aligned} x^2+y^2 &= 1,\\ e^x-y &= 0 \end {aligned} \tag {1} \] Discuss the selection of a starting point and the advantages and disadvantages of Newton’s method.

The Jacobi Method Let \(A\) be a strictly diagonally dominant \(n\times n\) matrix, i.e., \[ \lvert a_{ii}\rvert >\sum _{j\ne i}\lvert a_{ij}\rvert , \qquad i=1,\ldots ,n. \tag {2} \] Describe the Jacobi method applied to the linear system \[ Ax=b \tag {3} \] and prove that it converges for all initial approximations of the solution.

Ill Conditioning Define the condition number of a square matrix \(A\) and its significance for the solution of the linear system \[ Ax=b. \tag {5} \] Can you define a similar number that describes the conditioning of the eigenvalue problem \[ Ax=\lambda x? \tag {6} \]

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

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\)?

Consider a linear chain of \(2N\) atoms with the same spacing \(a\) and force constant \(\beta \) but of two different masses \(m,M\), \(M>m\), placed alternately. Label the light atoms by even integers \(2n\), \(n=0,\ldots ,N-1\) and the heavy atoms by odd integers \(2n-1\), \(n=1,\ldots ,N\). Denoting their displacements from equilibrium by the variables \(U_{2n}\) and \(V_{2n-1}\), Newton’s law of motion gives \[ \begin {aligned} m\ddot {U}_{2n} &= \beta [V_{2n-1}+V_{2n+1}-2U_{2n}],\\ M\ddot {V}_{2n-1} &= \beta [U_{2n}+U_{2n-2}-2V_{2n-1}]. \end {aligned} \] Determine the modes of vibration of the system and sketch the resulting pair of dispersion curves. [Assume periodic boundary conditions \(U_0=U_{2N}\) and \(V_1=V_{2N+1}\).]

Singular Value Decomposition The singular value decomposition (SVD) of a real \(m\times n\) matrix \(A\) (you may assume \(m\ge n\)) is defined by \[ A=U\Sigma V^T \] where \(U\) is \(m\times m\) and orthogonal, \(\Sigma \) is \(m\times n\) and diagonal, \(V\) is \(n\times n\) and orthogonal.

a) Prove that every matrix \(A\) has an SVD.

b) Use the SVD of \(A\) to prove the famous result from linear algebra that the null space of \(A\) is the orthogonal complement of the range of \(A^T\) (the adjoint of \(A\)).

c) For a nonsingular \(m\times m\) matrix \(A\), the solution to the system \(Ax=b\) can be expressed in terms of \(A\)’s SVD. Write \(b=\sum _{j=1}^m c_j u_j\) where \(u_j\) is the \(j^{\text {th}}\) column of \(U\). Then the solution \(x\) is \[ x=\sum _{j=1}^m \frac {c_j}{\sigma _j}v_j \] where \(v_j\) is the \(j^{\text {th}}\) column of \(V\) and \(\sigma _j\) is the \(j^{\text {th}}\) singular value. Verify this formula. Consider the special case \(b=u_1\). What perturbation \(\delta b\) in \(b\), with \(\lVert \delta b\rVert _2=1\), causes the largest change \(\delta x\) in the solution \(x\)?

Iterative Methods for Linear Systems Consider the fixed-point iteration \[ u^{(k+1)}=Tu^{(k)}+c \] for finding a solution of the problem \[ u=Tu+c \] where \(T\) is an \(m\times m\) real matrix and \(c\) is a real \(m\)-vector.

a) Show that the fixed point iteration will converge for an arbitrary initial guess \(u^{(0)}\) if and only if the spectral radius of \(T\), \(\rho (T)\), is less than 1.

b) 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 the Jacobi iterative method for this problem. In particular, 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 of \(\delta \) depend on \(h\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

Sensitivity Analysis Consider the linear system \(Ax=b\) where \(A\) is a square, non-singular matrix. A numerical solution \(\hat {x}\) actually satisfies \[ (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_1=\frac {\lVert \delta x\rVert }{\lVert \hat {x}\rVert } \] in terms of the condition number of \(A\) and the norms of \(A,\hat {x},\delta A,\delta b\), where \(\lVert \cdot \rVert \) is a suitable vector norm. Also, we can define another relative error \[ E_2=\frac {\lVert \delta x\rVert }{\lVert x\rVert }. \] Derive an upper bound for \(E_2\) in terms of the condition number of \(A\) and the norms of \(\delta A,A,\delta b\) and \(b\). When are these two errors close to each bound? What’s the significance of the condition number?

Iterative Methods for Linear Systems Consider the linear system \(Ax=b\) and a splitting of \(A:A=M-K\), with \(M\) nonsingular, and the iteration \[ 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.

Given an \((N-1)\times (N-1)\) \[ \setlength {\arraycolsep }{5pt} A= \begin {bmatrix} 2 & -1 & 0 & \cdots & 0\\ -1 & 2 & -1 & \ddots & 0\\ 0 & -1 & \ddots & \ddots & 0\\ \vdots & \ddots & -1 & 2 & -1\\ 0 & \cdots & 0 & -1 & 2 \end {bmatrix}, \qquad b=\frac {1}{N^2} \begin {pmatrix} f_1\\ f_2\\ \cdot \\ f_{N-2}\\ f_{N-1} \end {pmatrix}. \] Analyze the convergence properties of the Jacobi and Gauss-Seidel methods for this problem. For each method, 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\) for each method?

Sensitivity Analysis Consider the (square, non-singular) linear system \[ Ax=b. \tag {3} \] Suppose you compute (by any method) a numerical solution \(\hat {x}\) and you are concerned about the relative error \[ E=\frac {\lVert x-\hat {x}\rVert }{\lVert x\rVert } \tag {4} \] where \(\lVert \cdot \rVert \) denotes a suitable vector norm. Derive upper and lower bounds on \(E\) in terms of the residual \(r=b-A\hat {x}\) and the condition number of \(A\) with respect to \(\lVert \cdot \rVert \).

Singular Value Decomposition Suppose that \(A\) is a real \(m\times n\) matrix, \(m\ge n\).

a) 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.

b) 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 \(m\times m\) 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? \]

c) For a nonsingular matrix \(A\) which does have an SVD, express the relative distance from \(A\) to the nearest rank deficient matrix in terms of \(A\)’s singular values.

Iterative Methods for Linear Systems a) Consider the fixed-point iteration \[ u^{(k+1)}=Tu^{(k)}+c \] for finding a solution of the problem \[ u=Tu+c \] where \(T\) is an \(m\times m\) real matrix and \(c\) is a real \(m\)-vector. Show that the fixed point iteration will converge for an arbitrary initial guess \(u^{(0)}\) if and only if the spectral radius of \(T\), \(\rho (T)\), is less than 1.

b) 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 the Jacobi iterative method for this problem. In particular, 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 of \(\delta \) depend on \(h\)? In practice, would you use this method to solve the given problem? If so, explain why this is a good idea? If not, how would you solve it in practice?

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? \]

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?

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

Gaussian Elimination is LU Factorization. What is the LU factorization of a square matrix? How do you solve a linear system \(Ax=b\) by Gaussian Elimination? Show that for a general \(4\times 4\) system the matrix part of Gaussian Elimination is equivalent to computing the LU factorization of \(A\).

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.

a) 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?

b) What is the Singular Value Decomposition (SVD) of \(A\)? Express the condition number above in terms of \(A\)’s singular values.

c) When the condition number is very large, the matrix is called ill-conditioned. Sometimes people say that an ill-conditioned matrix is ‘nearly-singular’. Make this precise by expressing the (relative) distance between \(A\) and the nearest singular matrix to \(A\) in terms of the condition number of \(A\).

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)\). Imagine applying the Jacobi iterative method to solve this system of linear equations. Analyze the convergence properties of the Jacobi method for this problem. In particular, express the speed of convergence as a function of the discretization stepsize \(h\). You may find it helpful to recall that the \((N-1)\times (N-1)\) matrix \(M\) defined by \[ m_{ij}= \begin {cases} 2, & \text {if } i=j,\\ -1, & \text {if } \lvert i-j\rvert =1,\\ 0, & \text {otherwise} \end {cases} \] has eigenvalues \[ \mu _l=4\sin ^2\left (\frac {\pi l}{2N}\right ),\qquad l=1,2,\ldots ,N-1. \] (NOTE: This is a model problem designed to give you a simple linear system on which to analyze the Jacobi method; one would certainly not use an iterative method to solve a tri-diagonal system in practice!)

Finite-Element Methods Consider the problem: \[ -u''+u=f,\qquad \text {on }[0,1] \] with \[ u(0)=0\qquad \text {and}\qquad u(1)=1. \]

Derive a variational formulations of this problem (be sure to define the space on which your variational problem is to be solved). Formulate a Galerkin finite-element method for the variational problem using piecewise linear functions. Prove that the linear system which arises in the finite-element method has a solution. In what sense is the solution optimal?

a) Show that Jacobi’s method for solving \(Ax=b\) iteratively is convergent when \(A\) is strictly diagonally dominant. What does this imply about the spectral radius of another matrix obtained from \(A\).

b) Describe the successive over-relaxation (SOR) method with parameter \(\omega \) and it’s rate of convergence. Give a sufficient condition upon \(A\) and \(\omega \) for convergence.

c) Briefly compare the complexity of solving Poisson’s problem in three dimensions by direct and iterative methods. Assume the simplest methods, grid, differencing scheme, boundary conditions, etc.