高斯消元

高斯消元 定义

以下定义以方程组为例子：

\[\begin{cases} \notag 2*x_1+3*x_2+5*x_3=33\\\notag 9*x_1+2*x_2+7*x_3=18\\\notag 9*x_1+8*x_2+4*x_3=43 \end{cases} \]

系数矩阵

定义：将方程组的系数一一对应在矩阵内。

\[\begin{bmatrix} \notag 2&3&5\\\notag 9&2&7\\\notag 9&8&4 \end{bmatrix} \]

增广矩阵

定义：在系数矩阵的右侧加上一列，这一列是方程组等号右边的值。

\[\begin{bmatrix} \notag 2&3&5&33\\\notag 9&2&7&18\\\notag 9&8&4&43 \end{bmatrix} \]

行阶梯形矩阵:

定义：元素不全为零的行(非零行）的第一个非零元素所在列的下标随着行标的增大而严格增大(列标一定不小于行标)﹔元素全为零的行(如果有的话)必在矩阵的最下面几行。

即形如：

\[\begin{bmatrix} \notag 1&2&3&4\\\notag 0&5&6&7\\\notag 0&0&8&9\\\notag 0&0&0&10 \end{bmatrix} \begin{bmatrix} \notag 1&2&3&4\\\notag 0&5&6&7\\\notag 0&0&0&9\\\notag 0&0&0&0 \end{bmatrix} \]

矩阵的秩

定义：在线性代数中，一个矩阵\(A\)的列秩是\(A\)的线性独立的纵列的极大数，通常表示为$ r(A),rk(A)或rank A $。矩阵的行秩，列秩，秩都相等，等变换不改变矩阵的秩。

解释以下线性独立：

\[\begin{cases} 2*x_1+3*x_2=0\\\notag 4*x_1+6*x_2=3 \end{cases} \]

对于这个方程组，若 \((1)(2)\) 相消会消掉两个及其以上个未知数，则这两个方程线性不独立。

反之如：

\[\begin{cases} 2*x_1+3*x_2=0\\\notag 4*x_1+7*x_2=3 \end{cases} \]

这两个方程则线性独立。

当然这只是简易理解，若想详细了解请点击

可逆方阵

定义：矩阵\(A\)为\(n\)阶方阵，若存在\(n\)阶矩阵\(B\)，使得矩阵\(A,B\)的乘积为单位方阵，则称\(A\)为可逆阵，\(B\)为\(A\)的逆矩阵。若方阵的逆阵存在，则称为可逆矩阵或非奇异矩阵，且其逆矩阵唯一。即$ AB=I $

高斯消元 基本步骤

解方程组：

\[\begin{align} 2x+y-z&=8\\ -3x-y+2z&=-11\\ -2x+y+2z&=-3 \end{align} \]

可以将其增广矩阵写出：

\[\begin{bmatrix} \notag 2 &1 &-1 &8 \\ \notag -3 &-1 &2 &-11\\ \notag -2 &1 &2 &-3 \end{bmatrix} \]

将其转化为形如这样的行阶梯形矩阵：

\[\begin{bmatrix} \notag a&b&c&d\\ 0&e&f&g\\ 0&0&h&i \end{bmatrix} \]

转换步骤：

1. 选取第\(i\)列，取其绝对值最大的一行与第\(i\)行交换：

   绝对值最大的为第二排 \(-3\)，将第一行与第二行交换。

   \[\begin{bmatrix} \notag 2 &1 &-1 &8 \\ \notag -3 &-1 &2 &-11\\ \notag -2 &1 &2 &-3 \end{bmatrix} \iff \begin{bmatrix} \notag -3 &-1 &2 &-11\\ \notag 2 &1 &-1 &8\\ \notag -2 &1 &2 &-3 \end{bmatrix} \]

2. 第\(i\)列第\(i\)行之下的数变为\(0\):

   \[\begin{bmatrix} \notag -3 &-1 &2 &-11\\ \notag 2 &1 &-1 &8\\ \notag -2 &1 &2 &-3 \end{bmatrix} \iff \begin{bmatrix} \notag -3 &-1 &2 &-11\\ \notag 0 &\frac{1}{3} &\frac{1}{3} &\frac{2}{3}\\ \notag 0 &\frac{5}{3} &\frac{2}{3} &\frac{13}{3} \end{bmatrix} \]

   这个就像是线性方程组消元。

3. 重复以上过程直到\(i=n\)

   得到矩阵：

   \[\begin{bmatrix} \notag -3 &-1 &2 &-11\\ \notag 0 &\frac{5}{3} &\frac{2}{3} &\frac{13}{3}\\ \notag 0&0&\frac{1}{5}&-\frac{1}{5} \end{bmatrix} \]

   此时可以看作得到方程组：

   \[\begin{cases} \notag -3x-y+2z=-11\\ \frac{5}{3}y+\frac{2}{3}z=\frac{13}{3}\\ \frac{1}{5}z=-\frac{1}{5} \end{cases} \]

   所以带回可以得到：

   \[\begin{cases} \notag x=2\\ y=3\\ z=-1\\ \end{cases} \]

   经检验，正确。

   以上就是高斯消元的基本步骤，其实不太好将清楚，读者可以通过自行手动模拟或者查询资料得到熟悉。

时间效率

\(O(n^3)\)

关于有无解，有无穷解

显然，当在选取绝对值最大的数时，若选择出来的数为\(0\)，则证明此未知数是任意数，所以无解。（说法不完全正确，仅特质题目）

有解就比较简单，解出来了就有解。

但关于无解与无穷解的判定有些复杂，将在文末呈现。

CODE

#include<bits/stdc++.h>
#define ll long long
#define ld long double 
using namespace std;
const int maxn = 103;
const ld eps=1e-6;
int int_maxn=1e9;
ll ll_maxn=1e18;
inline ll read_int(){
    ll a=0,f=0,g=getchar();
    while(g<'0'||g>'9'){if(g=='-') f=1;g=getchar();}
    while('0'<=g&&g<='9') a=a*10+g-'0',g=getchar();
    return f ? -a : a;
}

inline void write(ll s,bool f){
    int top=0,a[40];
    if(s<0) s=-s,putchar('-');
    while(s) a[++top]=s%10,s/=10;
    if(top==0) a[++top]=0;
    while(top) putchar(a[top]+'0'),top--;
    if(f) putchar('\n');
}

const string NO="No Solution";
int n;
ld a[maxn][maxn];
const ld spe=1e-6;

inline bool Gauss(){
	for(int i=1;i<=n;i++){//消第几个元 
		int max_set=i;
		for(int e=i+1;e<=n;e++) if(fabs(a[e][i])>fabs(a[max_set][i])) max_set=e;
		if(fabs(a[max_set][i])<eps) return 1;//if the max number is zero ,it means that this x can be lots of number.this equation has endless solution. 
		for(int e=i;e<=n+1;e++) swap(a[i][e],a[max_set][e]);//将最大数所在行换到第i行
		for(int e=i+1;e<=n;e++){
			ld cs=a[e][i]/a[i][i];
			for(int j=i;j<=n+1;j++){
				a[e][j]-=cs*a[i][j];//基础矩阵变换 
			}
		}
	}
	return 0;
}

inline void Print(){
	for(int i=n;i>0;i--){
		for(int j=i+1;j<=n;j++){
			a[i][n+1]-=a[j][n+1]*a[i][j];
		}
		a[i][n+1]/=a[i][i];
		if(fabs(a[i][n+1])<eps) a[i][n+1]=0;
	}
	for(int i=1;i<=n;i++) printf("%.2Lf\n",a[i][n+1]);
}

inline void read(){
	n=read_int();
	for(int i=1;i<=n;i++)for(int e=1;e<=n+1;e++)a[i][e]=read_int();
	if(Gauss()) {cout<<NO;return;}
	Print();
}

int main (){
	read();
}

请注意以下几点：

1. 关于浮点数精度，我们认为当 fabs(a)<eps时，\(a=0\)。
2. if(fabs(a[i][n+1])<eps) a[i][n+1]=0; 是为了防止输出 \(0.00\)。
3. \(No~Solution\)的\(S\)是大写。

高斯-约旦消元法

高斯消元是先把矩阵转换成行阶梯矩阵，然后回代求出每个的解，高斯-约旦消元法在高斯消元的基础上进行修改，如果方程组存在唯一解则把系数矩阵消成对角线为\(1\)，其他项为\(0\)的形式，答案为：\(x_i=a[i][n+1]\)。

高斯-约旦消元法不需要回代求解，时间效率仍然为\(O(n^3)\)，但常数稍高。

即把矩阵转化为以下形式：

\[\begin{bmatrix} \notag 1&0&0&0&x_1\\ 0&1&0&0&x_2\\ 0&0&1&0&x_3\\ 0&0&0&1&x_4 \end{bmatrix} \]

模拟步骤

解方程组：

\[\begin{align} 2x+y-z&=8\\ -3x-y+2z&=-11\\ -2x+y+2z&=-3 \end{align} \]

可以将其增广矩阵写出：

\[\begin{bmatrix} \notag 2 &1 &-1 &8 \\ \notag -3 &-1 &2 &-11\\ \notag -2 &1 &2 &-3 \end{bmatrix} \]

当处理第\(i\)个元时，选择绝对值最大的那个放在第\(i\)行：

\[\begin{bmatrix} \notag 2 &1 &-1 &8 \\ \notag -3 &-1 &2 &-11\\ \notag -2 &1 &2 &-3 \end{bmatrix} \iff \begin{bmatrix} \notag -3 &-1 &2 &-11\\ \notag 2 &1 &-1 &8\\ \notag -2 &1 &2 &-3 \end{bmatrix} \]

将第\(i\)行的第\(i\)个数变为\(1\)：

\[\begin{bmatrix} \notag -3 &-1 &2 &-11\\ \notag 2 &1 &-1 &8\\ \notag -2 &1 &2 &-3 \end{bmatrix} \iff \begin{bmatrix} \notag 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 2 &1 &-1 &8\\ \notag -2 &1 &2 &-3 \end{bmatrix} \]

使第\(i\)列除第\(i\)行外全部变为\(0\)：

\[\begin{bmatrix} \notag 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 2 &1 &-1 &8\\ \notag -2 &1 &2 &-3 \end{bmatrix} \iff \begin{bmatrix} 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 0 &\frac{1}{3} &\frac{1}{3} &\frac{2}{3}\\ \notag 0 &\frac{5}{3} &\frac{2}{3} &\frac{13}{3} \end{bmatrix} \]

重复以上步骤，直到\(i=n\)时\(x_i=a[i][n+1]\)

我们不妨在进一步模拟，当处理到第二元时：

选择最大数字：

\[\begin{bmatrix} 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 0 &\frac{1}{3} &\frac{1}{3} &\frac{2}{3}\\ \notag 0 &\frac{5}{3} &\frac{2}{3} &\frac{13}{3} \end{bmatrix} \iff \begin{bmatrix} \notag 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 0 &\frac{5}{3} &\frac{2}{3} &\frac{13}{3}\\ \notag 0 &\frac{1}{3} &\frac{1}{3} &\frac{2}{3} \end{bmatrix} \]

将\(a[i][i]\)变为\(1\)：

\[\begin{bmatrix} \notag 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 0 &\frac{5}{3} &\frac{2}{3} &\frac{13}{3}\\ \notag 0 &\frac{1}{3} &\frac{1}{3} &\frac{2}{3} \end{bmatrix} \iff \begin{bmatrix} \notag 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 0 &1 &\frac{2}{5} &\frac{13}{5}\\ \notag 0 &\frac{1}{3} &\frac{1}{3} &\frac{2}{3} \end{bmatrix} \]

使第\(i\)列除第\(i\)行外全部变为\(0\)：

\[\begin{bmatrix} \notag 1&\frac{1}{3}&-\frac{2}{3}&\frac{11}{3}\\ \notag 0 &1 &\frac{2}{5} &\frac{13}{5}\\ \notag 0 &\frac{1}{3} &\frac{1}{3} &\frac{2}{3} \end{bmatrix} \iff \begin{bmatrix} \notag 1&0&-\frac{4}{5}&\frac{14}{5}\\ \notag 0 &1 &\frac{2}{5} &\frac{13}{5}\\ \notag 0 &0 &\frac{1}{5} &-\frac{1}{5} \end{bmatrix} \]

建议多多手动模拟

CODE

#include<bits/stdc++.h>
#define ll long long
#define ld long double 
using namespace std;
const int maxn = 103;
const ld eps=1e-6;
int int_maxn=1e9;
ll ll_maxn=1e18;
inline ll read_int(){
    ll a=0,f=0,g=getchar();
    while(g<'0'||g>'9'){if(g=='-') f=1;g=getchar();}
    while('0'<=g&&g<='9') a=a*10+g-'0',g=getchar();
    return f ? -a : a;
}

inline void write(ll s,bool f){
    int top=0,a[40];
    if(s<0) s=-s,putchar('-');
    while(s) a[++top]=s%10,s/=10;
    if(top==0) a[++top]=0;
    while(top) putchar(a[top]+'0'),top--;
    if(f) putchar('\n');
}

const string NO="No Solution";
int n;
ld a[maxn][maxn];
const ld spe=1e-6;

inline bool Gauss(){
	for(int i=1;i<=n;i++){//消第几个元 
		int max_set=i;
		for(int e=i+1;e<=n;e++) if(fabs(a[e][i])>fabs(a[max_set][i])) max_set=e;
		if(fabs(a[max_set][i])<eps) return 1;
		if(max_set^i) for(int e=i;e<=n+1;e++) swap(a[i][e],a[max_set][e]);
		for(int e=n+1;e>=i;e--) a[i][e]/=a[i][i];
		for(int e=1;e<=n;e++){
			if(i==e) continue; 
			ld tmp=a[e][i];
			for(int k=i;k<=n+1;k++){
				a[e][k]-=a[i][k]*tmp;
			}
		}
	}
	return 0;
}

inline void Print(){
	for(int i=1;i<=n;i++) printf("%.2Lf\n",fabs(a[i][n+1])<eps ? 0 : a[i][n+1]);
}

inline void read(){
	n=read_int();
	for(int i=1;i<=n;i++)for(int e=1;e<=n+1;e++)a[i][e]=read_int();
	if(Gauss()) {cout<<NO;return;}
	Print();
}

int main (){
	read();
}

关于方程解的判断

这个可以通过改进高斯——约旦消元算法得到。

难点在于无解和无数解的判定，常规高斯消元法和高斯·约旦消元法遇到这两种情况直接结束消元，那如何判断无解和无穷解呢?

根据上面的结论，首先要把增广矩阵消解成上三角型的行阶梯矩阵如果存在\(a[i][i]=0\)，时\(a[i][n+1]=0\)，时就可能是无穷解，如果存在一行\(a[i][i]=0\)，时\(a[i][n+1]!=0\)时就无解，但用正常的消解法这样并不能做出准确的判断，例如存在如下的增广矩阵:

\[\begin{bmatrix} \notag 1&2&3&4&5&6&10\\ 0&0&1&2&3&4&15\\ 0&0&0&2&6&7&30\\ 0&0&0&0&6&8&10\\ 0&0&0&0&0&1&30\\ 0&0&0&0&0&0&0 \end{bmatrix} \]

此矩阵不需消解就是一个行阶梯矩阵，其中\(a[5][5]==0\)，但\(a[5][7]= 30\)，判断方程组无解，实际上无穷解，如果我们把增广矩阵第一列与第六列进行交换，这样并不影响方程的解。

\[\begin{bmatrix} \notag 6&2&3&4&5&1&10\\ 4&0&1&2&3&0&15\\ 7&0&0&2&6&0&30\\ 8&0&0&0&6&0&10\\ 1&0&0&0&0&0&30\\ 0&0&0&0&0&0&0 \end{bmatrix} \]

经过消解后变成如下矩阵（为了便于观察，对实数取整)︰

\[\begin{bmatrix} \notag 8&0&0&0&6&0&10\\ 0&2&3&4&0&1&2\\ 0&0&1&2&0&0&10\\ 0&0&0&2&1&0&21\\ 0&0&0&0&-1&0&29\\ 0&0&0&0&0&0&0 \end{bmatrix} \]

此矩阵除了\(a[6][6]=0\)外其他\(a[i][i]!=0\)，所以方程组无穷解，所以说方程的顺序会影响解的判断，那如何在不改变变量位置情况下能正确对方程组进行判断呢?我们对高斯消元进行稍作改进，消解时把坡度一致的倒三角改成坡度不一的倒三角，具体做法如下只需在消解到第\(i\)行时，如果$i $ ~ $ n\(行的第\)i$列全为\(0\)时，常规做法是枚举第\(i+1\)行消解第\(i+1\)列，此时我们改变策略，我们保持消解行仍为第\(i\)行但消解第\(i+1\)列，如果第\(i+1\)列让然找不到非\(0\)的主元，接着再第\(i\)行\(i+2\)列寻找非零主元，依次类推直到找到非零主元为止，具体操作见代码。

CODE

#include<bits/stdc++.h>
#define ll long long
#define ld long double 
using namespace std;
const int maxn = 103;
const ld eps=1e-6;
int int_maxn=1e9;
ll ll_maxn=1e18;
inline ll read_int(){
    ll a=0,f=0,g=getchar();
    while(g<'0'||g>'9'){if(g=='-') f=1;g=getchar();}
    while('0'<=g&&g<='9') a=a*10+g-'0',g=getchar();
    return f ? -a : a;
}

inline void write(ll s,bool f){
    int top=0,a[40];
    if(s<0) s=-s,putchar('-');
    while(s) a[++top]=s%10,s/=10;
    if(top==0) a[++top]=0;
    while(top) putchar(a[top]+'0'),top--;
    if(f) putchar('\n');
}

const string NO="No Solution";
int n;
ld a[maxn][maxn];
const ld spe=1e-6;

int base=1;
inline bool Gauss(){
	for(int i=1;i<=n;i++){
		int max_set=base;
		for(int e=base+1;e<=n;e++) if(fabs(a[e][i])>fabs(a[max_set][i])) max_set=e;
		if(fabs(a[max_set][i])<eps) continue;
		if(max_set^base) for(int e=i;e<=n+1;e++) swap(a[base][e],a[max_set][e]);
		for(int e=n+1;e>=i;e--) a[base][e]/=a[base][i];
		for(int e=1;e<=n;e++){
			if(e==base) continue;
			for(int k=n+1;k>=i;k--){
				a[e][k]-=a[e][i]*a[base][k];
			}
		}
		base++;
	}
	return 0;
}

inline void Print(){
	if(base<=n){
		for(int i=base;i<=n;i++){
			if(fabs(a[i][n+1])>eps) {write(-1,1);return;}
		}
		write(0,1);
		return;
	}
	for(int i=n;i>0;i--){
		for(int j=i+1;j<=n;j++){
			a[i][n+1]-=a[j][n+1]*a[i][j];
		}
		a[i][n+1]/=a[i][i];
		if(fabs(a[i][n+1])<eps) a[i][n+1]=0;
	}
	for(int i=1;i<=n;i++) printf("x%d=%.2Lf\n",i,a[i][n+1]);
}

inline void read(){
	n=read_int();
	for(int i=1;i<=n;i++)for(int e=1;e<=n+1;e++)a[i][e]=read_int();
	Gauss();
	Print();
}

int main (){
	read();
}

各位大佬们，多多手动模拟啊！