n维球体的通式:
(x1-a1)^2 + (x2-a2)^2 + …… + (xn-an)^2 =R ^ 2
球心坐标就是(a1,a2,……,an)
给出n+1个坐标,均带入,发现不好解,因为是二次方程,而且R^2是未知数
由于所有方程右边都是R^2 我们可以两两联立方程,消去一次 变成多元一次方程 用高斯消元搞搞就行。
例如:
(x11-a1)^2 + (x12-a2)^2 + …… + (x12-an)^2 =R ^ 2
(x21-a1)^2 + (x22-a2)^2 + …… + (x2n-an)^2 =R ^ 2
联立前两个得到:
2*(x21-x11)*a1+2*(x22-x12)+……+ 2 * (x2n-x1n) =(x21 ^ 2 - x11 ^ 2)+(x22 ^ 2 - x12 ^ 2)+……+(x2n ^ 2 - x1n ^ 2)
刚好是高斯消元初始矩阵的一行
两两联立即可
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
#define ls (o<<1)
#define rs (o<<1|1)
#define pb push_back
const double PI= acos(-1.0);
const double eps=1e-8;
const int M = 20;
double a[M][M],b[M][M];
int n;
int main()
{
scanf("%d",&n);
for(int i=1;i<=n+1;i++)
for(int j=1;j<=n;j++)scanf("%lf",&b[i][j]);
for(int i=1;i<=n;++i)
for(int j=1;j<=n;++j)
a[i][j]=2*(b[i+1][j]-b[i][j]),a[i][n+1]+=b[i+1][j]*b[i+1][j]-b[i][j]*b[i][j];
//下面高斯消元解方程
for(int i=1;i<=n;++i)//枚举列(项)
{
int max=i;
for(int j=i+1;j<=n;++j)//选出该列最大系数
if(fabs(a[j][i])>fabs(a[max][i]))max=j;
if(i!=max)swap(a[i],a[max]);//交换
if(fabs(a[i][i])<eps)//最大值等于0则说明该列都为0,肯定无解
{
puts("No Solution");
return 0;
}
for(int j=1;j<=n;++j)//每一项都减去一个数(就是小学加减消元)
{
if(j==i)continue;
double temp=a[j][i]/a[i][i];
for(int k=i+1;k<=n+1;++k)
a[j][k]-=a[i][k]*temp;
}
}
for(int i=1;i<=n;++i)
printf("%.3lf ",a[i][n+1]/a[i][i]);
return 0;
}