如何用matlab三角矩阵分解？

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

matlab有多种LU分解程序下面算一种：

function [L,U]=myLU(A)

%实现对矩阵A的LU分解，L为下三角矩阵

A

[n,n]=size(A)

L=zeros(n,n)

U=zeros(n,n)

for i=1:n

L(i,i)=1

end

for k=1:n

for j=k:n

U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)')

end

for i=k+1:n

L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k)

end

end

用法，在控制台输入

A=[1 2 3 -4-3 -4 -12 132 10 0 -34 14 9 -13]

---