您的位置:首页 > 其它

zoj 3645 高斯消元

2013-07-22 08:38 211 查看
题意是:给12个方程,形如 (a0-x0)^2 +(a1-x1)^2+ ·········+ (a10-x10)^2=d^2;

利用嘴一个方程和上面11个方程相减,编程一次方程组 yoga高斯消元法解方程;

#include<cstdio>
#include<cmath>
#include<cstring>
#include<queue>
#include<iostream>
using namespace std;
#define FF  freopen("Input.txt","r",stdin)
#define mem(x,y) memset(x,y,sizeof(x))
const double eps=1e-12;
int n,m; //m个方程,n个变量
double Aug[15][15]; //增广矩阵
bool free_x[15]; //判断是否是不确定的变元
double x[15];    //解集
int sign(double x) { return (x>eps)- (x<-eps) ;}
/*
返回值:-1 无解,
0 有且紧一组解。
>=1 有多个解。
*/
int gauss()
{
int i,j,row,col,max_r;
for(row=0,col=0;row<m&&col<n;row++,col++)
{
max_r=row;
for(i=row+1;i<m;i++)
{
if(sign(fabs(Aug[i][col])-fabs(Aug[max_r][col]))>0)
max_r=i;
}
if(max_r!=row)
{
for(j=row;j<n+1;j++)
swap(Aug[max_r][j],Aug[row][j]);
}
if(sign(Aug[row][col])==0) { row--;continue;}
for(i=row+1;i<m;i++)
{
if(sign(Aug[i][col])==0) continue;
double ta=Aug[i][col]/Aug[row][col];
for(j=col;j<n+1;j++)
Aug[i][j]-=Aug[row][j]*ta;
}
}
for(i = row; i < m; i++) ///col=n存在0...0,a的情况,无解
{
if(sign(Aug[i][col]))
return -1;
}
if(row < n) //存在0...0,0的情况,有多个解,自由变元个数为n-row个
{
for(i = row-1; i >=0; i--)
{
int free_num = 0;   //自由变元的个数
int free_index;     //自由变元的序号
for(j = 0; j < n; j++)
{
if(sign(Aug[i][j])!=0 && free_x[j])
free_num++,free_index=j;
}
if(free_num > 1) continue; //该行中的不确定的变元的个数超过1个,无法求解,它们仍然为不确定的变元
//只有一个不确定的变元free_index,可以求解出该变元,且该变元是确定的
double tmp = Aug[i]
;
for(j = 0; j < n; j++)
{
if(sign(Aug[i][j])!=0 && j!=free_index)
tmp -= Aug[i][j]*x[j];
}
x[free_index] = tmp/Aug[i][free_index];
free_x[free_index] = false;
}
return n-row;
}//有且仅有一个解,严格的上三角矩阵(n==m)
for(i = n-1; i >= 0; i--)
{
double tmp = Aug[i]
;
for(j = i+1; j < n; j++)
if(sign(Aug[i][j])!=0)
tmp -= Aug[i][j]*x[j];
x[i] = tmp/Aug[i][i];
}
return 0;
}
int main()
{
FF;
int i,j,t;
double a[12][12];
scanf("%d",&t);
while(t--)
{
mem(Aug,0.0);
mem(x,0.0); mem(free_x,true);
for(i=0;i<12;i++)
for(j=0;j<12;j++)
scanf("%lf",&a[i][j]);
double sum=0;
for(i=0;i<11;i++) sum+=a[11][i]*a[11][i];
for(i=0;i<11;i++)
{
for(j=0;j<11;j++)
{
Aug[i][j]=2*( a[i][j]-a[11][j] );
Aug[i][11]+=a[i][j]*a[i][j];
}
Aug[i][11]+=-a[i][11]*a[i][11]+a[11][11]*a[11][11]-sum;
}
m=n=11;
gauss();
for(i=0;i<n;i++)
printf("%.2lf%c",x[i],i==n-1?'\n':' ');
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: