207. 球形空间产生器

题目链接

207. 球形空间产生器

有一个球形空间产生器能够在 $n$ 维空间中产生一个坚硬的球体。

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

注意： 数据保证有唯一解。

输入格式

第一行是一个整数 $n$。

接下来的 $n+1$ 行，每行有 $n$ 个实数，表示球面上一点的 $n$ 维坐标。

每一个实数精确到小数点后 $6$ 位，且其绝对值都不超过 $20000$。

输出格式

有且只有一行，依次给出球心的 $n$ 维坐标（$n$ 个实数），两个实数之间用一个空格隔开。

每个实数精确到小数点后 $3$ 位。

数据范围

$1≤n≤10$

输入样例：

2
0.0 0.0
-1.0 1.0
1.0 0.0

输出样例：

0.500 1.500

解题思路

高斯消元

设球心坐标为 $(x_1,x_2,\dots,x_n)$，依题意，有如下方程组（$n+1$ 个等式，$n+1$ 个未知数）：

$$\left\{\begin{array}{l} \left(a_{01}-x_{1}\right)^{2}+\left(a_{02}-x_{2}\right)^{2}+\cdots+\left(a_{0 n}-x_{n}\right)^{2}=R^{2} \\ \left(a_{11}-x_{1}\right)^{2}+\left(a_{12}-x_{2}\right)^{2}+\cdots+\left(a_{1 n}-x_{n}\right)^{2}=R^{2} \\ \cdots \\ \left(a_{n1}-x_{1}\right)^{2}+\left(a_{n2}-x_{2}\right)^{2}+\cdots+\left(a_{n n}-x_{n}\right)^{2}=R^{2} \end{array}\right.$$

相邻两个等式作差，得：

$$\left\{\begin{array}{l} 2\left(a_{0\,1}-a_{1\,1}\right) x_{1}+2\left(a_{0\,2}-a_{1\,2}\right) x_{2}+\cdots+2\left(a_{0\,n}-a_{1\, n}\right) x_{n}=(a_{0\,1}^{2}-a_{1\,1}^{2})+(a_{0\,2}^{2}-a_{1\,2}^{2}) +\dots+(a_{0\,n}^{2}-a_{1\,n}^{2})\\ 2\left(a_{1\,1}-a_{2\,1}\right) x_{1}+2\left(a_{1\,2}-a_{2\,2}\right) x_{2}+\cdots+2\left(a_{1\, n}-a_{2\, n}\right) x_{n}=(a_{1\,1}^{2}-a_{2\,1}^{2})+(a_{1\,2}^{2}-a_{2\,2}^{2}) +\dots+(a_{1\,n}^{2}-a_{2\,n}^{2})\\ \cdots \\ 2\left(a_{{n-1}\,{1}}-a_{n\,1}\right) x_{1}+2\left(a_{n-1\,2}-a_{n\,2}\right) x_{2}+\cdots+2\left(a_{n-1\, n}-a_{n\, n}\right) x_{n}=(a_{n-1\,1}^{2}-a_{n\,1}^{2})+(a_{n-1\,2}^{2}-a_{n\,2}^{2}) +\dots+(a_{n-1\,n}^{2}-a_{n\,n}^{2})\\ \end{array}\right.$$

此即为 $n$ 元一次方程组，可先将其转化为 $n$ 行 $n+1$ 列的矩阵形式，再进行高斯消元：

• 先构造成如下形式（上三角矩阵）：

$\begin{bmatrix} 1 & & & \\ 0 & 1 & & \\ 0&0 & 1 & \end{bmatrix}$

按对角线 $(r,c)$ 枚举，为了使除数不为 $0$，找绝对值最大的 $b[i][c]$，将对应的行换到 $r$ 行，将 $b[r][c]$ 化成 $1$，每次按列将 $(r,c)$ 下面列消为 $0$，

• 最后构造如下形式（简化阶梯型矩阵）：

$\begin{bmatrix} 1 &0 & 0 & \\ 0 & 1 & 0 & \\ 0 & 0& 1& \end{bmatrix}$

按列消，从最后一列往前，最终 $n+1$ 列即为答案

• 时间复杂度：$O(n^3)$

爬山法

即找到一个点，使该点到其他各点之间的距离相等，不妨先考虑一维两个点的情况：显然中点即为答案，此时中点到两个点的距离最小，引申到 $n$ 维空间，即寻找一个点使其到已知点的距离之和最小，对一维来说，$dx=(x-x0)^2$ 是凸函数，引申到 $n$ 维，由于凸函数的和仍为凸函数，所以，最终形成的函数是一个单峰，可用爬山法来求解：不用求解函数具体值，只要确定一个方向，每次往峰顶走

本题可先将中心作为备选答案，求出各个点到备选点的距离，以及各个距离的平均值，如果距离大于平均值的话，对应维度需要变化，且距离相对平均值越大的话变化就越大

• 时间复杂度：$O(随机)$

模拟退火

一般情况，模拟退火可以取代爬山法
爬山法关键在于确定方向，即通过本身及相关性质确定方向，而模拟退火关键在于跳点，即在答案范围内随机生成一个点，如果生成的点在相关性质下要比当前点的性质要强烈，且越强烈跳到该点的概率越大

本题球心要求备选点到其他各点的距离相等，可以将这些距离求出来，即要求这些距离最小，相当于方差最小，方差为 $0$ 时备选点即为答案

• 时间复杂度：$O(随机)$

代码

• 高斯消元

// Problem: 球形空间产生器
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/description/209/
// Memory Limit: 64 MB
// Time Limit: 1000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

// %%%Skyqwq
#include <bits/stdc++.h>
 
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
 
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
 
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
 
template <typename T> void inline read(T &x) {
    int f = 1; x = 0; char s = getchar();
    while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
    while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
    x *= f;
}

const int N=15;
double a[N][N],b[N][N];
int n;
void Gauss()
{
	for(int r=1,c=1;r<=n;r++,c++)
	{
		int t=r;
		for(int i=r+1;i<=n;i++)
			if(fabs(b[i][c])>fabs(b[t][c]))t=i;
		for(int i=c;i<=n+1;i++)swap(b[r][i],b[t][i]);
		for(int i=n+1;i>=c;i--)b[r][i]/=b[r][c];
		for(int i=r+1;i<=n;i++)
			for(int j=n+1;j>=c;j--)b[i][j]-=b[i][c]*b[r][j];
	}
	for(int j=n;j>1;j--)
		for(int i=j-1;i>=1;i--)
		{
			b[i][n+1]-=b[j][n+1]*b[i][j];
			b[i][j]=0;
		}
}
int main()
{
    scanf("%d",&n);
    for(int i=0;i<=n;i++)
    	for(int j=1;j<=n;j++)scanf("%lf",&a[i][j]);
    for(int i=1;i<=n;i++)
    	for(int j=1;j<=n;j++)
    	{
    		b[i][j]=2*(a[i-1][j]-a[i][j]);
    		b[i][n+1]+=a[i-1][j]*a[i-1][j]-a[i][j]*a[i][j];
    	}
	Gauss();
	for(int i=1;i<=n;i++)printf("%.3lf ",b[i][n+1]);
    return 0;
}

• 爬山法

// Problem: 球形空间产生器
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/description/209/
// Memory Limit: 64 MB
// Time Limit: 1000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

// %%%Skyqwq
#include <bits/stdc++.h>
 
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
 
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
 
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
 
template <typename T> void inline read(T &x) {
    int f = 1; x = 0; char s = getchar();
    while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
    while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
    x *= f;
}

const int N=15;
double d[N][N],dist[N],res[N],delta[N];
int n;
void cal()
{
	double avg=0;
	for(int i=0;i<=n;i++)
	{
		dist[i]=delta[i]=0;
		for(int j=1;j<=n;j++)dist[i]+=(d[i][j]-res[j])*(d[i][j]-res[j]);
		dist[i]=sqrt(dist[i]);
		avg+=dist[i]/(n+1);
	}
	for(int i=0;i<=n;i++)
		for(int j=1;j<=n;j++)delta[j]+=(d[i][j]-res[j])*(dist[i]-avg)/avg;
}
int main()
{
    scanf("%d",&n);
    for(int i=0;i<=n;i++)
    	for(int j=1;j<=n;j++)scanf("%lf",&d[i][j]),res[j]+=d[i][j]/(n+1);
    for(double T=1e4;T>1e-6;T*=0.99995)
    {
    	cal();
    	for(int i=1;i<=n;i++)res[i]+=delta[i]*T;
    }
    for(int i=1;i<=n;i++)printf("%.3lf ",res[i]);
    return 0;
}

• 模拟退火

// Problem: 球形空间产生器
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/description/209/
// Memory Limit: 64 MB
// Time Limit: 1000 ms
// 
// Powered by CP Editor (https://cpeditor.org)

// %%%Skyqwq
#include <bits/stdc++.h>
 
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
 
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
 
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
 
template <typename T> void inline read(T &x) {
    int f = 1; x = 0; char s = getchar();
    while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
    while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
    x *= f;
}

const int N=15;
double d[N][N],dist[N],res[N],cur[N],t[N],min_dist=1e9;
int n;
double rand(double l,double r)
{
	return (double)rand()/RAND_MAX*(r-l)+l;
}
double cal(double a[])
{
	double avg=0;
	for(int i=0;i<=n;i++)
	{
		dist[i]=0;
		for(int j=1;j<=n;j++)dist[i]+=(d[i][j]-a[j])*(d[i][j]-a[j]);
		avg+=dist[i]/(n+1);
	}
	double ret=0;
	for(int i=0;i<=n;i++)ret+=(dist[i]-avg)*(dist[i]-avg);
	if(min_dist>ret)
	{
		min_dist=ret;
		memcpy(res,cur,sizeof cur);
	}
	return ret;	
}

void simulate_anneal()
{
	for(double T=1e4;T>1e-6;T*=0.9999)
	{
		for(int i=1;i<=n;i++)t[i]=cur[i]+rand(-T,+T);
		double delta=cal(t)-cal(cur);
		if(exp(-delta/T)>rand(0,1))memcpy(cur,t,sizeof t);
	}
}
int main()
{
    scanf("%d",&n);
    srand(time(NULL));
    for(int i=0;i<=n;i++)
    	for(int j=1;j<=n;j++)
    	{
    		scanf("%lf",&d[i][j]);
    		cur[j]+=d[i][j]/(n+1);
    	}
    double s=clock();
    while(((double)clock()-s)/CLOCKS_PER_SEC<0.7)simulate_anneal();
    for(int i=1;i<=n;i++)printf("%.3lf ",res[i]);
    return 0;
}

__EOF__

本文作者：acwing_zyy
本文链接：https://www.cnblogs.com/zyyun/p/16145716.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！