matlab实现用不同方法求线性方程组AX=b的解 高斯消元法、列主消元法、杜立特尔分解法

代码如下

A=[2 1 1;1 3 2;1 2 2];
b=[4 5 6]';
fprintf("高斯消元法:\n")
x = GaussPivot(A,b)%%高斯消元法
fprintf("列主消元法:\n")
x = col_e(A,b)%%列主消元法
fprintf("杜利特尔消元法:\n")
x = LUfenjiefa(A,b)%%杜利特尔消元法

分别建立函数文件

注意文件名匹配

function x = GaussPivot(A,b)
%高斯
[m,n]=size(A);
if m~=n,error("系数矩阵须为方阵");end
nb=n+1;
Aug=[A b];
for k=1:n-1
    %选主元(max函数可以传回最大值和最大值索引)
    [big,i] = max(abs(Aug(k:n,k)));
    ipr = i+k-1;
    if ipr~=k
        Aug([k,ipr],:) = Aug([ipr,k],:);%交换：将绝对值大的数作主元
    end
    for i=k+1:n
        factor = Aug(i,k)/Aug(k,k);
        Aug(i,k:nb) = Aug(i,k:nb)-factor*Aug(k,k:nb);
    end
end
x = zeros(n,1);%创建解向量
x(n)=Aug(n,nb)/Aug(n,n);%先求出最后一个值
for i=n-1:-1:1
    x(i)=(Aug(i,nb)-Aug(i,i+1:n)*x(i+1:n))/Aug(i,i);
end
end

function x=LUfenjiefa(A,b)
n=length(b);
k=2;
X=A
Y=b
U(1,1:n)=A(1,1:n);
L(2:n,1)=A(2:n,1)/U(1,1);
for k=2:n
    U(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
    L(k+1:n,k)=(A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);
end
L
U
%用向前消去法解下三角方程组Ly=b
y=zeros(n,1);
y(1)=b(1);
for k=2:n
    y(k)=b(k)-L(k,1:k-1)*y(1:k-1);
end
y
%用回代法解上上角方程组Ux=y
x=zeros(n,1);
x(n)=y(n)/U(n,n);
for k=n-1:-1:1
    x(k)=(y(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
end
end

function [x] =col_e(A,b)
%col_e 列主元消去法(Out principal component elimination)计算线性方程组
%   INPUT:
%       A:方程系数矩阵a11,a12,a21...ann
%       b:方程的增广矩阵b1,b2...bn
%   OUTPUT:
%       x:A/b消元
[m,n]=size(A);  %A的行数和列数（A的大小）
x=zeros(n,1);
if m~=n
    error('A不是矩阵');
    return;
end
if m~=size(b)
    error('b的列数与A的行列数不同');
    return;
end
if rank(A)~=rank([A,b])
    error('A的矩阵的秩和增广矩阵的秩不同，方程不存在唯一的解');
    return;
end

for i=1:n-1
    [max_a,m]=max(abs(A(:,i)));%找到每列最大项，max_a为该项值，m为行数
    if max_a==0
        error('奇异方程');
        return;
    end
    
    if(m~=i)
        %调换矩阵行
        A([i,m],:)=A([m,i],:);
        b([i,m],:)=b([m,i],:);
    end
    for k=i+1:n
        for j=i+1:n
            A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);
        end
        b(k)=b(k)-b(i)* A(k,i)/A(i,i);
        A(k,i)=0;
    end
end
x(n)=b(n)/A(n,n);
for i=n-1:-1:1
    sum=0;
    for j=i+1:n
        sum=sum+A(i,j)*x(j);
    end
    x(i)=(b(i)-sum)/A(i,i);
end
disp(A);
end

建立好后直接运行就行