[浅谈] 高斯消元
\(\color{purple}\text{P3389 【模板】高斯消元法}\)
所谓高斯消元就是解个 \(n\) 元一次方程。
用矩阵记录每个方程的系数满足第 \(i\) 个方程：\(a[i][1]x_1+a[i][2]x_2+\dots +a[i][n]x_n=a[i][n+1]\)
然后从消元，一个一个项消元，如消除 \(i\) 项。先选定一个此项系数绝对值最大的方程（这样可以减小误差？），然后将他的未知数 \(i\) 的系数化为 \(1\) ，然后把剩下的方程与它相减消除未知数 \(i\) 。
最后一一带回字符串。实现的细节还得看代码。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
int read(){
int x=0,f=1;char c=getchar();
while(c>'9' || c<'0'){if(c=='-')f=-1;c=getchar();}
while(c>='0' && c<='9'){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
return x*f;
}
int n;
double mp[110][110],ans[110];
int main(){
n=read();
for(int i=1;i<=n;i++)
for(int j=1;j<=n+1;j++)
scanf("%lf",&mp[i][j]);
for(int i=1;i<=n;i++){
int rcd=i;
for(int j=i+1;j<=n;j++)
if(fabs(mp[rcd][i])<fabs(mp[j][i]))rcd=j;
if(fabs(mp[rcd][i])<eps){
printf("No Solution\n");
return 0;
}//找到系数最大的方程
if(i!=rcd)swap(mp[i],mp[rcd]);//换到当前行
double div=mp[i][i];
for(int j=i;j<=n+1;j++)
mp[i][j]/=div;//将此方程系数变为1
for(int j=i+1;j<=n;j++){
div=mp[j][i];
for(int k=i;k<=n+1;k++)
mp[j][k]-=mp[i][k]*div;//消元
}
}
ans[n]=mp[n][n+1];
for(int i=n-1;i>=1;i--){
ans[i]=mp[i][n+1];
for(int j=i+1;j<=n;j++)
ans[i]-=mp[i][j]*ans[j];
}
for(int i=1;i<=n;i++)printf("%.2lf\n",ans[i]);
return 0;
}
\(\color{purple}\text{P4035 [JSOI2008]球形空间产生器}\)
球心到所有点距离是一定的。我们把到第一个点的距离作为半径，球心到其他点的距离都等于这个半径：设球的圆心坐标用 \(b\) 表示：
\((a_1^2-2a_1b_1+b_1^2)+(a_2^2-2a_2b_2+b_2^2)=(a_{i1}^2-2a_{i1}b_1+b_1^2)+(a_{i2}^2-2_{i2}b_2+b_2^2)\)
化简：
\((2a_{i1}-2a_1)b_1+(2a_{i2}-2a_2)b_2=a_{i1}^2+a_{i2}^2-a_1^2-a_2^2\)
高斯消元即可。
点击查看代码
#include<bits/stdc++.h>
using namespace std;
const double eps=1e-8;
int read(){
int x=0,f=1;char c=getchar();
while(c>'9' || c<'0'){if(c=='-')f=-1;c=getchar();}
while(c>='0' && c<='9'){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
return x*f;
}
int n;
double mp[110][110],ans[110],tmp[50],now;
bool GSXY(){
for(int i=1;i<=n;i++){
int rcd=i;
for(int j=i+1;j<=n;j++)
if(fabs(mp[rcd][i])<fabs(mp[j][i]))rcd=j;
if(fabs(mp[rcd][i])<eps){
printf("No Solution\n");
return false;
}//找到系数最大的方程
if(i!=rcd)swap(mp[i],mp[rcd]);//换到当前行
double div=mp[i][i];
for(int j=i;j<=n+1;j++)
mp[i][j]/=div;//将此方程系数变为1
for(int j=i+1;j<=n;j++){
div=mp[j][i];
for(int k=i;k<=n+1;k++)
mp[j][k]-=mp[i][k]*div;//消元
}
}
ans[n]=mp[n][n+1];
for(int i=n-1;i>=1;i--){
ans[i]=mp[i][n+1];
for(int j=i+1;j<=n;j++)
ans[i]-=mp[i][j]*ans[j];
}
return true;
}
int main(){
n=read();
for(int i=1;i<=n;i++)scanf("%lf",&tmp[i]);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
scanf("%lf",&now);
mp[i][j]=2*now-2*tmp[j];
mp[i][n+1]+=now*now-tmp[j]*tmp[j];
}
}
GSXY();
for(int i=1;i<=n;i++)printf("%.3lf ",ans[i]);printf("\n");
return 0;
}