Due: 11:59:59pm Friday October 17th , 2014 AMATH 301 - Autumn 2014 Homework 2 VERY IMPORTANT: Remember to suppress all output and plotting before submitting your code. Exercise 1: Compute the eigenvalues of the following matrices by hand and store them in Matlab as a column vector [|λ1 |, |λ2 |]T for each matrix so that |λ1 | ≤ |λ2 |: 2 0 2 1 0 8 cos θ sin θ (a) (b) (c) (d) 0 4 2 1 2 0 − sin θ cos θ Answers: Should be written out as A1-A4.dat. Consider the following connectivity matrix and it’s graphical representation:     A=    0 1 1 0 0 0 1 0 1 0 0 0 1 1 0 1 1 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0         Note that aij is either 0 or 1 depending on if there is a connection between location i and location j. Use Matlab to get a diagonal matrix of the eigenvalues of A (call it D) and a matrix of eigenvectors for A (call it V). Use D to find the smallest magnitude eigenvalue (λ such that |λ| is the smallest). Answer: Should be written out as A5.dat. Use V to find the normalized eigenvector that corresponds to the smallest magnitude eigenvalue (keep the vector as a 6x1). Answer: Should be written out as A6.dat. A special property of connectivity matrices is that A2 shows which locations are 2 connections apart. This generalizes to An showing which locations are n connections apartment. Recall that because A = VDV−1 , powers of A can be efficiently computed by An = (VDV−1 )n = VDn V−1 . Use this to compute A3 to verify that Seattle and Charlotte are 3 connections apart. Answer: Should be written out as A7.dat. Also recall that Dn is simply a diagonal matrix with λn on the diagonal. Looking at the eigenvalues in D, what do you expect to happen to the magnitude of the values in An as n → ∞. Use 1 for “grow to infinity”, 2 for “shrink to zero”, and 3 for “neither grow nor shrink signficantly”. Answer: Should be written out as A8.dat. Think about what implications this would have if you were iterating with xk+1 = Axk . 1 AMATH 301 - Autumn 2014 Homework 2 Due: 11:59:59pm Friday October 17th , 2014 Exercise 2: The forces in the beams of a simple bridge can be modeled as Ax = b, where x = [F1 . . . F9 ]T ,        A=         0 0 −1 0 0 0 0 0 0  0 0 0 − sin(π/4) −1 0 0 0 0     0 sin(π/4) 1 sin(π/4) 0 0 0 0 0      0 0 0 0 0 1 sin(π/4) 0 0    1 − cos(π/4) 0 0 0 0 0 0 0 ,b =    1 0 0 0 0 − cos(π/4) 0 0 0      0 0 0 cos(π/4) 0 0 0 1 −1     0 0 0 0 0 0 1 −1 0 0 cos(π/4) 0 0 0 cos(π/4) 1 0 1 w1 w2 0 0 0 0 0 0 0               First, verify that A is not singular by computing the condition number. Answer: Should be written out as A9.dat. Now that we know there’s a solution, find a LU decomposition for A using Matlab. Using the permutation matrix P, determine if any rows of A were interchanged by Matlab to form the decomposition (use 1 for yes, 0 for no). Answer: Should be written out as A10.dat. The numbers w1 and w2 represent the mass of the cars at the locations 1 and 2 on the bridge. Solve for the forces (x), given w1 = 10 and w2 = 10, by first solving Ly = b and then Ux = y. Answer: Write out the column vector x as A11.dat . Let’s say that the beams on this bridge can only withstand forces up to 30 (compression or extension... positive or negative). Assuming that the mass at each location is equal (w1 = w2 ), find the minimum value w1 that will collapse this bridge to the nearest integer (note you are still using the same L, U, and P from before). You can use the infinity norm to look at the maximum force in the solution and automate this searching process with a loop (for or while). Answer: Should be written out as A12.dat. 2 AMATH 301 - Autumn 2014 Due: 11:59:59pm Friday October 17th , 2014 Homework 2 Now consider adding weight to location w1 (keep w2 = 10). Again find the minimum mass at w1 at which the bridge collapses. Answer: Should be written out as A13.dat. Finally, consider adding weight to location w2 while keeping w1 = 10. Find the minimum mass at w2 at which the bridge collapses. Answer: Should be written out as A14.dat. Although you can probably guess from the drawing, but what does this tell you about a weak point in the bridge? 3 AMATH 301 - Autumn 2014 Homework 2 Due: 11:59:59pm Friday October 17th , 2014 Exercise 3: Steady fluid flow in a pipe can modeled as Ax = b, where x = [x1 . . . x100 ]T represents the velocity at 100 spatial locations across the pipe,  A100x100      =     1 0 0 0 0 ... 0 1 −2 1 0 0 ... 0 0 1 −2 1 0 ... 0 .. . . . .. .. .. .. . . . . . .. . 0 ... 0 1 −2 1 0 0 ... 0 0 1 −2 1 0 0 0 0 0 0 1   0 −0.01 −0.01 −0.01 .. .            , b100x1 =          −0.01 0            While there a many ways to construct this matrix, one way is to start by constructing a 100 by 100 tri-diagonal matrix with 1 on the sub-diagonals and −2 on the main diagonal (use of Matlab’s diag function is recommended). Modify the first and last row of the tri-diagonal matrix to get A. Now try computing a LU decomposition and use the commands tic and toc to time the decomposition (tic; [L,U,P] = lu(A); toc). Note how much longer this takes than the bridge problem. Due to the nature of A (large and sparse), it is preferable to use iterative methods here over LU decomposition. First, let’s see how well Jacobi iteration handles this system. Note that the Jacobi interation can be written as xk+1 = Mxk + c. Check to make sure that this method will converge by computing the spectral radius ρ of M (ρ(M) = max |λj |). Answer: Should be written out as A15.dat. Solve for x using Jacobi iteration to a tolerance of rk = |xk+1 − xk |∞ < 1E-5. Use an initial guess of all zeros. Plot the values of x to verify that your iteration produced at least a smooth result. You also should plot the values of rk to ensure that the residual is going to zero. Save only the first hundred rk values as a 1x100 vector, and save your iteration count. Answers: Should be written out as A16-A17.dat. Now, let’s see how well Gauss-Seidel iteration handles the system. Again, Gauss-Seidel can be written as xk+1 = Mxk + c, with a different M and c than Jacobi iteration. Check to make sure that Gauss-Seidel iteration will converge by computing the spectral radius of M. Answer: Should be written out as A18.dat. Solve using Gauss-Seidel iteration to a tolerance of rk < 1E-5 and an initial guess of all zeros. Again, plot both x and rk to verify that your code is working (the solution x should be very similar to the Jacobi solution). Save only the first hundred rk values as a 1x100 vector, and save your iteration count. Answers: Should be written out as A19-A20.dat. Recall from class that Gauss-Seidel is considered an improvement to Jacobi. Do you see this in your spectral radii, rk plots, and iteration counts? 4