[高斯消元] BZOJ1013: [JSOI2008]球形空间产生器sphere

最新推荐文章于 2020-02-12 16:41:03 发布

原创最新推荐文章于 2020-02-12 16:41:03 发布 · 1.2k 阅读

0 ·

CC 4.0 BY-SA版权

102 篇文章

订阅专栏

高斯消元

4 篇文章

订阅专栏

本文介绍了一种通过解方程组来确定n维空间中n+1个点所构成球体球心坐标的算法。利用球心到各点距离相等特性，通过两方程相减的方式将原问题转化为线性方程组，并使用高斯消元法求解。

题意

给出n维空间中的n+1个点，保证这些点都在一个球体的表面。求这个n维球体的球心坐标。
给出两个定义：
球心：到球面上任意一点距离都相等的点。
距离：设两个n为空间上的点A, B的坐标为(a1, a2, …, an), (b1, b2, …, bn)，则AB的距离定义为：dist = sqrt( (a1-b1)^2 + (a2-b2)^2 + … + (an-bn)^2 )
n<=10

题解

n维,n+1个点想到什么？——解方程。
设球心坐标为( $x_1,x_2,x_3...x_n$ )， $a_{i,j}$ 为第i个点的j坐标。
根据球心到球面上任意一点距离都相等，可列出n+1个方程,其中第k个方程为：
$\sum_{i=1}^n ( x_i-a_{k,j} )^2=R^2$
注意观察式子，虽然x都是二次的，但是如果我们把两个方程相减：
$\sum_{i=1}^n ( x_i-a_{k1,j} )^2-( x_i-a_{k2,j} )^2=R^2-R^2$
$\sum_{i=1}^n (x_i-a_{k1,j}+x_i-a_{k2,j})*( x_i-a_{k1,j}-(x_i-a_{k2,j}))=R^2-R^2$
$\sum_{i=1}^n (2*x_i-a_{k1,j}-a_{k2,j})*(-a_{k1,j}+a_{k2,j})=0$
就变成一次了。
那么我们把第2~n+1个方程依次减第一个，就得到了n个方程个n未知数的线性方程。
高斯消元求解即可。

#include<cstdio>
#include<algorithm>
using namespace std;
const int maxn=25;
int n;
double a[maxn][maxn],K[maxn][maxn],ans[maxn];
int main(){
    freopen("bzoj1013.in","r",stdin);
    freopen("bzoj1013.out","w",stdout);
    scanf("%d",&n);
    for(int i=1;i<=n+1;i++)
     for(int j=1;j<=n;j++) scanf("%lf",&a[i][j]);
    for(int i=2;i<=n+1;i++)
     for(int j=1;j<=n;j++) K[i-1][j]+=(a[i][j]-a[1][j])*2, K[i-1][0]+=(a[i][j]-a[1][j])*(-a[i][j]-a[1][j]);

    for(int i=1;i<=n;i++){
        int where=0;
        for(int j=i;j<=n;j++) if (K[j][i]!=0) where=j;
        for(int j=0;j<=n;j++) swap(K[where][j],K[i][j]);
        for(int j=i+1;j<=n;j++){
            double t=K[j][i]/K[i][i];
            for(int k=0;k<=n;k++) K[j][k]-=t*K[i][k];
        }
    }
    for(int i=n;i>=1;i--){
        ans[i]=-K[i][0];
        for(int j=i+1;j<=n;j++) ans[i]-=K[i][j]*ans[j];
        ans[i]/=K[i][i];
    }
    for(int i=1;i<=n-1;i++) printf("%.3lf ",ans[i]); printf("%.3lf\n",ans[n]);
    return 0;
}