BZOJ-1013&洛谷P4035球形空间产生器sphere【SJOI2008】高斯消元|模拟退火

Time Limit: 1 Sec  Memory Limit: 162 MB

洛谷:时间限制1.00s  内存限制125.00MB

题目链接:https://www.lydsy.com/JudgeOnline/problem.php?id=1013

洛谷:https://www.luogu.com.cn/problem/P4035

Description

  有一个球形空间产生器能够在n维空间中产生一个坚硬的球体。现在,你被困在了这个n维球体中,你只知道球
面上n+1个点的坐标,你需要以最快的速度确定这个n维球体的球心坐标,以便于摧毁这个球形空间产生器。

Input

  第一行是一个整数n(1<=N=10)。接下来的n+1行,每行有n个实数,表示球面上一点的n维坐标。每一个实数精确到小数点
后6位,且其绝对值都不超过20000。

Output

  有且只有一行,依次给出球心的n维坐标(n个实数),两个实数之间用一个空格隔开。每个实数精确到小数点
后3位。数据保证有解。你的答案必须和标准输出一模一样才能够得分。

Sample Input

2
0.0 0.0
-1.0 1.0
1.0 0.0

Sample Output

0.500 1.500

HINT

  提示:给出两个定义:1、 球心:到球面上任意一点距离都相等的点。2、 距离:设两个n为空间上的点A, B

的坐标为$(a_{1}, a_{2}, …, a_{n}), (b_{1}, b_{2}, …, b_{n})$,则AB的距离定义为:$dist=sqrt{(a_{1}-b_{1})^{2}+(a_{2}-b_{2})^{2}+cdots +(a_{n}-b_{n})^{2}}$


n维空间确实不太好想,但我们可以先考虑二维空间,二维空间中给我们圆上的两个点,那么我们可以找到此两点连线的中垂线到这两点的距离一样,所以我们还需要第三个点来确定圆心。那么我们设圆心的坐标为$(x_{1}, x_{2}, …, x_{n})$,其半径设为$r$。那么总共n+1个方程,n+1个未知数。

然后我们来列个方程:$sum_{i=1}^{n}(x_{i}-a_{i})^{2}=r^2$,有n+1个方程,解多元方程的话我们会高斯消元,但这是线性的,这个多元元2次方程不是线性的,所以我们只能想法了。

我们将相邻两个式子相减就会得到n个不含r的方程:$sum_{j=1}^{n}(o_{j}-x_{i,j})^{2}-sum_{j=1}^{n}(o_{j}-x_{i+1,j})^{2}=0$

整理可得n个不含二次项的方程:$sum_{j=1}^{n}(o_{j}-x_{i,j})^{2}-(o_{j}-x_{i+1,j})^{2}=sum_{j=1}^{n}x_{i,j}^{2}-x_{i+1,j}^{2}+2o_{j}x_{i+1,j}-2o_{j}x_{i,j}$其中o表示圆心

移项得$sum_{j=1}^{n} 2o_{j}x_{i+1,j}-2o_{j}x_{i,j}=sum_{i=1}^{n}x_{i+1,j}^{2}-x_{i,j}^{2}$,即$o_{j}$的系数为$2(x_{i+1,j}-x_{i,j})$

用线性方程就可以解答了。那么我们可以用高斯消元来解决这题,实际上高斯消元就是把矩阵化成上三角矩阵而已。然后就是一个个减下来

当然这是一种做法,还有一种比较麻烦的做法就是用模拟退火,这个做法之所以比较麻烦就是调参的时候可能会WA掉或者T掉十来发。。。这题在网上似乎也并没有完美用模拟退火A了的解答。

模拟退火的具体讲解:https://www.cnblogs.com/no-true/p/9737193.html  

https://blog.csdn.net/daaikuaichuan/article/details/81381875

随机一个点当圆心,然后和每个点比较,求出平均距离r,如果到这个点的距离大于r,说明离这个点远了,就给圆心施加一个向这个点的力;

若小于r,说明近了,就施加一个远离这个点的力。所有点比较完后,把假设的圆心按合力方向移动一个距离,距离和当前温度(T)有关。时间越久,温度越低

T要从一个比较大的数(为了开始时ans能尽快移动到圆心附近)到小于精度的数(保证答案精度),迭代次数尽量多(就是不TLE的情况下T每次乘的数尽量接近1)

为了是结果更精确,我们的距离不开方。

以下是AC代码:

#include <bits/stdc++.h>
using namespace std;

const double esp=1e-8;

double g[20][20],mat[20][20];

double pw(double a)
{
    return a*a;
}

int main()
{
    int n;
    scanf ("%d",&n);
    for (int i=1; i<=n+1; i++)
        for (int j=1; j<=n; j++)
            scanf ("%lf",&g[i][j]);
    double s;
    for (int i=1; i<=n; i++){
        s=0;
        for (int j=1; j<=n; j++){
            mat[i][j]=2*(g[i+1][j]-g[i][j]);
            s+=pw(g[i+1][j])-pw(g[i][j]);
        }
        mat[i][n+1]=s;
    }
    for (int i=1; i<=n; i++){
        int now=i;
        for (int j=i+1; j<=n; j++) 
            if (fabs(mat[j][i])>fabs(mat[now][i])) now=j;
        if (now!=i) 
            swap(mat[now],mat[i]);
        for (int k=i+1; k<=n; k++){
            double p=mat[k][i]/mat[i][i];
            for (int j=i; j<=n+1; j++)
                mat[k][j]-=mat[i][j]*p;
        }    
    }
    for (int i=n; i>=1; i--){
        for (int j=i+1; j<=n; j++)
            mat[i][n+1]-=mat[j][n+1]*mat[i][j];
        mat[i][n+1]/=mat[i][i];
    }
    for (int i=1; i<=n; i++)
        printf("%.3f%c",mat[i][n+1],i==n?'
':' ');
    return 0;
}
View Code

 模拟退火洛谷AC代码:(BZOJ上会T)

#include <bits/stdc++.h>
using namespace std;

const double inf=10000;
const double esp=1e-5;
const double lim=0.9999;

double p[20][20],ans[20],f[20],dis[20];

double dist(double *a,double *b,int n)
{
    double sum=0;
    for (int i=1; i<=n; i++)
        sum+=(a[i]-b[i])*(a[i]-b[i]);
    return sum;
}

int main()
{
    int n;
    scanf ("%d",&n);
    for (int i=1; i<=n+1; i++)
        for (int j=1; j<=n; j++)
            scanf ("%lf",&p[i][j]);
    for (int i=1; i<=n; i++) ans[i]=0;
    for (double T=inf; T>=esp; T*=lim){
        double avg=0;
        for (int i=1; i<=n+1; i++){
            dis[i]=dist(ans,p[i],n);
            avg+=dis[i];
        }
        avg/=n+1;
        for (int i=1; i<=n; i++) f[i]=0;
        for (int i=1; i<=n+1; i++)
            for (int j=1; j<=n; j++)
                f[j]+=(dis[i]-avg)*(p[i][j]-ans[j])/avg;
        for (int i=1; i<=n; i++)
            ans[i]+=f[i]*T;   
    }
    for (int i=1; i<=n; i++)
        printf("%.3f ",ans[i]);
}
View Code
路漫漫兮
原文地址:https://www.cnblogs.com/lonely-wind-/p/12217466.html