LU ¤À¸Ñªk

 

¦pªG§Ú­Ì¥i¥H§â½u©Ê¥N¼Æ°ÝÃDªº Ax = b ¤¤ªº¯x°} A ¤À¸Ñ¦¨ A = L U ¡A¨ä¤¤ L ¬O¤U¤T¨¤¯x°} (lower triangular matrix)¡BU ¬O ¤W¤T¨¤¯x°} (upper triangular matrix)¡A«h­ì¥»¨D¸Ñ Ax = b ªº°ÝÃD´N¥i¥H¤À¦¨¨â¶¥¬q¨Ó¨D¸Ñ¡A¤À§O¬O Ly = b ¤Î Ux = y¡A¥¦­Ì³£¥i¥H¾A¥Î«e¤@­Ó³æ¤¸¤¶²Ð¹Lªº back substitution ¤èªk¡A¥ý±o¥X y ¦A±o¥X x¡A¦p¦¹§â¸Ñ¨D±o¡C­È±oª`·N¥B¦³½ìªº¬O¡A¹ï©ó¥ô¦óµ¹©wªº¤è¯x°} A¡A¬O§_¤@«ß³£¯à³Q¤À¸Ñ¦¨ L ©M U¡H

 

«Øºc©ÊªºÃÒ©ú

­nÃÒ©ú¨ä¦s¦b¡A§Ú´N®Ú¾Ú©Òµ¹ªº¸ê®Æ§â¸Óª«¥ó°µ¥X¨Ó¡A¥Ñ©ó³o¬O¤@¯ë©Ê¦ÓµL¯S®í©Êªº¡A¦]¦¹ÃÒ©ú¤F¸Óª«¥óªº¦s¦b©Î¸Ó¯S©Êªº¦¨¥ß¡C

¦b¥»³æ¤¸¤¤¡A§Ú­Ì­nÃÒ©ú¥ô¦ó¥ô¤è°} A ³£¥i¥H³Q¤À¸Ñ¼g¦¨ A = L U ¨â­Ó¤è¯x°}ªº§Î¦¡¡A¦b§@ªk¤W¡A§Ú­Ì´N¥ý¤jÁx¦a¼g¤U¦p½Ò¥» (2.3.2) ¦¡ªº¦¡¤l¡G

LU_alpha_beta

¦b¤£¥¢¤@¯ë©Êªº±¡ªp¤U¡A§Ú­Ì¥i¥H¥O ©Ò¦³¹ï¨¤½u¤Wªº aii ³£µ¥©ó 1¡C³o¼Ë°µ¯uªº¨S°ÝÃD¶Ü¡H¯uªº¥i¥u¦³¤@²Õ L »P U ¥i¥H±o±o¨ì¶Ü¡H§Ú­Ì¥ý¤£­n©È³Â·Ð¡A¸ÕµÛ¤@­Ó¤@­Ó¨Ó¬Ý aij »P bjk ­¼¦b¤@°_¦p¦ó§Î¦¨ aik ¡C³Ì²³æªº¬O²Ä¤@­Ó¡A¤w¸g»¡¦n¤F aii ³£µ¥©ó 1¡A¦]¦¹ b11 = a11¡C¦A©¹¤U¬Ý¤@­Ó¡A¬O a21b11 = a21¡A¦Ó b11 ­è¤~¤v¸g±oª¾¬O a11¡A¦]¦¹ a21 ´N¦Û°Êª¾¹D¤F¡C¦P²z¡A a31 »P a41 ¤]³£¬OÃþ¦üªº¤èªk»´©ö±oª¾¡C¦A¨Ó¡A§Ú­Ì¬Ý U ¤è°}²Ä¤G¦C­¼¦V L ¤è°}¨º¥|¦Cªºµ²ªG¡A¡]°O±o a ªº 11¡B21¡B31¡B41 ³£¤v¸gª¾¹D¤F¡^­¼ L ªº²Ä¤@¦C·|¹G¥X b12 ªº­È¡C ­¼¨ì L ªº²Ä¤G¦C®É¡A¥Ñ©ó a22 ¬° 1¡A´N½T©w¤F b22 ¡C¦b b12¡Bb22 ³£¬O¤vª¾ªº±¡ªp¤U¡A­¼²Ä¤T¦C·|©w¥X a32¡A¦P²z­¼²Ä¥|¦C·|©w¥X a42 ¡A¦p¦¹ L ¤è°}ªº²Ä¤G¦æ¤]©w§¹¤F¡C¦A±µ¤U¨Ó¬O®³ U Á«°}ªº²Ä¤T¦æ­¼¨ì L ¤è°}ªº¥|­Ó¦æ¡A¨Ì§Ç©w¥X b13¡Bb23¡Bb33¡Ba43¡C³Ì«á¡A¬Û¦Pªºµ¦²¤·|¨Ì§Çµ¹¥X b14¡Bb24¡Bb34¡Bb44¡A¦p¦¹©Ò¦³ L ¤Î U ¤è°}ªº¤¸¯À´N³£°ß¤@½T©w¤F¡C§Ú­Ì¦]¦¹«Øºc¦¡¦aÃÒ©ú¤F¡]¦b 4x4 ªº¨Ò¤l¤¤¡^L U = A ªº¥i¦æ©Ê¡C

Á`¨¥¤§¡A­Y¥H L »P U ¼g¦b¤@°_ªººò±K§Î¦¡¨Ó¬Ý

LU_combined

¤W¦¡¤¸¯À¨Ì§Ç©w¥Xªº¶¶§Ç¦p¤U¡G

1 5 9 13
2 6 10 14
3 7 11 15
4 8 12 16

 

 

°Æµ{¦¡ªº¨Ï¥Î

¥»¸`©Ò¥Î¨ìªº°Æµ{¦¡¦³¨â­Ó¡A¤À§O¬O§@ LU ¤À¸Ñªº ludcmp¡A¥H¤Î°µ°t¦X¨ä§@¤Ï¦V¥N¦^ªº lubksb ¡C ³o¨â­Ó°Æµ{¦¡ªº©I¥s¤Þ¼Æ (argument) ¦w±Æ»P½Ò¤å¤º»¡©ú¤À§O¦p¤U¡G

SUBROUTINE ludcmp(a,n,np,indx,d)
INTEGER n,np,indx(n),NMAX
REAL d,a(np,np),TINY
PARAMETER (NMAX=500,TINY=1.0e-20) Largest expected n, and a small number.
Given a matrix a(1:n,1:n), with physical dimension np by np, this routine replaces it by the LU decomposition of a rowwise permutation of itself. a and n are input. a is output, arranged as in equation (2.3.14) above; indx(1:n) is an output vector that records the row permutation effected by the partial pivoting; d is output as ¡Ó1 depending on whether the number of row interchanges was even or odd, respectively. This routine is used in combination with lubksb to solve linear equations or invert a matrix.
INTEGER i,imax,j,k

ludcmp »Ý­n±q¥Dµ{¦¡°t¦X«Å§iªº¬O a(np,np)¡A¥H¤Î¶Ç¤J¾ã¼Æ­È n¡A´N³o¼Ë¦Ó¤w¡A³Ñ¤Uªº¤Þ¼Æ indx ¤Î d «h¬O¬°¤F¤§«á»Î±µ lubksb ¸Ñ½u©Ê¥N¼Æ°ÝÃD¤§¥Î¡]ludcmp ¥u§@¯x°} A ªº¤µ¤À¸Ñ¡^¡C

SUBROUTINE lubksb(a,n,np,indx,b)
INTEGER n,np,indx(n)
REAL a(np,np),b(n)
Solves the set of n linear equations A ¡P X = B. Here a is input, not as the matrix A but rather as its LU decomposition, determined by the routine ludcmp. indx is input as the permutation vector returned by ludcmp. b(1:n) is input as the right-hand side vector B, and returns with the solution vector X. a, n, np, and indx are not modified by this outine and can be left in place for successive calls with different right-hand sides b. This routine takes into account the possibility that b will begin with many zero elements, so it is efficient for use in matrix inversion.
INTEGER i,ii,j,ll
REAL sum

lubksb ¤@¦¸¥u§@¤@­Óª½¦æªº b ¦V¶q¡A¤£¹³¬O Gauss-Jordan ¨º¼Ë°µ column-augmented ¤§¦h­« b ªº²Õ¦X¡C ³o¬O¨S¦³¥²­nªº¡A¦]¬°­ìÃD¥Øµ¹ªº¯x°} A ·~¸g¤À¸Ñ¬° L »P U¡A¨D¸Ñ§Ö¤S¦³®Ä²v¡C ©I¥s¹L lubksb ¤§«á¡Aµª®× x ·|¦s¦b b ªº¦ì¸m¶Ç¥X¨Ó¡C

¦pªG­n¨D¸Ñ A ªº¤Ï¯x°} A-1¡A§Ú­Ì¥i¥H¦b¥Dµ{¦¡¤¤¦w±Æ¦h­Ó b ¦V¶q¡A¨C¤@­Ó¨Ì§Ç³£¬Oºc¦¨³æ¦ì¯x°}ªº¤@ª½¦æ¡A¦p¦¹¦h¦¸©I¥s lubksb ¸Ñ¥X¦U¹ïÀ³ªº¦æ¦V¶q x ¨Ó¥[¥H²Õ¦X¡A´N·|±o¨ì¤Ï¯x°} A-1¡C¸Ô²Ó§@ªk¡]µ{¦¡½X¡^½Ò¤å¤¤¦³¡C

­Y­n¨D¦æ¦C¦¡­È¡A«h¦b¥Dµ{¦¡¤¤¼g­Ó²³æªº°j°é¡A§â¤À¸Ñ«á¶Ç¥X¤§¯x°} a ¨ä¤W¤§¹ï¨¤½u¤Wªº¤¸¯À¡]§Y bii ³s­¼°_¨Ó¡A´N¬O det(U)¡A³o¬O¦]¬°¤T¨¤¯x°}ªº¦æ¦C¦¡¯S§O®e©ö¡A´N¬O¹ï¨¤¤¸¯Àªº³s­¼¿n¡A³oùئ³¥Î¨ì det(A) = det(LU) = det(L)det(U) = 1.det(U) = det(U)¡C¡^¡C ¸Ô²Ó§@ªkªºµ{¦¡½X¥ç¥i¨£½Ò¤å¡C³oùتº­«ÂI¬O¡AU ¯x°}ªº¹ï¨¤¤¸¯À¿n¡AÁÙ­n¦A­¼¤W§Q¥Î ludcmp ¦^¶Çªº d ­È¡A¤~¬O³Ì«á¥¿½Tªº¦æ¦C¦¡­È¡]¥Ñ©ó pivoting ¤§¦æ¦C¥æ´«©Ò³y¦¨ªº¥¿­t®t²§¬O¦s¦b d ùØ­±¡^¡C

 

½Ò¤å¸É¥Rª¾ÃÑ

¦pªG§A»Ý­n³B²z¥H¤Uªº°ÝÃD¡Gµ¹©w A ¤Î B¡A¨D A-1B¡A«h½Ò¥»§@ªÌ´£¿ô§AÀ³¸Ó§â³o­Ó°ÝÃD·í§@ AX = B ¦Ó¨D¸Ñ X ¨Ó³B²z¡]¦]¬° X ¥¿¬O A-1B¡^¡A¦Ó¤£¬O¥ý¥Î LU ¤À¸Ñªk¨D¤Ï¯x°} A-1¡AµM«á¦A¥h°µ¾ã®M¯x°}­¼ªk¡A¦p¦¹¤ñ°_ª½±µ¸Ñ X ¶O®É¥B¤£¦p¨äºë·Ç¡C

 

¥Dµ{¦¡½d¨Ò

program lu_main
implicit none
integer n, np, i, j
parameter (np=10)
integer indx(np)
real a(np,np), b(np), d
character*40 filename1, filename2

write (*,*)'Please tell me the dimension n of the matrix A(n,n):'
read (*,*) n

write (*,*) 'Please give the file name that contains A :'
read (*,*) filename1
write (*,*) 'Please give the file name that contains b :'
read (*,*) filename2

open (unit=20, file=filename1)
do i=1,n
  read (20,*) ( a(i,j), j=1,n )
enddo
close(20)

open (unit=20, file=filename2)
do i=1,n
  read (20,*) b(i)
enddo
close(20)

call ludcmp(a,n,np,indx,d)

call lubksb(a,n,np,indx,b)

write (*,*) 'The solution matrix is :'
do i=1,n
  write (*,*) b(i)
enddo

do i=1,n
  d = d * a(i,i)
end do
write(*,*) 'The determinant is', d

end