工程计算：LU分解（Doolittle及Cholesky）

LU分解

目的：解决采用高斯消元法大量重复计算导致的效率低下的问题。
高斯消元过程其实就是一个矩阵的三角分解过程：
由系数矩阵$A$的求解过程：
$$A^{(k+1)}=L_k{A}^{(k)}$$
我们可以总结出：
$$A=LU$$
即LU分解（Doolittle分解）
其中，$L为：单位下三角矩阵$ $U为：上三角矩阵$。
定理：若A的所有顺序主子式不为零，则A存在唯一的LU分解。
LU分解流程：

举个“栗子”–Doolittle

代码：

A = [3 2 -1 2
     1 4 0 2
     2 1 2 -1
     1 1 -1 3];

b=[-2
    1
    3
    4];

%初始化
n= size(A,1);
L=zeros(n);%单位矩阵
U=eye(n);%零矩阵
%计算
%对应从第U的第1行
for r=1:n
   for j=r:n
       U(r,j)=A(r,j);
       for k=1:r-1
           U(r,j)=U(r,j)-L(r,k)*U(k,j);
       end
   end
   
   for i=r:n
       L(i,r) =A(i,r);
       for k = 1:r-1
           L(i,r)=L(i,r)-L(i,k)*U(k,r);
       end
       L(i,r)=L(i,r)/U(r,r);
   end
end
L
U

%LY = b
y = zeros(n,1);
for i= 1:n
    y(i)=b(i);
    for j = 1:i-1
        y(i) = y(i)-L(i,j)*y(j);
    end
end

%UX= y
x = zeros(n,1);
x(n)=y(n)/U(n,n);
for i = n:-1:1
    x(i)=y(i);
    for j=i+1:n
        x(i)=x(i)-U(i,j)*x(j);
    end
    x(i)=x(i)/U(i,i);
end
x

我们也可以使用sum（）函数来代替

       L(i,r) =A(i,r);
       for k = 1:r-1
           L(i,r)=L(i,r)-L(i,k)*U(k,r);
       end

这样可以减少循环嵌套，降低时间复杂度。

A = [3 2 -1 2
     1 4 0 2
     2 1 2 -1
     1 1 -1 3];

b=[-2
    1
    3
    4];
[n,n]=size(A);
L=zeros(n);
U=zeros(n);
x=zeros(n,1);
y=zeros(n,1);
for r=1:n
    for i= r:n
        U(r,i)=A(r,i)-sum(L(r,1:r-1).*U(1:r-1,i)');
        L(i,r)=(A(i,r)-sum(L(i,1:r-1).*U(1:r-1,r)'))/U(r,r);
    end
end
L,U
for i=1:n
    y(i)=b(i)-sum(L(i,1:i-1).*y(1:i-1)');
end
for j=n:-1:1
     x(j)=(y(j)-sum(U(j,j+1:n).*x(j+1:n)'))/U(j,j);
end
x

运行结果：

L =
1.0000 0 0 0
0.3333 1.0000 0 0
0.6667 -0.1000 1.0000 0
0.3333 0.1000 -0.2593 1.0000

U =
3.0000 2.0000 -1.0000 2.0000
0 3.3333 0.3333 1.3333
0 0 2.7000 -2.2000
0 0 0 1.6296

y =
-2.0000
1.6667
4.5000
5.6667

x =
-0.5909
-1.3409
4.5000
3.4773

举个“栗子”–Cholesky

对称矩阵的分解法

代码：

A=[5 -4 1 0
   -4 6 -4 1
   1 -4 6 -4
   0 1 -4 5];

b=[2
   -1
   -1
   2];

n = size(A,2);
d = zeros(n,1);
L = eye(n);

for k=1:n
    d(k) = A(k,k);
    for j=1:k-1
        d(k) = d(k) - L(k,j)*d(j)*L(k,j);
    end
    for i = k+1:n
        L(i,k) = A(i,k);
        for j=1:k-1
            L(i,k) = L(i,k) - L(i,j)*d(j)*L(k,j);
        end
        L(i,k)=L(i,k)/d(k);
    end
end
L
D = diag(d);
L*D*L'
A

%Lz =b
z=zeros(n,1);
for i=1:n
   z(i)=b(i);
   for j=1:i-1
       z(i) = z(i) - L(i,j)*z(j);
   end
end
 
%Dy = z
y=zeros(n,1);
for i=1:n
    y(i)=z(i)/d(i);
end

%L'x = y
x=zeros(n,1);
for i=n:-1:1
    x(i)=y(i);
    for j=i+1:n
        x(i)=x(i)- L(j,i)*x(j);
    end
end
x

运行结果：

L =
1.0000 0 0 0
-0.8000 1.0000 0 0
0.2000 -1.1429 1.0000 0
0 0.3571 -1.3333 1.0000

ans =
5 -4 1 0
-4 6 -4 1
1 -4 6 -4
0 1 -4 5

A =
5 -4 1 0
-4 6 -4 1
1 -4 6 -4
0 1 -4 5

x =
1.0000
1.0000
1.0000
1.0000