Matlab 中部分旋转的高斯消元法

Ilja Kostenko 2023年1月30日
  1. 在 Matlab 中使用高斯消元法
  2. 在 Matlab 中使用旋转作为辅助功能
Matlab 中部分旋转的高斯消元法

本文将帮助读者了解如何在 Matlab 中使用高斯消元法。

在 Matlab 中使用高斯消元法

高斯方法是求解线性代数方程 (SLA) 系统的经典方法。

这是一种顺序排除变量的方法;当使用基本变换时,方程组被简化为等价的三角形系统。

所有系统变量都按顺序定位,从最后一个(按编号)开始。

Gauss-Jordan 算法在第一步将第一行除以系数 A[1,1]

然后该算法将第一行添加到具有这样系数的其余行中,使得它们在第一列中的系数变为零 - 为此,当将第一行添加到第 i 行时,有必要将其乘以 -A[i,1]

对于矩阵 A 的每个操作(除以一个数字并与另一行相加),用向量 b 执行相应的操作;从某种意义上说,它的行为就像是矩阵 A 的第 m+1 列。

结果,在第一步结束时,矩阵 A 的第一列将变为单列(即,它将在第一行包含一个,在其余行包含零)。

同理,执行算法的第二步;现在只考虑第二列和第二行。首先,第二行除以 A[2,2],然后从所有其他行中减去系数,例如重置矩阵 A 的第二列。

依此类推,直到我们处理矩阵 A 的所有行或列。如果 n= m,那么从算法的构造中可以明显看出,矩阵 A 将变成奇异矩阵,这正是我们所需要的。

在 Matlab 中使用旋转作为辅助功能

当然,上述方案是不完整的。它仅在每个第 i 步元素 A[i, i] 不为零时才有效 - 否则,我们将无法通过将第 i 行添加到当前列中的剩余系数来将其归零.

有一个选择参考元素的过程(总而言之,旋转)以使算法在这种情况下工作。它包括这样一个事实,即矩阵的行和列被重新排列,以便所需的元素 A[i, i] 具有非零数。

请注意,在计算机上,行的重排比列的重排要容易得多:毕竟,在交换某些两列时,你需要记住这两个变量交换了位置,以便在还原答案时,正确还原哪个 answer 属于哪个变量。

重新排列行时,无需执行此类附加操作。

换句话说,在使用部分旋转启发式执行 Gauss-Jordan 算法的第 i 阶段之前,需要在索引从 i 到 n 的元素中找到第 i 列的最大模,并交换第 i 行包含此元素的行。

function x = pivot(~, ~)
% Program to Solve Ax = b using Gaussian Elimination Method with partial-
% pivoting (Row-exchange)
%==========================================================================

% OUTPUT:
%==========================================================================

% Solution vector, x

%==========================================================================
A=[1 2 3; 4 5 6];
b=[1 -2 5];
[m,n] = size(A); % checking the size of matrix


%==========================================================================
% Initialization
x = zeros(m,1);
l = zeros(m,m-1);

%==========================================================================
% Main Program
%==========================================================================

% Reducing Matrix A to upper triangular form

for k = 1:m-1
    % =========Performing Partial-pivoting=================================
        for p = k+1:m
            if (abs(A(k,k)) < abs(A(p,k)))
                A([k p],:) = A([p k],:);
                  b([k p]) = b([p k]);
            end
        end
    % =====================================================================
% =========Setting Matrix A to triangular form=================================
    for i = k+1:m
        l(i,k) = A(i,k)/A(k,k);
        for j = k+1:n
            A(i,j) = A(i,j)-l(i,k)*A(k,j);
        end
        b(i) = b(i)-l(i,k)*b(k);
    end
end

for k = 1:m-1
    for i = k+1:m
        A(i,k) = 0;
    end
end

%==========================================================================
% Backward substitution of Matrix A
%==========================================================================
x(m) = b(m)/A(m,m);

for i = m-1:-1:1
    sum = 0;
    for j = i+1:m
        sum = sum + A(i,j)*x(j);
    end
    x(i) = (b(i)- sum)/A(i,i);
end

%==========================================================================
end

输出:

ans =

    -3
     2