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

最新推荐文章于 2021-05-24 00:31:31 发布

BFNJ

最新推荐文章于 2021-05-24 00:31:31 发布

阅读量4.5k

点赞数 1

分类专栏：计算方法笔记文章标签： matlab

本文链接：https://blog.csdn.net/zhulinzhulinlin/article/details/78187251

版权

计算方法笔记专栏收录该内容

5 篇文章 1 订阅

订阅专栏

矩阵分解

定理：

设A如果为非奇异矩阵，则必存在排列矩阵p以及单位下三角矩阵L和非奇异上三角矩阵U,使得

$P A = L U$ $PA=LU$

原理:根据矩阵乘法推导

Doolittle分解( Matlab程序)

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
    end
end

%LUx=b
%Ly=b先解y,再解ux=y,再解x
y=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

注意：

上面为自己根据理解编写的程序，还需改进，如：可以用二维数组存放增广矩阵A(1:n;1:n+1)，在矩阵A基础上进行处理，这样系数矩阵的Doolittle分解与解方程组 $Ly=b$ 可以同时进行，最后增广矩阵最后一列就是 $y$ 的值，再解 $UX=y$ 即完成了计算。

Crout分解( Matlab程序)

针对特殊系数矩阵（对角占优的三对角矩阵）,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);