Mathematics 172 – Numerical Analysis 2 Lab/Problem Set 5 April 8, 1999

Goals

Let A = −L + D − U be the decomposition of the matrix A into its diagonal, and negative upper- and lower-triangular parts. Today’s lab will illustrate the comparisons between the Jacobi:

x^{(n)} = D^{−1}(L + U)x^{(n−1)} + D^{−1}b

and Gauss-Seidel:

x^{(n)} = (D − L)^{−1}Ux^{(n−1)} + (D − L)^{−1}b

methods we talked about in class on Tuesday. We will also investigate the possibility of attaining even faster convergence for iterative methods for Ax = b, by considering the “parametrized family” of methods containing Gauss-Seidel:

(1)  x^{(n)} = (D − ωL)^{−1}((1 − ω)D + ωU)x^{(n−1)} + (D − ωL)^{−1}ωb

(Gauss-Seidel is the case ω = 1). These methods go by the name of SOR (“successive over-relaxation”) methods.

Further Features of the linalg package

Maple’s linalg package contains a number of functions for matrix arithmetic operations, including the following operations that will be useful today:

1) multiply – matrix multiplication

2) scalarmul – multiply a matrix times a scalar

3) inverse – compute the inverse of a square matrix (This is done via the usual method: to find A^{−1}, we reduce the n × 2n matrix (A|I) to echelon form using Gauss-Jordan elimination. In the resulting matrix (I|B), the second block is B = A^{−1}.)

4) norm – compute the matrix norm (Maple knows the 1-, 2- , and ∞ norms, as well as the Frobenius norm as in the exercise from Problem Set 4.

5) eigenvals – compute the eigenvalues of a matrix (real and complex). We will learn more about the numerical methods Maple uses for problems like this in the last unit of the course.

Look at the on-line help page for the linalg package for information about the syntax of these commands. If you type in the entries of a matrix as real numbers in Maple, all arithmetic will be done in floating point. Be sure you do that in all calculations today.

Important Note: Maple has special evaluation rules for matrices. This means that if you assign a array result to a variable A, then tell Maple to print out the value of that variable, you will just see the variable name as output, not the matrix you want. To see the matrix,

use eval(A);

Here are sample commands for a Maple session illustrating their use. Suppose we want to compare the rates of convergence for Jacobi and Gauss-Seidel iteration for the system with A given by

A = ( \begin{matrix}6 & −2 & −1 \\ −2 & 7 & −2 \\ −1 & −2 & 8 \end{matrix} )

We begin by defining the L, D, U matrices. First,

                DD:=array([[6.0,0,0],[0,7.0,0],[0,0,8.0]]);

(don’t use D – that’s reserved for the differentiation operator for functions in Maple!), and similarly for L and U. Then, for example, using the commands above, we can compute the Jacobi G-matrix as follows:

                    GJac := multiply(inverse(DD),L+U);

The 1-norm of this matrix can be computed as follows:

                               norm(GJac,1);

You can also get the more refined information of the spectral radius ρ(G_{Jac}) from the output of:

                             eigenvals(GJac);

(How?) The Gauss-Seidel G-matrix is

                     GGS := multiply(inverse(DD-L),U);

You can compute its matrix norms and spectral radius for G_{GS} in the same way you did for G_{Jac} above. How do they compare to the norms and spectral radius for G_{Jac} for this matrix?

Lab Questions

A) Consider the matrix:

A = ( \begin{matrix}10 & −2 & −1 & 0 \\ −2 & 10 & −5 & 0 \\ −1 & −5 & 8 & −2 \\ 0 & 0 & −2 & 5 \end{matrix} )

1) Note that A is a symmetric matrix, so by the spectral theorem from linear algebra, all of its eigenvalues should be real. Check that in fact all of them are strictly positive. We say a symmetric matrix A is positive-definite if all of its eigenvalues are > 0.

2) Compute the Jacobi and Gauss-Seidel G-matrices for this A and compare their 1, 2, ∞-norms and spectral radii. Which iterative method will converge faster for a system Ax = b for this coefficient matrix A? (Note that the theorem we quoted in class comparing the Jacobi and Gauss-Seidel methods does apply here – this result is also Theorem 7.22 in the text.)

3) Find a 4 × 4 matrix A′ with a′_{ii} > 0 all i, a′_{ij} ≤ 0 if i ≠ j for which 1 < ρ(G′_{Jac}) < ρ(G′_{GS}). (This is the second case in Theorem 7.22 from the text.)

4) Now consider the one-parameter family of matrices for SOR iteration:

G_{ω} = (D − ωL)^{−1}((1 − ω)D + ωU)

Using Maple, determine the interval of ω-values for which ∥G_{ω}∥_{∞} < 1, and estimate the value of ω that minimizes ∥G_{ω}∥_{∞}. (Suggestion: You can plot ∥G_{ω}∥_{∞} as a function of ω in Maple(!).)

5) Now, try to estimate the value of ω that minimizes the spectral radius of G_{ω}. (It is harder, though certainly not impossible, to generate a plot in this case; if you don’t see how to do that, either “trial and error” or some kind of systematic search will be necessary.) Is this the same ω as in part 3? Explain.

6) Perform 4 iterations of Gauss-Seidel and SOR (see equation (1)) above) using your optimal ω from part 4 to solve the system Ax = b, where b = (1, 0, 1, 0)^{t}. Discuss your results.

Comment: As you might guess, the choice of the optimal ω for SOR iteration is something of a “black art”. There is one case, though, in which the optimal ω can be determined exactly. If A is positive definite and tridiagonal, then it can be proved that the ω minimizing the spectral radius of G_{ω} is

ω_{opt} = \frac{2}{1 + \sqrt{1 − (ρ(G_{Jac}))^{2}}}

where ρ(G_{Jac}) is the spectral radius of the corresponding Jacobi G-matrix. This information may be useful for the next problem!

B) Problem 12 from Section 7.3 in the text. Use Gauss-Seidel and/or SOR iteration to find approximate solutions accurate to a tolerance of 10^{−5} in the ∞-norm. Of course, one of the main questions you will need to address is: When do you want to stop iterating? As part of your solution explain the stopping criteria you used. Implementation note: You DON’T want to enter a 50 × 50 or 100 × 100 matrix entry by entry(!) Since most of the entries are zero anyway, it will be simpler to set up a matrix of all zeroes, then fill in the non-zero ones to get A (or L, D, U). The initialization:

                  B := array(1..50,1..50,sparse);

sets up a 50 × 50 matrix of all zeroes. Then you can change an entry to something else when necessary. B[i,j] is the Maple notation for the entry in row i and column j in B. For instance, the command B[3,4] := .4; would set the entry in row 3 and column 4 in B to .4. Even now, think about how you can enter all the entries for L, D, U before setting them manually one by one!!

Assignment

Individual write-ups due in class on April 15.