如何用matlab三角矩阵分解？

Question

如题所述

Answer 1

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<dabs(a(ic,i))) then

   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