CS111 Homework 3 Solution

1. Consider the code for the temperature problem in the file cs111/temperature.py, especially the routine make b() that creates the right-hand side b. Let k = 100.
1.1. How many elements are there in b?
1.2. Considering all possible choices for the temperatures on the boundary, what is the smallest number of elements of b that could possibly be nonzero? Explain why in English.
1.3. Considering all possible choices for the temperatures on the boundary, what is the largest number of elements of b that could possibly be nonzero? Explain why in English.
2. The temperature problem models our cabin in the woods in two dimensions, but most modern scientific simulations are done in three dimensions. Here you will create the matrix that corresponds to a 3-D version of the temperature problem. The “cabin” is now the unit cube. As before, we will discretize the interior by dividing it into k points in each dimension, but now there are a total of k3 points rather than k2.
The partial differential equation still leads to the approximation that the temperature at any given point is the average of the temperatures at the neighboring points, but now there are 6 neighbors, with 2 in each of the three dimensions. The matrix A for the 3D version of the temperature problem expresses the fact that, in a 3D k-by-k-by-k grid, each interior point has a temperature that is the average of its 6 neighbors (left, right, up, down, in, out). The diagonal elements of A are all equal to 6, and the off-diagonal elements are either 0 or −1. Most of the rows of A have 7 nonzeros.
Using make A(k) from cs111/temperature.py as a model, write a routine make A 3D(k) that returns the k3-by-k3 matrix A. Like make A(k), your routine should create A as a sparse matrix using scipy.sparse.csr matrix(), after making a table “triples” of the nonzero elements of A.
Here below, for your debugging, is the correct matrix for k = 2. I converted it to a dense array for printing—you should also print it out as a sparse matrix (that line is commented out below), and indeed for k > 2 it’s going to be too large to see what’s going on in the dense matrix anyway.
[In:]
k = 2
A = make_A_3D(k) print(’k:’, k) print(’dimensions:’, A.shape) print(’nonzeros:’, A.nnz) #print(’A as sparse matrix: ’, A) print(’A as dense matrix: ’, A.toarray())
[Out:] k: 2
dimensions: (8, 8) nonzeros: 32
A as dense matrix:
[[ 6. -1. -1. 0. -1. 0. 0. 0.]
[-1. 6. 0. -1. 0. -1. 0. 0.]
[-1. 0. 6. -1. 0. 0. -1. 0.]
[ 0. -1. -1. 6. 0. 0. 0. -1.]
[-1. 0. 0. 0. 6. -1. -1. 0.]
[ 0. -1. 0. 0. -1. 6. 0. -1.]
[ 0. 0. -1. 0. -1. 0. 6. -1.]
[ 0. 0. 0. -1. 0. -1. -1. 6.]]
To complete a realistic simulation you would also write a routine make b 3D(k) to compute the right-hand side b. For this problem, you don’t have to do that.
3. Do problem 2.3 on pages 32–33 of Chapter 2 of the NCM book, showing the numpy code you use and its output. Note: To understand intuitively what the problem means by “assume that joint 1 is rigidly fixed both horizontally and vertically and that joint 8 is fixed vertically,” think of the truss as a (2-dimensional) drawbridge across a river, with the left end being a hinge and the right end lying on the ground.
4. The condition number of a matrix A is a measure of how close to being singular it is. The definition is κ(A) = ||A|| ||A−1||,
with the convention κ(A) = ∞ if A is singular. Each different matrix norm gives its own definition of condition number; usually we’ll consider the 2-norm condition number. The condition number is always at least 1. (I’ll prove that in class.) A matrix is well-conditioned if its condition number is not too big (think 103 or less), and ill-conditioned if its condition number is very large (think 106 or more). In numpy, the condition number is computed by npla.cond().
In this problem you’ll explore the relationship between condition number and the behavior of LU factorization.
4.1. Consider the linear system
,
for some α < 1. Clearly the solution is (x0,x1)T = (−1,1)T.
For each value of α = 10−4, 10−8, 10−16, 10−20, try to solve this system twice, using the routine cs111.LUsolve(), first with pivoting = True and then with pivoting = False. For each solution, print: the condition number of A; the computed x; the norm of the error; and the relative residual norm returned by cs111.LUsolve(). Show your numpy code and its output.
As α decreases, does A become ill-conditioned or does it remain well-conditioned? Describe and explain the trends in the relative residual norm for both the non-pivoting and pivoting solutions. 4.2. Consider the linear system
,
for some α < 1. Again, the solution is (x0,x1)T = (−1,1)T.
For each value of α = 10−4, 10−8, 10−16, 10−20, try to solve this system twice, using the routine cs111.LUsolve(), first with pivoting = True and then with pivoting = False. This time, some of the tries will fail to return an answer. For each try, print the condition number of A. If an answer was returned, print: the computed x; the norm of the error; and the relative residual norm returned by cs111.LUsolve(). Show your numpy code and its output.
As α decreases, does A become ill-conditioned or does it remain well-conditioned? Describe and explain the trends in the relative residual norm for the successful tries. Explain why the other tries were unsuccessful.