题目链接: 洛谷P3389高斯消元【模板】
提醒:在cnblog博客或是luogu博客看此文更优哦,戳这里
现在有一个三元方程:
2x+y+z=1;
6x+2y+z=-1;
-2x+2y+z=7;
很显然如果你有兴趣解这个方程组的解为
x=-1.00
y=2.00
z=1.00
如果这道题出在数学试卷上相信你一定会做,三元一次方程组嘛! 那么这里就写个程序模拟你的运算过程.
为了书写方便,这里化为矩阵的写法:
x y z val
line1: 2 1 1 1
line2: 6 2 1 -1
line3:-2 2 1 7
用等式1的x消去等式2、3的第一个x
x y z val
line1: 2 1 1 7
line2: 0 -1 -2 -4
line3: -2 2 1 7(line1*(-3)+line2)
-->
x y z val
line1: 2 1 1 7
line2: 0 -1 -2 -4
line3: 0 3 2 8 (line1*1+line3)
接下来我们拿line2第二个数消掉下面的一个数
x y z val
line1: 2 1 1 7
line2: 0 -1 -2 -4
line3: 0 0 -4 -4 (line2*3+line3)
这样我们就得到一个左下角一块全是0的矩阵,
我们将之称之为上三角矩阵 就是这个东东
x y z
line1: 2 1 1
line2: 0 -1 -2
line3: 0 0 -4
(去掉val)
那么我们还是看这个矩阵
x y z val
line1: 2 1 1 7
line2: 0 -1 -2 -4
line3: 0 0 -4 -4
(加上val)
- 对于line3我可以轻易的知道z=(-4)/(-4)=1.00
- 那么知道z的值那么可以推出line2中的y,y=(-4-(-2)*1)/(-1)=2其实就是把z=1代入line2啦
- 知道y,z的值代入line我们就可以计算出x=-1啦(方法和上面一样)
An another question:
怎么判断是不是无解或是多组解呢?
对于处理处的上三角矩阵,如果全是0(含val)那么多组解
如果有一行全是0(不含val)但是val的值不为0那么无解
这里有个问题必须要说明,对于主元为0的情况因为除数为0,可能存在错误,我们采用交换行的方式避免错误,如
HACK!数据 Ghastlcon 2018-07-17 11:46
HACK.in
3
1 1 1 3
1 1 2 4
3 2 2 7
HCAK.out
1.00
1.00
1.00
这组数据,我们尝试分析:
x y z val
line1 1 1 1 3
line2 1 1 2 4
line3 3 2 2 7
--->
x y z val
line1 1 1 1 3
line2 0 0 -1 -1 (line1-line2)
line3 0 -1 -1 -2 (line1*(-3)+line3)
我们发现line2和line3位置颠倒,换一下刚好是上三角矩阵,那就换下吧
x y z val
line1 1 1 1 3
line2 0 -1 -1 -2
line3 0 0 -1 -1
然后按照上述方案进行处理 如果是多元的从上倒下依次交换像插入排序一样,就可以保证最后处理出的矩阵是上三角矩阵。
那么我们可以模拟上述过程,写出一个高斯消元的模板:
# include <bits/stdc++.h> using namespace std; const int MAXN=105; double eps=0.000001; //精度推荐1e-4到1e-10之间 double a[MAXN][MAXN],val[MAXN]; int n; bool flag; void swapline(int i)//交换i行和下面的行,保证i行及以上有序 { int j=i+1; while (fabs(a[j][i]) < eps && j<=n) ++j; if (j>n) { printf("No Solution "); flag=true; return; } if (j<=n) {swap(a[i],a[j]);swap(val[i],val[j]);} } int main() { scanf("%d",&n); for (int i=1;i<=n;i++) { for (int j=1;j<=n;j++) scanf("%lf",&a[i][j]); scanf("%lf",&val[i]); } for (int i=1;i<=n-1;i++) { double num=a[i][i]; if (fabs(num)<eps) swapline(i); if (flag) return 0; for (int j=i+1;j<=n;j++) { if (fabs(a[j][i])<=eps) continue;//当此行主元a[i][i]=0时交换行 double w=-a[j][i]/num; a[j][i]=0; for (int k=i+1;k<=n;k++) a[j][k]+=a[i][k]*w; val[j]+=val[i]*w; } }//求上三角矩阵 for (int i=1;i<=n;i++) { flag=true; for (int j=1;j<=n;j++) if (fabs(a[i][j])>eps) { flag=false;break; } if (flag) { printf("No Solution "); return 0; //判断是否无解或有多组解 } } val[n]=val[n]/a[n][n];//求出第n个元 for (int i=n-1;i>=1;i--) { double sum=0; for (int j=n;j>=i+1;j--) sum=sum+a[i][j]*val[j]; val[i]=(val[i]-sum)/a[i][i]; }//代入已知元求出未知的一个元 for (int i=1;i<=n;i++) printf("%.2lf ",val[i]); return 0; }