Chapter 11 Matrix is Close to Singular or Badly Scaled

Matlab solution to Prob. 11.11 (ii): a=[1 4 9 16;4 9 16 25;9 16 25 36;16 25 36 49] a = 1 4 9 16 4 9 16 25 9 16 25 36 16 25 36 49 b=[30 54 86 126] b = 30 54 86 126 b=b' b = 30 54 86 126 x=a\b Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.092682e-018. x = 1.1053 0.6842 1.3158 0.8947 x=inv(a)*b Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.092682e-018. x = 0 0 0 0 cond(a) ans = 4.0221e+017 11.12 Program Linsimp Use IMSL Implicit None Integer::ipath,lda,n,ldfac Parameter(ipath=1,lda=3,ldfac=3,n=3) Integer::ipvt(n),i,j Real::A(lda,lda),Rcond,Res(n) Real::Rj(n),B(n),X(n) Data A/1.0,0.5,0.3333333,0.5,0.3333333,0.25,0.3333333,0.25,0.2/ Data B/1.833333,1.083333,0.783333/ Call linsol(n,A,B,X,Rcond) Print *, 'Condition number = ', 1.0E0/Rcond Print * Print *, 'Solution:' Do I = 1,n Print *, X(i) End Do End Program Subroutine linsol(n,A,B,X,Rcond) Implicit none Integer::n, ipvt(3) Real::A(n,n), fac(n,n), Rcond, res(n) Real::B(n), X(n) Call lfcrg(n,A,3,fac,3,ipvt,Rcond) Call lfirg(n,A,3,fac,3,ipvt,B,1,X,res) End 11.13 Option Explicit Sub TestChol() Dim i As Integer, j As Integer Dim n As Integer Dim a(10, 10) As Single n = 3 a(1, 1) = 6: a(1, 2) = 15: a(1, 3) = 55 a(2, 1) = 15: a(2, 2) = 55: a(2, 3) = 225 a(3, 1) = 55: a(3, 2) = 225: a(3, 3) = 979 Call Cholesky(a(), n) 'output results to worksheet Sheets("Sheet1").Select Range("a3").Select For i = 1 To n For j = 1 To n ActiveCell.Value = a(i, j) ActiveCell.Offset(0, 1).Select Next j ActiveCell.Offset(1, -n).Select Next i Range("a3").Select End Sub Sub Cholesky(a, n) Dim i As Integer, j As Integer, k As Integer Dim sum As Single For k = 1 To n For i = 1 To k - 1 sum = 0 For j = 1 To i - 1 sum = sum + a(i, j) * a(k, j) Next j a(k, i) = (a(k, i) - sum) / a(i, i) Next i sum = 0 For j = 1 To k - 1 sum = sum + a(k, j) ^ 2 Next j a(k, k) = Sqr(a(k, k) - sum) Next k End Sub 11.14 Option Explicit Sub Gausseid() Dim n As Integer, imax As Integer, i As Integer Dim a(3, 3) As Single, b(3) As Single, x(3) As Single Dim es As Single, lambda As Single n = 3 a(1, 1) = 3: a(1, 2) = -0.1: a(1, 3) = -0.2 a(2, 1) = 0.1: a(2, 2) = 7: a(2, 3) = -0.3 a(3, 1) = 0.3: a(3, 2) = -0.2: a(3, 3) = 10 b(1) = 7.85: b(2) = -19.3: b(3) = 71.4 es = 0.1 imax = 20 lambda = 1# Call Gseid(a(), b(), n, x(), imax, es, lambda) For i = 1 To n MsgBox x(i) Next i End Sub Sub Gseid(a, b, n, x, imax, es, lambda) Dim i As Integer, j As Integer, iter As Integer, sentinel As Integer Dim dummy As Single, sum As Single, ea As Single, old As Single For i = 1 To n dummy = a(i, i) For j = 1 To n a(i, j) = a(i, j) / dummy Next j b(i) = b(i) / dummy Next i For i = 1 To n sum = b(i) For j = 1 To n If i <> j Then sum = sum - a(i, j) * x(j) Next j x(i) = sum Next i iter = 1 Do sentinel = 1 For i = 1 To n old = x(i) sum = b(i) For j = 1 To n If i <> j Then sum = sum - a(i, j) * x(j) Next j x(i) = lambda * sum + (1# - lambda) * old If sentinel = 1 And x(i) <> 0 Then ea = Abs((x(i) - old) / x(i)) * 100 If ea > es Then sentinel = 0 End If Next i iter = iter + 1 If sentinel = 1 Or iter >= imax Then Exit Do Loop End Sub 11.15 As shown, there are 4 roots, one in each quadrant. 8 g (−0.618,3.236) 4 (1, 2) 0 -4 -2 0 -4 2 (1.618, −1.236) f (−2, −4) -8 It might be expected that if an initial guess was within a quadrant, the result wouls be the root in the quadrant. However a sample of initial guesses spanning the range yield the following roots: 6 3 0 -3 -6 (-2, -4) (-0.618,3.236) (-0.618,3.236) (1,2) (-0.618,3.236) (-0.618,3.236) (-0.618,3.236) (1,2) (1,2) (1.618, -1.236) (1.618, -1.236) (1.618, -1.236) (-2, -4) (-2, -4) (1.618, -1.236) (1.618, -1.236) (-2, -4) (-2, -4) (-2, -4) (1.618, -1.236) -6 -3 0 3 (-0.618,3.236) (-0.618,3.236) (1.618, -1.236) (1.618, -1.236) (-2, -4) 6 We have highlighted the guesses that converge to the roots in their quadrants. Although some follow the pattern, others jump to roots that are far away. For example, the guess of (-6, 0) jumps to the root in the first quadrant. This underscores the notion that root location techniques are highly sensitive to initial guesses and that open methods like the Solver can locate roots that are not in the vicinity of the initial guesses. 11.16 x = transistors y = resistors z = computer chips System equations: 3 x + 3 y + 2 z = 810 x + 2 y + z = 410 2 x + y + 2 z = 490 3 3 2    Let A = 1 2 1  and B = 2 1 2  810  410    490  Plug into Excel and use two functions- minverse mmult Apply Ax = B x = A-1 * B Answer: x = 100, y = 110, z = 90 11.17 As ordered, none of the sets will converge. However, if Set 1 and 3 are reordered so that they are diagonally dominant, they will converge on the solution of (1, 1, 1). Set 1: 8x + 3y + z = 12 2x + 4y – z = 5 −6x +7z = 1 Set 3: 3x + y − z = 3 x + 4y – z = 4 x + y +5z =7 At face value, because it is not diagonally dominant, Set 2 would seem to be divergent. However, since it is close to being diagonally dominant, a solution can be obtained by the following ordering: Set 3: −2x + 2y − 3z = −3 2y – z = 1 −x + 4y + 5z = 8 11.18 Option Explicit Sub TriDiag() Dim i As Integer, n As Integer Dim e(10) As Single, f(10) As Single, g(10) As Single Dim r(10) As Single, x(10) As Single n = 4 e(2) = -1.2: e(3) = -1.2: e(4) = -1.2 f(1) = 2.04: f(2) = 2.04: f(3) = 2.04: f(4) = 2.04 g(1) = -1: g(2) = -1: g(3) = -1 r(1) = 40.8: r(2) = 0.8: r(3) = 0.8: r(4) = 200.8 Call Thomas(e(), f(), g(), r(), n, x()) For i = 1 To n MsgBox x(i) Next i End Sub Sub Thomas(e, f, g, r, n, x) Call Decomp(e, f, g, n) Call Substitute(e, f, g, r, n, x) End Sub Sub Decomp(e, f, g, n) Dim k As Integer For k = 2 To n e(k) = e(k) / f(k - 1) f(k) = f(k) - e(k) * g(k - 1) Next k End Sub Sub Substitute(e, f, g, r, n, x) Dim k As Integer For k = 2 To n r(k) = r(k) - e(k) * r(k - 1) Next k x(n) = r(n) / f(n) For k = n - 1 To 1 Step -1 x(k) = (r(k) - g(k) * x(k + 1)) / f(k) Next k End Sub 11.19 The multiplies and divides are noted below Sub Decomp(e, f, g, n) Dim k As Integer For k = 2 To n e(k) = e(k) / f(k - 1) f(k) = f(k) - e(k) * g(k - 1) Next k End Sub '(n – 1) '(n – 1) Sub Substitute(e, f, g, r, n, x) Dim k As Integer For k = 2 To n r(k) = r(k) - e(k) * r(k - 1) '(n – 1) Next k x(n) = r(n) / f(n) ' 1 For k = n - 1 To 1 Step -1 x(k) = (r(k) - g(k) * x(k + 1)) / f(k) '2(n – 1) Next k End Sub Sum = 5(n-1) + 1 They can be summed to yield 5(n – 1) + 1 as opposed to n3/3 for naive Gauss elimination. Therefore, a tridiagonal solver is well worth using. 1000000 100000 Tridiagonal Naive Gauss 10000 1000 100 10 1 1 10 100