SolnLinearSys.doc

PAGE PAGE 2 Numerical Solution of Linear Systems Elimination Methods As an example consider the system EMBED Equation.3 or EMBED Equation.3 . In order to solve we can: Multiply a row by a number, Add a row to another row, and Interchange two rows. First construct the augmented matrix EMBED Equation.3 . Next make new EMBED Equation.3 zero by performing EMBED Equation.3 and then make the new EMBED Equation.3 zero by , then EMBED Equation.3 . Now make new EMBED Equation.3 zero by EMBED Equation.3 then EMBED Equation.3 . Now make new EMBED Equation.3 zero by EMBED Equation.3 and then make new EMBED Equation.3 zero by EMBED Equation.3 then EMBED Equation.3 . Now make EMBED Equation.3 zero by EMBED Equation.3 then EMBED Equation.3 . Finish by EMBED Equation.3 then EMBED Equation.3 . This is EMBED Equation.3 . Backward Substitution From the point where we had an upper triangular matrix we could have solved by moving up the matrix, that is EMBED Equation.3 . Gauss and Gauss-Jordon Elimination Methods Form the augmented matrix EMBED Equation.3 nx(n+1) Set EMBED Equation.3 Interchange rows to make EMBED Equation.3 the largest for rows i, i+1,?n Create zeros under row i by subtracting EMBED Equation.3 from row j for j=i+1, i+2,?, n Set i = i+1 If i < n then got to step 3 Solve remaining upper triangular matrix by backward substitution. That is Set EMBED Equation.3 For i=n-1 to 1 with steps -1 do EMBED Equation.3 . Iterative Methods These are excellent methods for parallel computing. For example, consider EMBED Equation.3 . Note that the diagonal elements dominate (i.e., for each equation EMBED Equation.3 ). Then we could use the idea from fixed point iteration for finding roots of an equation, or EMBED Equation.3 . That is if EMBED Equation.3 then let EMBED Equation.3 . Where EMBED Equation.3 has zero elements off the diagonal and contains the diagonal elements of EMBED Equation.3 on its diagonal while EMBED Equation.3 has zero elements on the diagonal and contains the other elements of EMBED Equation.3 off the diagonal. Then EMBED Equation.3 or EMBED Equation.3 and EMBED Equation.3 . Then set EMBED Equation.3 and next iterate using EMBED Equation.3 which will converge if EMBED Equation.3 . This method is sometimes called Jacobi iteration. If a scalar processor (not vector or parallel) then it would help to substitute the new EMBED Equation.3 immediately. This is Gauss-Sidel iteration. LU Decomposition If we could decompose EMBED Equation.3 in EMBED Equation.3 (Bars and arrows have been dropped for the remaining discussion) so that EMBED Equation.3 where EMBED Equation.3 is lower triangular and EMBED Equation.3 is upper triangular (with ones on the diagonal) then we could write EMBED Equation.3 Let EMBED Equation.3 then EMBED Equation.3 . In matrix form the last equation looks like EMBED Equation.3 can be readily solved for EMBED Equation.3 by forward substitution. That is Set EMBED Equation.3 For i=2 to n with a step of 1 Set EMBED Equation.3 . Then since EMBED Equation.3 or EMBED Equation.3 can be readily solved by backward substitution Set EMBED Equation.3 For i=n-1 to 1 with a step of -1 Set EMBED Equation.3 . To complete the algorithm we need to decompose EMBED Equation.3 by finding L and U from EMBED Equation.3 , or EMBED Equation.3 Then in order from row 1 to n EMBED Equation.3 successively will find all of the EMBED Equation.3 and EMBED Equation.3 . The algorithm is: 1) EMBED Equation.3 for EMBED Equation.3 2) EMBED Equation.3 for EMBED Equation.3 3) For EMBED Equation.3 to EMBED Equation.3 step 1 For EMBED Equation.3 to EMBED Equation.3 step 1 EMBED Equation.3 For j=i+1 to EMBED Equation.3 step 1 EMBED Equation.3 The method also economizes storage since the elements of L and U can be stored back in the original matrix A EMBED Equation.3 Note that once the LU decomposition is completed it can be saved and used t solve for additional right hand sides Fortran code to perform a solution by LU decomposition follows. program solve C C Solve a*x=b C implicit none integer NDIM ! Size of arrays parameter (NDIM=20) C C Declare variables C real*8 a(NDIM,NDIM) ! Array of coefficients real*8 lu(NDIM,NDIM) ! LU decomposition of a real*8 x(NDIM) ! Unknowns real*8 b(NDIM) ! Right hand side integer n ! Number of unknowns integer i,j ! Counters C C Read the array & RHS C read (*,*) n write(6,'(2(A,I3))') 'n = ',n,' NDIM = ',NDIM write(6,'(//)') do i=1,n do j=1,n read (*,*) a(i,j) write (6,'(A,I3,A,I3,A,1PE15.8)') $ ' a(',i,',',j,') = ', a(i,j) end do end do write(6,'(//)') do i=1,n read (*,*) b(i) write(6,'(A,I3,A,1PE15.8)') 'b(',i,') = ',b(i) end do C C Solve the system C call ludecomp(a,lu,n,NDIM) C C Write the LU deconposition C L is from the diagonal down C U is above the diagonal and has ones on the diagonal C write(6,'(//)') do i=1,n do j=1,n write (6,'(A,I3,A,I3,A,1PE15.8)') $ 'lu(',i,',',j,') = ',lu(i,j) end do end do call solvelu(lu,b,x,n,NDIM) C C Write the solution C write(6,'(//)') write(6,'(A)') 'By back substitution:' do i=1,n write(6,'(A,I3,A,1PE15.8)') 'x(',i,') = ',x(i) end do C C Directly find the solution C call matrix_solve(a,b,x,n,NDIM) write(6,'(//)') write(6,'(A)') 'By direct solution:' do i=1,n write(6,'(A,I3,A,1PE15.8)') 'x(',i,') = ',x(i) end do C C Normal stop C stop end subroutine matrix_solve(a,b,x,n,NDIM) C C Solve a set of simultaneous equations by LU decomposition wo. pivoting C implicit none C C Declare variables C real*8 a(NDIM,NDIM) ! Coefficient matrix real*8 b(NDIM) ! Right hand side vector real*8 x(NDIM) ! Unknowns integer n ! Number of equations integer NDIM ! Declared array size real*8 lu(NDIM,NDIM) ! LU decomposition of a C C Decompose & back substitute C call ludecomp(a,lu,n,NDIM) call solvelu(lu,b,x,n,NDIM) C C Normal return C return end subroutine solvelu(lu,b,x,n,NDIM) C C Solve lower*upper*x=b for x C implicit none C C Declare variables C integer NDIM ! Declared size of arrays integer n ! Number of unknowns real*8 lu(NDIM,NDIM) ! LU decomposed matrix real*8 b(NDIM) ! Right hand side real*8 x(NDIM) ! Unknows real*8 sum ! Running sum integer i,k ! Counters C C Run down U to solve intermediate values C x(1)=b(1)/lu(1,1) do i=2,n sum=0.D0 do k=1,i-1 sum=sum+lu(i,k)*x(k) end do x(i)=(b(i)-sum)/lu(i,i) end do C C Run up L to get the final solution C do i=n-1,1,-1 sum=0.D0 do k=i+1,n sum=sum+lu(i,k)*x(k) end do x(i)=x(i)-sum end do C C Normal return C return end subroutine ludecomp(a,lu,n,NDIM) C C Decompose a into lower*upper WITHOUT pivoting C implicit none real*8 lu(NDIM,NDIM) ! LU decomposition of matrix a real*8 a(NDIM,NDIM) ! Input matrix a real*8 sum ! Running sum integer NDIM ! Size of arrays integer n ! Number of unknowns integer i,j,k ! Counters C C Find First column C do i=1,n lu(i,1)=a(i,1) end do C C Find remaining part of column 1 C do j=2,n lu(1,j)=a(1,j)/lu(1,1) end do C C Loop over remaining columns C do j=2,n do i=j,n sum=0.D0 do k=1,j-1 sum=sum+lu(i,k)*lu(k,j) end do lu(i,j)=a(i,j)-sum end do do i=j+1,n sum=0.D0 do k=1,j-1 sum=sum+lu(j,k)*lu(k,i) end do lu(j,i)=(a(j,i)-sum)/lu(j,j) end do end do C C Normal return C return end Determinant Evaluation Some rules that follow from the definition for the evaluation of a determinant are: 1) Adding a multiple of one row to another does not change the value of a determinant. 2) Inverting two rows changes the sign 3) Multiplying a row by a constant will multiply the determinant by the same constant. 4) The determinant of the transpose is equal to the determinant of the original matrix. Example: EMBED Equation.3 EMBED Equation.3 , EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 Inverse Matrix Evaluation Since we are looking for the inverse, first form an augmented matrix with the identity matrix on the right. Consider the example below EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 EMBED Equation.3 Note that the inverse of a symmetric matrix is symmetric. Checking: EMBED Equation.3 . Checks.