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

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

文章对应视频讲解：(Gauss)高斯消去法、列主元、全主元

一、问题描述

对于线性方程组
$$Ax=b,A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right),b=\left(\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right)$$

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

二、Gauss消去法

思路：

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

1. (耿直版)消元法

$$\{\begin{array}{l}{a}_{11}^{(1)}{x}_1+a_{12}^{(1)}{x}_2+\dots +a_{1n}^{(1)}{x}_n=b_1^{(1)}\\ a_{21}^{(1)}{x}_1+a_{22}^{(1)}{x}_2+\dots +a_{2n}^{(1)}{x}_n=b_2^{(1)}\\ \cdots \\ a_{n1}^{(1)}{x}_1+a_{n2}^{(1)}{x}_2+\dots +a_{nn}^{(1)}{x}_n=b_n^{(1)}\end{array}$$

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

若$a_{11}^{(1)}\ne 0$，令$l_{i1}=\frac{a_{i1}^{(1)}}{a_{11}^{(1)}}$，

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

此时 $A^{(1)}\to A^{(2)},b^{(1)}\to b^{(2)}$
$$\begin{gathered} 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) \end{gathered}$$

若$a_{kk}^{(k)}\ne 0$，令$l_{ik}=\frac{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)}\ne 0$ 若为0，不能使用

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

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

2. 列主元

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

3. 全主元

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

三、北太天元 or matlab 实现

耿直版

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-1 if A1(i,i) != 0 % 对每一行进行操作 for k = i+1:1:n lik = A1(k,i)/A1(i,i); % 算子 A1(k,i:n1) = A1(k,i:n1)-lik*A1(i,i:n1); % T=A1 % 显示运算步骤 end else t = 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:n lik = 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-1 A = 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:n lik = 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)