【计算方法笔记】矩阵分解求解线性方程组

一般有Doolittle分解和Crout分解。

• Doolittle分解

  将系数矩阵A分解为单位下三角阵(L)和上三角矩阵(U)

• Crout分解

  将系数矩阵A分解为单位上三角阵和下三角矩阵

原理:根据矩阵乘法推导

Doolittle分解，分解的矩阵保存到原系数矩阵，减少内存

clear;A=[5,7,9,10;6,8,10,9;7,10,8,9;5,7,6,5];[n,cl]=size(A);b=[26;18;22;9];% disp(b(4));列向量行向量都可以用数组方式（）访问for i=1:n    for j=1:n        if(i<=j)            A(i,j)=A(i,j)-A(i,1:i-1)*A(1:i-1,j);            %L单位下三角阵        else             A(i,j)=(A(i,j)-A(i,1:j-1)*A(1:j-1,j))/A(j,j);            %U上三角阵        end    endend%LUx=b%Ly=b先解y,再解ux=y,再解xy=zeros(4,1);for i=1:n    y(i,1)=b(i,1)-A(i,1:i-1)*y(1:i-1,1);end%利用matlab矩阵可以用冒号迭代，减少代码数量x=zeros(4,1);for i=n:-1:1    x(i,1)=(y(i,1)-A(i,i+1:n)*x(i+1:n,1))/A(i,i);end

针对特殊系数矩阵（对角占优的三对角矩阵）,Crout分解，此为追赶法。

方便计算记住：
L对角元(i,i)=A对角元(i,i)-A对角元左边(i,i-1)*A对角元上面(i-1,i)
U次对角元(i,i+1)=A(i,i+1)/A对角元(i,i)

A=[2,1,0,0;1,4,1,0;0,1,4,1;0,0,1,2];b=[-3;6;14;2];[n,cl]=size(A);x=zeros(4,1); %将A分解成二对角元的下三角矩阵L和二对角元的上三角矩阵U,若U对角元为1则为Crout分解 A(1,2)=A(1,2)/A(1,1); for i=2:n     A(i,i)=A(i,i)-A(i,i-1)*A(i-1,i);     if(i<=n-1)         A(i,i+1)=A(i,i+1)/A(i,i);     end end %解LUx=b,先解Ly=b,再解Ux=b; x(1)=b(1)/A(1,1); for i=2:n     x(i)=(b(i)-A(i,i-1)*x(i-1))/A(i,i); end x(n)=x(n); for i=n-1:-1:1     x(i)=x(i)-A(i,i+1)*x(i+1); end format rat;%以分数展示结果 %format解除分数显示 disp(A); disp(x);