MATLAB 高斯消元法求解线性方程组：代码解析与步骤

本文将详细解析一段使用 MATLAB 高斯消元法求解线性方程组的代码，并解释代码中每个步骤的含义。

代码

clear; 
clc;

% 定义系数矩阵 A 和常数向量 b
A = [10 -7 0 1; -3 2.099999 6 2; 5 -1 5 -1; 2 1 0 2];
b = [8; 5.900001; 5; 1];

% 获取系数矩阵 A 的行数和列数，并初始化未知数向量 x 和置换向量 P
[n, n] = size(A); 
x = zeros(n, 1); 
P = [1:n];

% 构造增广矩阵 Aug
Aug = [A, b];

% 进行高斯消元法计算
for k = 1:n-1
    [piv, r] = max(abs(Aug(k:n, k))); % 找列主元所在子矩阵的行 r
    r = r + k - 1; % 列主元所在大矩阵的行
    if r > k
        Aug([k, r], :) = Aug([r, k], :); 
        P([k, r]) = P([r, k]);
    end
    if Aug(k, k) == 0
        error('对角元出现 0');
    end
    % 把增广矩阵消元成为上三角
    for p = k+1:n
        Aug(p, :) = Aug(p, :) - Aug(k, :) * Aug(p, k) / Aug(k, k);
    end
end

% 解上三角方程组
A = Aug(:, 1:n); 
b = Aug(:, n+1);
x(n) = b(n) / A(n, n);
for k = n-1:-1:1
    x(k) = (b(k) - A(k, n:-1:k+1) * x(n:-1:k+1)) / A(k, k);
end

% 输出置换向量 P 和解向量 x，以及系数矩阵 A 的行列式 detA
P
x
detA = det(A)

代码解析

清空命令窗口和工作区
- clear;: 清空工作区，删除所有变量。
- clc;: 清空命令窗口，删除所有显示内容。
定义系数矩阵 A 和常数向量 b
- 系数矩阵 A 代表线性方程组的系数，常数向量 b 代表方程组的常数项。
获取系数矩阵 A 的行数和列数，并初始化未知数向量 x 和置换向量 P
- [n, n] = size(A);: 获取系数矩阵 A 的行数和列数，并将其分别赋值给变量 n。
- x = zeros(n, 1);: 初始化未知数向量 x 为一个 n 行 1 列的零向量。
- P = [1:n];: 初始化置换向量 P 为一个从 1 到 n 的行向量，表示初始的行号顺序。
构造增广矩阵 Aug
- Aug = [A, b];: 将系数矩阵 A 和常数向量 b 横向拼接，形成增广矩阵 Aug。
进行高斯消元法计算
- 循环语句 for k = 1:n-1 控制消元过程，从第一列开始逐列进行消元操作。
- [piv, r] = max(abs(Aug(k:n, k)));: 找出从当前行 k 到最后一行 n 中，第 k 列元素绝对值最大的元素，并将其所在的行号记录在变量 r 中。
- r = r + k - 1;: 将子矩阵的行号 r 转换为整个矩阵 Aug 中的行号。
- if r > k: 如果列主元所在的行号 r 大于当前行号 k，则交换两行，使列主元位于当前行。
- Aug([k, r], :) = Aug([r, k], :);: 交换增广矩阵 Aug 中的第 k 行和第 r 行。
- P([k, r]) = P([r, k]);: 交换置换向量 P 中的第 k 个元素和第 r 个元素，记录行交换操作。
- if Aug(k, k) == 0: 判断当前行对角元是否为 0。如果为 0，则无法进行消元操作，抛出错误信息。
- 循环语句 for p = k+1:n: 从下一行 k+1 开始，逐行进行消元操作，将当前行 k 下方的所有元素消为 0。
- Aug(p, :) = Aug(p, :) - Aug(k, :) * Aug(p, k) / Aug(k, k);: 将第 p 行减去第 k 行乘以一个系数，将第 p 行第 k 列元素消为 0。
解上三角方程组
- A = Aug(:, 1:n);: 将增广矩阵 Aug 的前 n 列提取出来，作为系数矩阵 A。
- b = Aug(:, n+1);: 将增广矩阵 Aug 的第 n+1 列提取出来，作为常数向量 b。
- x(n) = b(n) / A(n, n);: 使用最后一个方程，直接求解最后一个未知数 x(n)。
- 循环语句 for k = n-1:-1:1: 从倒数第二个方程开始，逐个求解未知数，直到第一个方程。
- x(k) = (b(k) - A(k, n:-1:k+1) * x(n:-1:k+1)) / A(k, k);: 利用当前方程和已求解的未知数，求解当前未知数 x(k)。
输出置换向量 P 和解向量 x，以及系数矩阵 A 的行列式 detA
- P: 输出置换向量 P，记录了行交换操作的顺序。
- x: 输出解向量 x，包含了线性方程组的解。
- detA = det(A): 计算系数矩阵 A 的行列式，并输出结果。

总结

这段代码使用 MATLAB 实现高斯消元法求解线性方程组，通过逐行解析代码，我们能够理解代码的具体实现过程，并掌握高斯消元法的基本原理和步骤。