其实并没有想象中那么复杂。
就是一个简单的解方程过程。
众所周知,(n)个不等价的(n)元一次方程可以确定一组解(或者确定方程无解)。
高斯消元就是为了求得这组解。
以一个四元一次方程组
(2x_1+x_2-5x_3+2x_4=12)
(-2x_1-x_2+x_3+x_4=14)
(x_1+6x_2-3x_3+x_4=20)
(3x_1+x_2+4x_3-7x_4=21)
为例。
方便操作,我们将系数放入矩阵,常数项放在这个矩阵右侧。
(left[ egin{array}{cccc} 2 & 1 & -5 & 2\-2 & -1 & 1 & 1\1 & 6 & -3 & 1\3 & 1 & 1 & -1\ end{array} ight])(left[ egin{array}{c}12\14\20\21end{array} ight])
首先,我们用第1行使后3行的第1个系数为0。
对于第2行,我们直接把它加上第1行;
对于第3行,我们把它加上第1行的(-frac{1}{2})倍;
对于第4行,我们把它加上第1行的(-frac{3}{2})倍;
于是我们得到了
(left[ egin{array}{cccc} 2 & 1 & -5 & 2\0 & 0 & -4 & 3\0 & frac{11}{2} & -frac{1}{2} & 0\0 & -frac{1}{2} & frac{17}{2} & -4\ end{array} ight])(left[ egin{array}{c}12\26\14\3end{array} ight])
然后,我们用第2行使后2行的第2个系数为0。
但是我们发现第2行第2个系数本身为0。
这时我们要交换它和之后某一第2个系数不为0的行,如果不存在则方程无解/有无数组解。
因为第3行难算我们交换第2,4行,得到
(left[ egin{array}{cccc} 2 & 1 & -5 & 2\0 & -frac{1}{2} & frac{17}{2} & -4\0 & frac{11}{2} & -frac{1}{2} & 0\0 & 0 & -4 & 3\ end{array} ight])(left[ egin{array}{c}12\3\14\26end{array} ight])
然后正常工作,矩阵变为
(left[ egin{array}{cccc}2&1&-5&2\0&-frac{1}{2}&frac{17}{2}&-4\0&0&93&-44\0&0&-4&3end{array} ight])(left[ egin{array}{c}12\3\47\26end{array} ight])
数没选好的后果即将出现
然后用第3行使最后一行的第3个系数变为0。
然后为了方便计算(没方便到哪里去),交换第3,4行。
(left[ egin{array}{cccc}2&1&-5&2\0&-frac{1}{2}&frac{17}{2}&-4\0&0&-4&3\0&0&93&-44end{array} ight])(left[ egin{array}{c}12\3\26\47end{array} ight])
然后...
(left[ egin{array}{cccc}2&1&-5&2\0&-frac{1}{2}&frac{17}{2}&-4\0&0&-4&3\0&0&0&frac{147}{4}end{array} ight])(left[ egin{array}{c}12\3\26\frac{1303}{2}end{array} ight])
好在我们结束了解方程。
根据第4行,我们知道了(x_4=frac{2606}{147});
知道了(x_4),根据第3行,我们知道了(x_3);
知道了(x_3,x_4),根据第3行,我们知道了(x_2);
知道了(x_2,x_3,x_4),根据第3行,我们知道了(x_1)。
由于太难算已放弃
解完了。
模板题代码:
#include<bits/stdc++.h>
using namespace std;
struct equation{
double n[110],x;
}e[110];
int n;
double ans[110];
void Work(int x,int y){
double v=-e[y].n[x]/e[x].n[x];
for(int i=1;i<=n;i++)e[y].n[i]+=e[x].n[i]*v;
e[y].x+=e[x].x*v;
}
void Solve(){
for(int i=1;i<=n;i++){
for(int j=i;j<=n;j++){
if(e[j].n[i]){
swap(e[j],e[i]);
break;
}
}
if(!e[i].n[i])puts("No Solution"),exit(0);
for(int j=i+1;j<=n;j++)Work(i,j);
}
ans[n]=e[n].x/e[n].n[n];
double k=0;
for(int i=n;i>=1;i--,k=0){
for(int j=i+1;j<=n;j++)k+=ans[j]*e[i].n[j];
ans[i]=(e[i].x-k)/e[i].n[i];
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++)scanf("%lf",&e[i].n[j]);
scanf("%lf",&e[i].x);
}
Solve();
for(int i=1;i<=n;i++)printf("%.2lf
",ans[i]);
return 0;
}
练习题:
[JSOI2008]球形空间产生器
->Luogu
->BZOJ
->题解