高斯消元 学习笔记
数学上，高斯消元法（或译：高斯消去法），是线性代数规划中的一个算法，可用来为线性方程组求解。
——百度百科
说实话，我不相信这是高斯发明的。感觉像是个小学生都学过的加减消元法。
它的时间复杂度与方程个数、未知数个数有关，一般来讲，是 \(O(n^3)\)。
概念1：线性方程组
有多个未知数，且每个未知数的次数皆为一次，这样的多个未知数所组成的方程组被称为 线性方程组。
其形式为：
概念2：增广矩阵
我们将上面线性方程组的所有系数记录在一个矩阵里，就叫此矩阵为增广矩阵。
概念3：主元
算法实现中的当前未知数称为主元。
高斯消元法：
例：
首先，进行消元操作：
将 (1)/3 ，将 \(x\) 的系数化为 \(1\)：\(x+\frac{2}{3}y+\frac{1}{3}z=2 \ \ \ (1)'\)。
再令 \((2)-2\times(1)'\) 、\((3)-4\times(1)'\)，消去 \((2),(3)\) 的 \(x\)。
得到：
然后，令 \((2)'\) 除以 \(\frac{2}{3}\) ，使得 \(y\) 的系数成为 \(1\)，得到：
\(y+2z=0 \ \ (2)''\)
再令 \((3)'-\frac{-14}{3}\times(2)''\)，使得 \(（3）\) 的 \(y\) 被消去，得到：
由 \((3)''\) 可以清晰地得到 \(z=-1\)，将其带入 \((2)''\)便可快速地解出 \(y=2\) 。
最后将二者带回 \((1)'\) 即可清晰快速地得到 \(x=1\) ，所以此线性方程组的解为：
增广矩阵模拟：
初态：
消去 \(x\)，得到：
消去 \(y\)，可以得到：
回带解得：
判断无单解情况
- 无解
消元完毕后，若有一行系数全部为 0，但最右端常数项不为 0，即无解，例：
- 多解 （ 前提：并非无解 ）
消元完毕后，若有一行系数全部为 0，但最右端常数也为 0，此时存在多解，例：
例题：
[SDOI2006] 线性方程组
题目描述
已知 \(n\) 元线性一次方程组。
\[\begin{cases} a_{1, 1} x_1 + a_{1, 2} x_2 + \cdots + a_{1, n} x_n = b_1 \\ a_{2, 1} x_1 + a_{2, 2} x_2 + \cdots + a_{2, n} x_n = b_2 \\ \cdots \\ a_{n,1} x_1 + a_{n, 2} x_2 + \cdots + a_{n, n} x_n = b_n \end{cases} \]请根据输入的数据，编程输出方程组的解的情况。
输入格式
第一行输入未知数的个数 \(n\)。
接下来 \(n\) 行，每行 \(n + 1\) 个整数，表示每一个方程的系数及方程右边的值。输出格式
如果有唯一解，则输出解（小数点后保留两位小数）。
如果方程组无解输出 \(-1\)；
如果有无穷多实数解，输出 \(0\)；
算法实现：
首先，将输入数据存成 double 类型的增广矩阵，然后进行模拟：
-
枚举两个变量 \(r,c\) 分别表示当前行和当前主元，找到主元系数不为 0 的第 \(t\) 行
-
交换当前行和第 \(t\) 行
-
把当前行主元系数变为 1
-
把从第 \(r+1\) 行道第 \(n\) 行的当前主元的系数变为 0
-
最后，从第 \(n\)行开始往上解出答案。
-
若循环完后的 \(r<=n\) ，则不存在单解，具体请在代码中看。
#include<bits/stdc++.h>
using namespace std;
#define int long long
int n;
int f;
double a[101][101];
signed main()
{
cin>>n;
for(int i=0;i<n;i++)
{
for(int j=0;j<=n;j++)
{
cin>>a[i][j];
}
}
int r=0,c=0;
for(;c<n;c++)
{
int t=r;
for(int i=r;i<n;i++)
{
if(fabs(a[i][c])>fabs(a[t][c]))
{
t=i;
}
}
if(fabs(a[t][c])<1e-6)
{
continue;
}
for(int i=c;i<=n;i++)
swap(a[r][i],a[t][i]);
/* 主元变为 1 */
for(int i=n;i>=c;i--)
{
a[r][i]/=a[r][c];
}
/* 清理后面的子矩阵 */
for(int i=r+1;i<n;i++)
{
if(fabs(a[i][c])>1e-6)
for(int j=n;j>=c;j--)
{
a[i][j]-=a[i][c]*a[r][j];
}
}
r++;
}
/* 回带三角矩阵求解 */
for(int i=n-1;i>=0;i--)
{
for(int j=i+1;j<n;j++)
{
a[i][n]-=a[i][j]*a[j][n];
// a[i][j]=0;
}
// a[i][i]=1;
}
/* 判断不存在单解的情况 */
if(r<n)
{
for(int i=r;i<n;i++)
{
if(fabs(a[i][n])>1e-6)
{
cout<<-1;return 0;
}
}
cout<<0;return 0;
}
for(int i=0;i<n;i++)
{
printf("x%d=%.2lf\n",i+1,a[i][n]);
}
}
