解线性方程组——直接解法：(Gauss)高斯消去法、列主元、全主元

解线性方程组——直接解法：(Gauss)高斯消去法、列主元、全主元 | 北太天元

一、问题描述

对于线性方程组
$Ax=b,\quad A=\begin{pmatrix} a_{11} & a_{12} &\cdots &a_{1n}\\ a_{21} & a_{22} &\cdots &a_{2n}\\ \vdots & \vdots & &\vdots\\ a_{n1} & a_{n2} &\cdots &a_{nn}\\ \end{pmatrix},\quad b=\begin{pmatrix} b_1\\b_2\\ \vdots\\ b_n \end{pmatrix}$

A为非奇异矩阵，求向量 $x$

二、Gauss消去法

思路：

系数矩阵逐步消元得到上三角矩阵，然后使用回代算法求解。

1. (耿直版)消元法

$\begin{cases} a_{11}^{(1)}x_1+a_{12}^{(1)}x_2+...+a_{1n}^{(1)}x_n=b_1^{(1)}\\ a_{21}^{(1)}x_1+a_{22}^{(1)}x_2+...+a_{2n}^{(1)}x_n=b_2^{(1)}\\\cdots \\ a_{n1}^{(1)}x_1+a_{n2}^{(1)}x_2+...+a_{nn}^{(1)}x_n=b_n^{(1)} \end{cases}$

简记为 $A^{(1)}x=b^{(1)}$ ，其中 $A^{(1)}=A,\quad b^{(1)}=b$

若 $a_{11}^{(1)}\neq0$ ，令 $l_{i1}=\dfrac{a_{i1}^{(1)}}{a_{11}^{(1)}}$ ，

令 $l_{i1}$ 乘第 $1$ 行加到第 $i$ 行 $(i=2,3,\cdots,n)$ 并保留第 $1$ 个方程。

此时 $A^{(1)}\rightarrow A^{(2)},b^{(1)}\rightarrow b^{(2)}$
$a_{ij}^{(2)}=a_{ij}^{(1)}-l_{i1}a_{1j}^{(1)}(i,j=2,3,\cdots,n)\\ \ \\ b_{i}^{(2)}=b_{i}^{(1)}-l_{i1}b_{1}^{(1)}(i=2,3,\cdots,n)$

若 $a_{kk}^{(k)}\neq0$ ，令 $l_{ik}=\dfrac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$ ，
令 $l_{ik}$ 乘第 $k$ 行加到第 $i$ 行 $(i=k+1,k+2,\cdots,n)$
$a_{ij}^{(k+1)}=a_{ij}^{(k)}-l_{ik}a_{kj}^{(k)}(i,j=k+1,k+2,\cdots,n)$ $b_{i}^{(k+1)}=b_{i}^{(k)}-l_{ik}b_{k}^{(k)}(i=k+1,k+2,\cdots,n)$
在这里插入图片描述

缺点

$a_{kk}^{(k)}\neq0$ 若为0，不能使用

$\left|a_{kk}^{(k)}\right|\approx0$ 或相对其他元素较小时，舍入误差增加，可能导致误差剧增

进行每一步消元前，先找到列中最大元素所在的那一行，并与原来的行交换位置

2. 列主元

进行每一步消元前，先找到列中最大元素所在的那一行，并与原来的行交换位置

3. 全主元

进行每一步消元前，先找到剩余矩阵中最大的元素，做对应的行交换与列交换，
由于进行了列交换，最后要把 $X$ 的位置对应回去

三、北太天元源代码

耿直版

function [X,Ae,be] = gsem_base(A,b)
% Gauss消去法耿直版
% A : 系数矩阵
% b : 右端常数 列向量
% X : 求得的解向量
% Ae: 得到的上三角矩阵
% be: 对应的右端项
% 使用要求，A1(k,k)在计算过程中不出现0
% 缺点：当A1(k,k)较小或近似于0时，计算误差可能增大很多
%
%   Version:            1.0
%   last modified:      09/09/2023
A1 = [A,b];%  A和b的增广矩阵
n = length(A);
n1 = length(A1);
% 系数矩阵 化为上三角的过程如下
t = 0;
for i = 1:n-1if A1(i,i) != 0% 对每一行进行操作for k = i+1:1:nlik = A1(k,i)/A1(i,i);    % 算子A1(k,i:n1) = A1(k,i:n1)-lik*A1(i,i:n1);
%            T=A1 % 显示运算步骤endelset = 1;fprintf('出现零元素，不能使用该函数');break;end
end
% 得到上三角后 进行回代
if t == 0
Ae = A1(:,1:n);
be = A1(:,n1);
X  = reg_utm(Ae,be);
end
end

保存为 gsem_base.m文件

列主元

function [X,Ae,be] = gsem_column(A,b)
% 列主元消去法
% A : 表示系数矩阵
% b : 表示右端常数 [列向量]
% X : 求得的解向量
% Ae: 得到的上三角矩阵
% be: 对应的右端项
% 避免了Gauss消去法中出现0的情况
%
%   Version:            1.0
%   last modified:      09/09/2023
A1 = [A,b];%  A和b的增广矩阵
n = length(A);
n1 = length(A1);
% 系数矩阵 化为上三角的过程如下
for i = 1:n-1%[a,s]= max(A1(i:n,i))[~,s]= max(A1(i:n,i)); % s表示对应的位置,由于a用不到，用~替代s = s+i-1; % 由于s是相A对位置，故有此步A1([i,s],:) = A1([s,i],:); % 对应行交换位置% 对每一行进行操作for k = i+1:1:nlik = A1(k,i)/A1(i,i);    % 算子A1(k,i:n1) = A1(k,i:n1)-lik*A1(i,i:n1);
%        T=A1 % 显示运算步骤end
end
% 得到上三角后 进行回代
Ae = A1(:,1:n);
be = A1(:,n1);
X = reg_utm(Ae,be);
end

保存为 gsem_column.m文件

全主元

写全主元时,
注意,进行列的交换后, 对应解的位置发生了互换,最后要换回来.

function [X,Ae,be] = gsem_complete(A,b)
% 全主元消去法
% A : 系数矩阵
% b : 右端常数 [列向量]
% X : 求得的解向量
% Ae: 得到的上三角矩阵
% be: 对应的右端项
% 比列主元更加稳定了
%
%   Version:            1.0
%   last modified:      09/09/2023
A1 = [A,b];    %A和b的增广矩阵
n = length(A);
n1 = length(A1);
X = zeros(n,1);
t = 1:n; % t 用于修正解的位置与原方程不匹配
% 系数矩阵化为上三角的过程如下
for i = 1:n-1A = A1(i:n,i:n);[m,a,b] = max_loc(A); %这里不用A1，避免误判b中元素% a,b是相对的位置a = a+i-1;b = b+i-1;A1([i,a],:) = A1([a,i],:); % 对应行交换位置A1(:,[i,b]) = A1(:,[b,i]); % 对应列交换位置[t(i),t(b)]= deal(t(b),t(i));% 对每一行进行操作for k = i+1:1:nlik = A1(k,i)/A1(i,i);    % 算子A1(k,i:n1) = A1(k,i:n1)-lik*A1(i,i:n1);
%        T=A1 % 显示运算步骤end
end
% 得到上三角后进行回代
Ae = A1(:,1:n);
be = A1(:,n1);
X([t])= reg_utm(Ae,be); %使解匹配原方程位置
end

保存为gsem_complete.m文件
其中用到了我写的 max_loc.m

function [m,a,b] =max_loc(A)
% 获取矩阵中最大元素的位置
% m：矩阵中最大元素的值
% [a,b] 最大值在矩阵中的位置
%
%   Version:            1.0
%   last modified:      05/16/2023[r,l] = size(A);if r == 1    % 行向量a = 1;[m,b] = max(A);elseif l == 1    % 列向量b = 1;[m,a] = max(A);else[M,I] = max(A);[m,b] = max(M); % b = 最大元素所在的列a = I(b); % a = 最大元素所在的行end
end

保存为max_loc.m文件

四、简单测试一下

% 线性方程组的例子
%   last modified:      09/09/2023
clc;clear all;
%%
A = [10 1 2 3 4; 1 9 -1 2 -3; 2 -1 7 3 -5; 3 2 3 12 -1; 4 2 3 4 15];
b = [12; -27; 14; -17; 10];X2 = gsem_column(A,b)
X1 = gsem_base(A,b)
X3 = gsem_complete(A,b)