Part 0. 引子
小学,我们学习了一元一次方程。其标准形式如下:
ax=b
- 当 a=0 时方程有唯一解。
- 当 a=0,b=0 时方程无解。
- 当 a=0,b=0 时方程有无穷解。
到了初中,我们又学习了二元一次方程组。其标准形式如下:
{a1x+b1y=c1a2x+b2y=c2
将两个式子化作两个一次函数,则方程解情况如下:
- 当两直线相交的时候,方程有唯一解。
- 当两直线平行的时候,方程无解。
- 当两直线重合的时候,方程有无数解。
我们是怎么解这个方程的?
加减消元,带入消元。
说白了,高斯消元就是这样解 n 元一次方程。
Part 1. 高斯-约当消元法
首先,我们假设有一个方程组长这个样子:
⎩⎪⎪⎪⎨⎪⎪⎪⎧a11x1+a12x2+...+a1nxn=c1a21x1+a22x2+...+a2nxn=c2...an1x1+an2x2+...+annxn=cn
要想求解怎么办呢?
说实话,如果手算,这个问题并不困难。我们利用加减消元,尝试挨个消元,(后面的过程我也没法描述),就差不多了。
然而我们要用计算机解决问题,也就是说我们需要找到一种通用解法。
其中一种解法是求根公式。基本上所有可以解的方程都可以求出求根公式,但是对于这样一个勾诗东西我们肯定不愿意算这么一个硕大无比的求根公式——而且你还得把他敲到代码上。
那么我们就类比。类比啥呢?我们解二元一次方程组的时候用到了两种方式:加减消元,代数消元。
- 代数消元。这种方式需要你先利用一个方程把一个未知数表示出来,然后再带进所有的方程里面。想想都头大。
- 加减消元。这种方式好像直接化成一般形式开算就完了。嗯,很不错。
于是我们就首推加减消元。
为了简化计算,我们把方程组变成矩阵的样子:
⎣⎢⎢⎡a11⋮am1⋯⋱⋯a1n⋮amnb1⋮bn⎦⎥⎥⎤
这个东西正好对应方程的系数。
现在我们考虑这个玩意怎么才能得到方程的解。显然,当一个矩阵的某一行是这个样子的时候:
[0⋯010⋯0b]
这个方程就能确定一个正整数解。
我们只需要对这些矩阵进行初等变换就可以得到这样的结果。好了,现在我们的目的就变成了把一个矩阵变成只有对角线上是 1,其他地方全是 0 的矩阵,然后最后一列就是方程的解。
大致上,我们可以把对求解每一列的过程分成三步:
- 找到这一列的(非零的)最小值。把这一行的每个数都除以这个数,使得这个最小值所在位置变成 1。
- 把每一行都除以当前列的数(如果这一列是 0 则跳过)
- 用这一行减去最小值那一列。此时除了最小值所在行,其余行当前列的数都变成 0。
重复这个操作每一列,得到答案。
不过……等等。方程一定有唯一解吗?答案是否定的。除非题目有描述,最好判断这种情况。
- 每一行都有一个 1 和一车 0。完美。
- 全是 0,但最后一个数却不是 0。别告诉我 0 等于一个非零的数。无解。
- 整一行找不出一个非零的数。好吧,你随便代数,反正能满足。
当然了,这个方程可能不是整数解,所以要注意浮……点……数……精……度……问……题……
Part 2. Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
| #include<iostream> #include<cstdio> #include<cmath> using namespace std; typedef double db; const int N=110; const db eps=1e-7; db a[N][N]; int n;
void gauss(){ for(int i=1;i<=n;i++){ int mx=0; for(int j=i;j<=n;j++){ if(fabs(a[i][j])>fabs(a[i][mx]))mx=j; } if(mx!=i)swap(a[i],a[mx]); if(fabs(a[i][i])<eps)continue; for(int j=n+1;j>=i;j--)a[i][j]/=a[i][i]; for(int j=1;j<=n;j++){ if(j==i)continue; for(int k=i+1;k<=n+1;k++)a[j][k]-=a[i][k]*a[j][i]; } } } int main(){ cin>>n; for(int i=1;i<=n;i++){ for(int j=1;j<=n+1;j++)cin>>a[i][j]; } gauss(); for(int i=1;i<=n;i++)if(fabs(a[i][i])<eps){cout<<"No Solution"<<endl;return 0;} for(int i=1;i<=n;i++)printf("%.2f\n",a[i][n+1]); }
|
Part 3. 例题
咕咕咕。