数值分析3-解线性方程组的高斯消去法、LU分解法及列主元消去法的matlab程序和调试方法

大家好，欢迎来到IT知识分享网。

对于形如Ax=b的线性方程组，在线性代数中是通过求逆的方式求解的，即x=A^-1b,而在数值分析中，解线性方程组的方法是通过直接法或者迭代法来实现的，今天写的三个程序都属于直接法，分别为高斯消去法、LU分解法和列主元消去法。
还有全主元消去，平方根法，追赶法，其中平方根法，追赶法的思路都高斯消去法差不多，稍微修改一下程序就行，而全主元消去法相较于列主元消去法而言，判断大小和换行换列会更麻烦一些，用循环语句能用到吐，没啥需求的话我实在是懒得写；
高斯消去法，就是线性代数中通过把系数矩阵化为行阶梯矩阵然后求解的方法，而LU分解法和高斯消去法十分相似，只不过是回代过程有所不同，所以我把这两个程序放在一起写了。程序中用到了较多的循环语句，要仔细分析一下，不然写的过程中很容易出错。
首先是高斯消去法：
代码块1：高斯消去法

clear;
clc;
fprintf('高斯消去法解线性方程组:\n')
n=3;
A=[1 2 3;2 5 2;3 1 5];
b=[14 18 20];                 %注意，这里我用的S是行向量，计算时注意转置
%n=input('请输入系数矩阵A的阶数:n=');
%A=input('请输入系数矩阵A:A=');
%b=input('请输入结果向量b:b=');%b为行向量，计算时需要转置
%消元计算---->产生行阶梯矩阵：
A1=A;                         %为了保留原始参数，以下用A1存储产生的上三角矩阵
b1=b;
for j=1:n                     %从第一列到第n列
    for i=j+1:n               %每一行从主元素的下侧第一个元素开始
        if A1(j,j)~=0         %每次的主元素都不可为零，这是能使用高斯消去法的前提
            M(i,j)=A1(i,j)/A1(j,j);
            for k=j:n         %第i行j列开始计算到第i行n列
                A1(i,k)=A1(i,k)-M(i,j)*A1(j,k);
            end
                b1(i)=b1(i)-M(i,j)*b1(j);   
        else 
            error('求解出错，高斯消去法不适用于此方程组')
        end
    end
end
%回代求解：
 x(n)=b1(n)/A1(n,n);    %先求出解向量的最后一个元素，用于回代
 i=n-1;                 %从n-1行开始，回代求解
while i>0               %从第n-1行到第一行，依次回代
     s=0;               %对s赋初值0，
     for j=i+1:n
         s=A1(i,j)*x(j)+s;%计算：（已知解向量×系数矩阵对应元素）的和
     end
     x(i)=(b1(i)-s)/A1(i,i);%回代求解，解向量中已知元素+1.用于回代求x(i-1)
     i=i-1;
end
fprintf('经过消元回代，此方程组的解为:\n')
disp(x);

%可以用inv(A)*b'直接求解验证计算
%inv(A)*b'

对于有的线性方程组是无法使用高斯消去法求解的，能使用高斯消去法的前提是系数矩阵A的顺序主子式不为0，这也是系数矩阵能LU分解的充分不必要条件，还有，如果系数矩阵A在计算过程中的主元素绝对值过小的的话，得到的结果误差会很大，这也是数值计算中要避免的一种情况，即：绝对值较小的数做除数；这种线性方程组可以考虑使用列主元消去法或者全主元消去法来计算。（但是matlab计算精度很高，就硬算也没啥问题）

所谓LU分解法，其实就是高斯消去法的一种变形，对于顺序主子式不为零的n阶矩阵A，都可以分解成A=L*U的形式，其中L为对角元素为1的下三角矩阵，U为上三角矩阵；因此Ax=b可以等价为：Ux=y，Ly=b，分成两步来回代求解，其代码如下所示：
代码块2：LU分解法解线性方程组

clear;
clc;
fprintf('LU分解法解线性方程组:\n')
n=3;
A=[1 2 3;2 5 2;3 1 5];
b=[14 18 20];                 %注意，这里我用的是行向量，计算时注意转置
%n=input('请输入系数矩阵A的阶数:n=');
%A=input('请输入系数矩阵A:A=');
%b=input('请输入结果向量b:b=');%b为行向量，计算时需要转置
k=1;
%消元计算---->产生上三角矩阵：
A1=A;                         %为了保留原始参数，以下用A1存储产生的上三角矩阵
%LU分解：
for j=1:n                     %从第一列到第n列
    for i=j+1:n               %每一行从主元素的下侧第一个元素开始
        if A1(j,j)~=0         %每次的主元素都不可为零，这是能使用高斯消去法的前提
            M(i,j)=A1(i,j)/A1(j,j);
            for k=j:n         %第i行j列开始计算到第i行n列
                A1(i,k)=A1(i,k)-M(i,j)*A1(j,k);
            end  
        else 
            error('求解出错，系数矩阵无法LU分解')
        end
    end
end
for i=1:n
    M(i,i)=1;%M即为所求的下三角矩阵L，上三角矩阵U=A1
end   
L=M;
U=A1;
%回代求解：
y(1)=b(1);
for i=2:n
    s=0;
    for k=1:i-1
        s=s+L(i,k)*y(k);
    end 
    y(i)=b(i)-s;    
end
x(n)=y(n)/U(n,n);
i=n-1;
while i>0
    s=0;
    for k=i+1:n
        s=s+U(i,k)*x(k);
    end
    x(i)=(y(i)-s)/U(i,i);
    i=i-1;
end
fprintf('此方程组的解为:\n')
disp(x);
%用inv(A)*b'直接求解验证计算

列主元消去法是高斯消去法的一种改进方法，它可以避免高斯消去法中主元素过小导致的误差；有利必有弊，列主元消去法计算量要大一些，因为他要有一个比较大小和重新排序换行的问题，这里我写的比较麻烦，用了好多个循环，我自己看的都要晕了，代码块如下：
代码块3：列主元消去法解线性方程组：

%列主元消去法解线性方程组
clear;
clc;
fprintf('列主元消去法解线性方程组:\n')
n=3;
%A=[1 2 3;2 5 2;3 1 5];
%b=[14 18 20];                 %注意，这里我用的S是行向量，计算时注意转置
n=input('请输入系数矩阵A的阶数:n=');
A=input('请输入系数矩阵A:A=');
b=input('请输入结果向量b:b=');%b为行向量，计算时需要转置
%消元计算---->产生上三角矩阵：
A1=A;                         %为了保留原始参数，以下用A1存储产生的上三角矩阵
b1=b;
for j=1:n                     %从第一列到第n列
    for i=j+1:n               %每一行从主元素的下侧第一个元素开始
        for k=i:n%这个循环选出第j列最大的元素做主元
            if abs(A1(k,j))>abs(A1(j,j))%if语句，别忘了绝对值，ps：冒泡排序，真滴好用
                for k1=j:n              %用来交换第k行和第j行
                    a1=A1(j,k1);
                    a2=A1(k,k1);
                    A1(j,k1)=a2;
                    A1(k,k1)=a1;
                end
                t=b1(k);
                b1(k)=b1(j);
                b1(j)=t;%别忘了结果向量，也要相应的交换行
            end
        end
%上述增补的部分就是列主元消去的精髓，下面的部分和高斯消去一样
        if A1(j,j)~=0         %每次的主元素都不可为零，这是能使用高斯消去法的前提
            M(i,j)=A1(i,j)/A1(j,j);
            for k=j:n         %第i行j列开始计算到第i行n列
                A1(i,k)=A1(i,k)-M(i,j)*A1(j,k);
            end
                b1(i)=b1(i)-M(i,j)*b1(j);   
        else 
            error('求解出错，高斯消去法不适用于此方程组')
        end
    end
end
%回代求解：
 x(n)=b1(n)/A1(n,n);    %先求出解向量的最后一个元素，用于回代
 i=n-1;                 %从n-1行开始，回代求解
while i>0               %从第n-1行到第一行，依次回代
     s=0;               %对s赋初值0，
     for j=i+1:n
         s=A1(i,j)*x(j)+s;%计算：（已知解向量×系数矩阵对应元素）的和
     end
     x(i)=(b1(i)-s)/A1(i,i);%回代求解，解向量中已知元素+1.用于回代求x(i-1)
     i=i-1;
end
fprintf('经过消元回代，此方程组的解为:\n')
disp(x);

%其实这样子写会计算时多判断好多次，不够简洁，

列主元消去比高斯消去增加了16-29行，就是用来判断大小和排序的环节，看着这一排的for和if，真滴让人头大呀！
代码中用s来实现累加器的功能，用for来实现变量递增的循环，用while实现变量递减的循环；
我是直接定义了矩阵A，结果向量b，也可以把我直接给出的矩阵和结果向量注释掉，把下方我写的自己输入A,b的语句放出来。
写程序的过程中谁也不可能说永远能一步到位，一下子就写出正确的来，所以出现错误后，可以按下图方式打点，运行后按步进方式，一步一步的运行并观测数据，排查错误出现在那一步；


应该没什么要注意的了，因为我测试的数据比较少，不知道哪里还有没发现的bug，欢迎大家交流呀~~~