怎样用LU分解法解线性方程组

怎样用LU分解法解线性方程组

Ax=B，改写成Ly=B，Ux=y的方程组。就相当于将A=LU分解成了两个矩阵。称为矩阵A的三角分解，或LU分解。如果L为单位下三角阵，则叫Doolittle分解，若U为单位上三角阵，则叫Crout分解。只要A的各顺序主子式不为零，则A可唯一分解成一个单位下三角阵L与一个上三角阵U的乘积。
•设Ax=b，A=LU，则Ax=LUx=b
于是令Ux=y，则Ly=b
这样原来方程能化为两个简单方程组

%A*x=b
A=[1 1;1 -1]
b=[2;3]
[L,U] = lu(A)
x = U\(L\b)

Ax=B，改写成Ly=B，Ux=y的方程组。就相当于将A=LU分解成了两个矩阵。称为矩阵A的三角分解，或LU分解。如果L为单位下三角阵，则叫Doolittle分解，若U为单位上三角阵，则叫Crout分解。只要A的各顺序主子式不为零，则A可唯一分解成一个单位下三角阵L与一个上三角阵U的乘积。
•设Ax=b，A=LU，则Ax=LUx=b
于是令Ux=y，则Ly=b
这样原来方程能化为两个简单方程组

下面是LU分解的Fortran子程序 希望可以有所帮助

!求解au=b,u
!n表示为方程维数
subroutine lu(a,b,n,u)
implicit real(8) (a-h,o-z)
real(8)::a(n,n),b(n),u(n),a_bak(n,n),b1(n),aL(n,n),aU(n,n),y(n)
!exchange rows
do i=1,n
    tmpMax=0.d0
do ic=i,n
if(tmpMax    tmpMax=dabs(a(ic,i))
   i_rec=ic
 endif
enddo
    if(i_rec.ne.i) then
do jc=i,n
tmp=a(i,jc)
a(i,jc)=a(i_rec,jc)
a(i_rec,jc)=tmp
enddo
tmp=b(i)
b(i)=b(i_rec)
b(i_rec)=tmp
        endif 
!decomposition
 do j=i,n
    tmp=0.d0
     do k=1,i-1
tmp=tmp+aL(i,k)*aU(k,j)
     enddo
    aU(i,j)=a(i,j)-tmp
    tmp=0.d0
     do k=1,i-1
tmp=tmp+aL(j,k)*aU(k,i)
     enddo
     aL(j,i)=(a(j,i)-tmp)/aU(i,i)
 enddo
enddo
!find answer
do i=1,n
    tmp=0.d0
     do j=1,i-1
     tmp=tmp+aL(i,j)*y(j)
     enddo
    y(i)=b(i)-tmp
enddo
do i=n,1,-1
    tmp=0.d0
   do j=i+1,n
    tmp=tmp+aU(i,j)*u(j)
   enddo
    u(i)=(y(i)-tmp)/aU(i,i)
enddo
end