\(Q:\)高斯消元是什么?听起来好高级啊??
\(A:\)二元一次方程组解过吗?那就是高斯消元。
- 高斯消元(\(Gauss\)),是一种用来求解线性方程组的方法,在\(OI\)竞赛中广泛使用。
首先对高斯消元做一些准备:
\(Q:\)什么是线性方程组?
\(A:\)鸡兔同笼方程组 由\(M\)个\(N\)元一次方程所构成的方程组。
为了简化表达,做出如下定义:
系数矩阵
把线性方程组未知数的系数写成一个\(M\times N\)的矩阵,即为“系数矩阵”。
增广矩阵
把线性方程组等号右边的常数写成一个\(M\times 1\)的矩阵,拼在系数矩阵右边,即为一个大小为\(M\times (N+1)\)的“增广矩阵”。
\(Example:\)
假如现有一方程组:
\[\begin{cases}x_1+x_2+4x_3=17\\5x_1+x_2+4x_3=25\\8x_1+x_2+0x_3=19\end{cases}\]
其系数矩阵即为:
\[\begin{bmatrix}1&1&4\\5&1&4\\8&1&0\end{bmatrix}\]
(???)
增广矩阵为:
\[\begin{bmatrix}1&1&4&17\\5&1&4&25\\8&1&0&19\end{bmatrix}\]
\(Q:\)这些很简单啊,接下来要干什么呢?
\(A:\)接着我们对方程组进行求解,包括如下\(3\)个步骤:
对其中一个方程乘以一个非\(0\)常数
用一个方程加上另一个方程
交换两个方程
称以上操作为矩阵的“初等行变换”。
\(Q:\)这些毫无意义的语言 都是什么意思啊?
以一个二元一次方程组为例:
\[\begin{cases}x_1+9x_2=21\quad(1)\\9x_1+x_2=29\quad(2)\end{cases}\]
增广矩阵:
\[\begin{bmatrix}1&9&21\\9&1&29\end{bmatrix}\]
- 对其中一个方程乘以一个非\(0\)常数:
\((1)*=-9:\)
\[\begin{cases}-9x_1+-81x_2&=-189\quad(1)\\9x_1+x_2&=29\quad(2)\end{cases}\]
增广矩阵:
\[\begin{bmatrix}-9&-81&-189\\9&1&29\end{bmatrix}\]
- 用一个方程加上另一个方程
\((1)+=(2):\)
\[\begin{cases}-80x_2&=-160\\9x_1+x_2&=29\end{cases}\]
\[\begin{cases}x_2&=2\\9x_1+x_2&=29\end{cases}\]
增广矩阵:
\[\begin{bmatrix}0&1&2\\9&1&29\end{bmatrix}\]
- 交换两个方程
\(Swap((1),(2)):\)
\[\begin{cases}9x_1+x_2&=29\\x_2&=2\end{cases}\]
增广矩阵:
\[\begin{bmatrix}9&1&29\\0&1&2\end{bmatrix}\]
怎么样,是不是很简单啊?
初等行变换后的矩阵被称为“阶梯型矩阵”(是不是很像?并没有)。
系数部分称为“上三角矩阵”。
最后从下到上更新一遍即可求出方程组的解。
另外,若在求解过程中出现了\(0=a\)的情况,则方程组无解。
若过程中无法消元(找不到一个方程,\(x_1\sim x_{i-1}\)系数为\(0\),\(x_i\)系数不为\(0\)),则解不唯一,称此类\(x_i\)为自由元。
时间复杂度 \(O(n^3)\)
空间复杂度 \(O(n^2)\)
接下来是实际应用。
这就是裸的模板题了。
时间复杂度 \(O(n^3)\)
空间复杂度 \(O(n^2)\)
代码:
#include <cstdio>
#include <algorithm>
inline double Abs(double x){return x>=0?x:-x;}
int n;
double a[105][105];
const double Eps=1e-7;//精度视情况而定,不要太小或太大
int main()
{
scanf("%d",&n);
for(int i=1;i<=n;++i)
for(int j=1;j<=n+1;++j)
scanf("%lf",&a[i][j]);//共同构成增广矩阵
for(int i=1;i<=n;++i)//进行n次消元,消去第i项未知数
{
for(int j=1;j<=n;++j)//找一个x[i]系数不为0的方程进行消元
if(Abs(a[j][i])>Eps)
for(int k=1;k<=n+1;++k)
std::swap(a[i][k],a[j][k]);
if(Abs(a[i][i])<=Eps)return puts("No Solution"),0;//出现自由元,无解
//消去其他方程x[i]系数
for(int j=1;j<=n;++j)
if(i!=j)
{
double Rate=a[j][i]/a[i][i];//相对比例
for(int k=1;k<=n+1;++k)a[j][k]-=a[i][k]*Rate;//消元
}
}
for(int i=1;i<=n;++i)printf("%.2f\n",a[i][n+1]/a[i][i]);
return 0;
}
稍微转换一下题目
设已知点为\(a_i\),求一个点\((x_1,x_2,\dots,x_n)\),使得存在常数\(Dis\),满足
\[\forall i\in[1,n+1],\sum_{1\le j\le n}(a_{i,j}-x_j)^2=Dis\]
这样就有了\(n+1\)个方程,但发现无法求解。
考虑构造\(n\)个\(n\)元一次方程组进行高斯消元。
用相邻的两个方程做差,得:
\[\sum_{1\le j\le n}(a_{i,j}-x_j)^2-\sum_{1\le j\le n}(a_{i+1,j}-x_j)^2=Dis-Dis\]
\[\sum_{1\le j\le n}(a_{i,j}-x_j)^2-(a_{i+1,j}-x_j)^2=0\]
\[\sum_{1\le j\le n}2(a_{i,j}-a_{i+1,j})x_j=\sum_{1\le j\le n}(a_{i,j}^2-a_{i+1,j}^2)\]
高斯消元即可。
时间复杂度 \(O(n^3)\)
空间复杂度 \(O(n^2)\)
代码:
#include <cmath>
#include <cstdio>
#include <algorithm>
const double Eps=1e-8;
int n;
double p[15][15],b[15],c[15][15];
int main()
{
scanf("%d",&n);
for(int i=1;i<=n+1;++i)
for(int j=1;j<=n;++j)
scanf("%lf",&p[i][j]);
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
{
c[i][j]=(p[i][j]-p[i+1][j])*2;
b[i]+=p[i][j]*p[i][j]-p[i+1][j]*p[i+1][j];//这里把系数和常数分开放
}
for(int i=1;i<=n;++i)
{
for(int j=i;j<=n;++j)
if(fabs(c[j][i])>Eps)
{
for(int k=1;k<=n;++k)std::swap(c[i][k],c[j][k]);
std::swap(b[i],b[j]);
}
for(int j=1;j<=n;++j)
{
if(i==j)continue;
double Rate=c[j][i]/c[i][i];
for(int k=i;k<=n;++k)c[j][k]-=c[i][k]*Rate;
b[j]-=b[i]*Rate;
}
}
for(int i=1;i<=n;++i)
printf("%.3f%c",b[i]/c[i][i],i==n?'\n':' ');
return 0;
}
难题就不继续探讨了,毕竟我也不会
参考资料:
各省省选题目
《算法竞赛进阶指南》 -- 李煜东