==================================== CHAPTER 2 PROBLEMS GERALD AND WHEATLEY, 6TH EDITION ==================================== ------------------ 17. The final augmented matrix is: [ 1 0 0 | 22/9 ] [ 0 1 0 | 3 ] [ 0 0 1 | 11/9 ] ------------------ 22. Here are the factors. You may either use eqns 2.20 and 2.21 to get these, or use the 'slicing' algorithm that I showed in class, which is demonstrated in example 2.3 on page 139 of the textbook. [1 0 0][1 1 -2 ] [4 -6 0][0 1 -1.5] [3 -4 3][0 0 1 ] L U As the book says, you first solve for Ly=b to get y=[3 7/6 11/9], then solve Ux=y to get [22/9 3 11/9]. To solve the system for a NEW b vector (right-hand side), you just repeat the Ly=b and Ux=y solves, using the new b. Both of these are cheap solves because the matrices involved are triangular. I get y=[2 3/2 5/3] for the Ly=b solve and then x=[4/3 4 5/3] from the Ux=y solve. ------------------ 31. a. Ill-conditioned matrices are "close to" singular, in the sense that their determinants are tiny. (A singular matrix has a zero determinant.) If you plug an ill-conditioned matrix into an Ax=b solver, the answer will be very sensitive to numerical error. I used Mathematica to compute the determinant of H_4: In[8]:= H4 = {{1, 1/2, 1/3, 1/4},{1/2, 1/3, 1/4, 1/5},{1/3, 1/4, 1/5, 1/6}, {1/4, 1/5, 1/6, 1/7}} 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Out[8]= {{1, -, -, -}, {-, -, -, -}, {-, -, -, -}, {-, -, -, -}} 2 3 4 2 3 4 5 3 4 5 6 4 5 6 7 In[9]:= Det[H4] 1 Out[9]= ------- 6048000 This is tiny, so the matrix is ill conditioned. b. I got [1.11 0.228 1.95 0.797] for this. The right answer is [1 1 1 1], so this answer is pretty terrible. The reason is that the matrix is ill conditioned; see part a, above. The mechanism is that the matrix's ticklish nature magnifies the effects of the finite precision arithmetic. c. I got [0.988 1.42 -0.428 2.10] for this. Note how different this is from the answer in part b. This, again, is because of the near-singular nature of the matrix, which magnifies the error introduced by the rounding to three digits. Since rounding error is different from chopping error, the part b and part c answers are very different. ------------------ 32. You should get the same answer --- 232 --- for both parts (a) and (b). Note that L only has 1s on the diagonal if you use the formal Gaussian elimination algorithm (subtracting a_i1/a_11 times the first equation from the ith equation, and so on). If you use some other variant --- say, multiplying one row by 2 and the other by -5 and adding them --- the corresponding value on the diagonal of L is NOT 1! c. To compute the inverse of this matrix, you augment it with the identity: [2 5 -1 | 1 0 0] [1 6 4 | 0 1 0] [7 -4 2 | 0 0 1] ...and perform Gauss-Jordan, obtaining: [1 0 0 | 0.1207 -0.0258 0.1121] [0 1 0 | 0.1121 0.0474 -0.0388] [0 0 1 | -0.1983 0.1853 0.0302] Here's how to do this in Mathematica: In[4]:= B = {{2,5,-1,1,0,0},{1,6,4,0,1,0},{7,-4,2,0,0,1}} Out[4]= {{2, 5, -1, 1, 0, 0}, {1, 6, 4, 0, 1, 0}, {7, -4, 2, 0, 0, 1}} In[5]:= RowReduce[B] 7 3 13 13 11 9 Out[5]= {{1, 0, 0, --, -(---), ---}, {0, 1, 0, ---, ---, -(---)}, 58 116 116 116 232 232 23 43 7 > {0, 0, 1, -(---), ---, ---}} 116 232 232 ------------------ 39. a. 1-norm = 17.45 2-norm = 10.912 infinity-norm = 10.0 c. 1-norm = 22 2-norm = 18.841 infinity-norm = 23 ------------------ 46. Condition(A) = 55,228 for the problem 42 matrix This is definitely an ill-conditioned matrix! ------------------ 50. If I solve problem 42 using 10-digit accuracy, I get x = (1592.61, -631.911, -493.62) If I use 3 digits, though, I get x = (-118, 47.1, 37.0). This is very different from the 10-digit answer because the matrix is ill conditioned and thus acts as a magnifier of any numerical approximation. The iterative improvement procedure is shown on pp159-160 of the text. I'll treat the 10-digit x as the "right answer," and try to improve on the 3-digit answer: x-bar = (-118, 47.1, 37.0) Multiply A by x-bar and subtract that from b to get r = (-1.463, 0.434, -1.5633) Then use your Ax=b solver to solve the system Ae=r, which gives you an approximation to the error: e-bar = (1710, -679, -529) ..... using 6 digit precision And add that into the old x-bar to get the improved answer: improved x = (1562, -632, -492) This is a much better answer. You can repeat this process over and over again to improve matters further. Note that if you use 3-digit precision to compute e-bar, you don't get as much improvement. ------------------ 55. If you just blindly solve the first eqn for the first variable, the second eqn for the second variable, and so on, you run into trouble in this example because the third eqn does not contain the third variable. This means that you have to swap them around first. It's a good idea to try to get the |biggest| coefficients on the diagonal (a.k.a. "diagonal dominance"), so here's how I'd order the eqns: 9x + 4y + z = -17 x + 6y = 4 x - 2y - 6z = 14 Then the rearrangement goes as follows: x = (-17 - 4y - z)/9 y = (4 - x)/6 z = (14 - x + 2y)/6 The corresponding Jacobi iteration sequence is: x_{n+1} = (-17 - 4y_n - z_n)/9 y_{n+1} = (4 - x_n)/6 z_{n+1} = (14 - x_n + 2y_n)/6 The corresponding Gauss-Seidel iteration sequence is: x_{n+1} = (-17 - 4y_n - z_n)/9 y_{n+1} = (4 - x_{n+1})/6 z_{n+1} = (14 - x_{n+1} + 2y_{n+1})/6 If I start Jacobi from (0,0,0), the next guess is (-17/9, 2/3, -14/6). Here's a Scheme function that implements this iteration sequence: (define (prob55-jacobi x y z n) (if (> 1 n) '() (let* ((new-x (/ (- -17 (* 4 y) z) 9.0)) (new-y (/ (- 4 x) 6.0)) (new-z (/ (- 14 x (* -2 y)) -6.0))) (cons (list new-x new-y new-z) (prob55-jacobi new-x new-y new-z (-1+ n)))))) And here are the first ten iterates: (pp (prob55-jacobi 0 0 0 10)) ((-1.8888888888888888 .6666666666666666 -2.3333333333333335) (-1.9259259259259258 .9814814814814815 -2.8703703703703702) (-2.006172839506173 .9876543209876543 -2.9814814814814814) (-1.996570644718793 1.001028806584362 -2.9969135802469133) (-2.0008001828989483 .9994284407864655 -2.999771376314586) (-1.9997713763145863 1.0001333638164913 -2.9999428440786464) (-2.000065623465258 .9999618960524311 -3.000006350657928) (-1.999982359283533 1.0000109372442096 -2.9999982359283535) (-2.000005057005387 .9999970598805888 -3.000000705628659) (-1.9999986148770772 1.000000842834231 -2.9999998627944273)) This is headed towards (-2,1,-3), which is the solution to that set of equations. Adapting this for Gauss-Seidel involves changing some of the local variable names: (define (prob55-gauss-seidel x y z n) (if (> 1 n) '() (let* ((new-x (/ (- -17 (* 4 y) z) 9.0)) (new-y (/ (- 4 new-x) 6.0)) (new-z (/ (- 14 new-x (* -2 new-y)) -6.0))) (cons (list new-x new-y new-z) (prob55-gauss-seidel new-x new-y new-z (-1+ n)))))) Note how much faster this converges: (pp (prob55-gauss-seidel 0 0 0 10)) ((-1.8888888888888888 .9814814814814815 -2.9753086419753085) (-1.9945130315500685 .9990855052583448 -2.998780673677793) (-1.999729038595065 .9999548397658442 -2.999939786354459) (-1.9999866191898796 .9999977698649799 -2.9999970264866396) (-1.9999993392192532 .9999998898698755 -2.9999998531598338) (-1.9999999673688522 .9999999945614754 -2.999999992748634) (-1.9999999983885854 .9999999997314308 -2.9999999996419078) (-1.9999999999204239 .9999999999867373 -2.9999999999823164) (-1.9999999999960703 .9999999999993451 -2.9999999999991265) (-1.9999999999998057 .9999999999999676 -2.999999999999957)) ------------------ 57. Same system as above, but a different rearrangement. 9x + 4y + z = -17 x + 6y = 4 x - 2y - 6z = 14 To rearrange for relaxation, put all terms of each eqn on the LHS and then divide through by the negative of the largest coeff: (9x + 4y + z + 17)/-9 = R_1 (x + 6y - 4)/-6 = R_2 (x - 2y - 6z - 14)/6 = R_3 This yields -x - 4/9y - z/9 - 17/9 = R_1 -x/6 - y + 2/3 = R_2 x/6 - y/3 - z - 7/3 = R_3 Plugging in (0,0,0), I get R_1 = -1.89, R_2 = -.67, and R_3 = 2.33. The biggest one is R_3, so I go back to the third equation and change z so as to zero out that residual: 0/6 - 0/3 - z - 7/3 = 0 so z = -7/3, and the next guess is (0,0,-7/3), and the corresponding residuals are R_1 = -1.63, R_2 = -.67, and R_3 = 0. The biggest residual is now R_1, so I go back to the first equation and use x to zero it out: -x - 4/9(0) - (-7/3)/9 - 17/9 = R_1 .....and so on. This too converges to (-2, 1 -3). ------------------ 62. The Jacobian is [ yz-2x xz+2y xy ] [ y x -2z ] [ e^x -e^y 1 ] For the starting guess given in the problem --- (1,1,1) --- eqn 2.42 becomes: [-1 3 1][delta-x] [-0.34] [ 1 1 -2][delta-y] = - [-0.09] [ e -e 1][delta-z] [ 0.59] I used Mathematica to solve this for the deltas: In[1]:= A = {{-1,3,1,0.34},{1,1,-2,0.09},{2.72,-2.72,1,-0.59}} In[2]:= RowReduce[A] Out[2]= {{1, 0, 0, -0.103622}, {0, 1, 0, 0.0951969}, {0, 0, 1, -0.0492126}} So the next guess is (1 - 0.103622, 1 + 0.0951969, 1 - 0.0492126) or (0.896378, 1.0951969, 0.9507874). For that (x,y,z) value, eqn 2.42 becomes [-0.751457 3.042659 0.981710 ][delta-x] [-1.0639e-2] [ 1.095197 0.896378 -1.901574 ][delta-y] = - [-1.2286e-2] [ 2.450711 -2.989771 1 ][delta-z] [1.72663e-3] ...which has the solution (0.00585362 0.00515511, -0.000659562), so the next guess is (.90223162, 1.10035201, .950127838), and so on. Note that the deltas get smaller at every step; this is a sign that the method is converging -- in this case, to (0.90223, 1.10035, 0.95013). (You can save some time by only re-evaluating the Jacobian every n^th step, by the way, where n is chosen "heuristically".) ------------------ 81. If I had n processors to throw at this problem, I would probably allocate one to each Jacobi equation and compute each element (x,y,z) of the new guess in parallel. To use n^2 processors, I would assign each one to one of the *terms* in the Jacobi iteration sequence (the 4 y_n etc. above).