qwq
首先看到这个题,感觉就应该从列方程入手。
我们设给定的点的坐标矩阵是 x x x,然后球心坐标 a 1 , a 2 . . . . a n a_1,a_2....a_n a1,a2....an
根据欧几里得距离公式,对于一个 n 维 空 间 n维空间 n维空间的第 i i i个点,他距离球心的距离可以表示为 ∑ j = 1 n ( x i j − a [ j ] ) 2 = r 2 \sum_{j=1}^n (x_{ij}-a[j])^2 = r^2 j=1∑n(xij−a[j])2=r2
通过 n + 1 n+1 n+1个点,我们可以轻松列出来 n + 1 n+1 n+1个方程,但是可惜不是线性。我们考虑该怎么优化这个过程
考虑到相邻的两个方程的右边都是相等的,我们不妨将相邻两个方程进行减法,得到n个新的方程,假设
i
和
i
+
1
i和i+1
i和i+1进行减法
∑
j
=
1
n
(
x
i
j
−
a
[
j
]
)
2
−
∑
j
=
1
n
(
x
i
+
1
j
−
a
[
j
]
)
2
=
0
\sum_{j=1}^n (x_{ij}-a[j])^2 - \sum_{j=1}^n (x_{i+1j}-a[j])^2 = 0
j=1∑n(xij−a[j])2−j=1∑n(xi+1j−a[j])2=0
∑
j
=
1
n
x
i
j
2
−
∑
j
=
1
n
x
i
+
1
j
2
+
2
×
∑
(
x
i
+
1
j
−
x
i
j
)
∗
a
j
=
0
\sum_{j=1}^n x_{ij}^2 - \sum_{j=1}^n x_{i+1j}^2 + 2\times \sum(x_{i+1j}-x_{ij})*a_j = 0
j=1∑nxij2−j=1∑nxi+1j2+2×∑(xi+1j−xij)∗aj=0
其中
x
x
x是已知量,我们稍微移项,就可以得到最后的柿子啦
∑
(
x
i
+
1
j
−
x
i
j
)
∗
a
j
=
∑
j
=
1
n
x
i
+
1
j
2
−
∑
j
=
1
n
x
i
j
2
2
\sum(x_{i+1j}-x_{ij})*a_j = \frac{\sum_{j=1}^n x_{i+1j}^2 - \sum_{j=1}^n x_{ij}^2}{2}
∑(xi+1j−xij)∗aj=2∑j=1nxi+1j2−∑j=1nxij2
那么经过一番变化,现在我们有n个未知数,n个方程,正好可以解出来每个未知数的解。
直接上高斯消元即可
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<queue>
#include<map>
#include<set>
#define mk make_pair
#define ll long long
using namespace std;
inline int read()
{
int x=0,f=1;char ch=getchar();
while (!isdigit(ch)) {if (ch=='-') f=-1;ch=getchar();}
while (isdigit(ch)) {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return x*f;
}
const int maxn = 510;
const double eps = 1e-8;
double a[maxn][maxn];
int n,m;
double x[maxn];
double b[maxn][maxn];
double sum[maxn];
void gauss()
{
int k=1;
for (int i=1;i<=n;i++)
{
int now = k;
while (now<=n && fabs(a[now][i])<=eps) now++;
if (now==n+1) continue;
for (int j=1;j<=n+1;j++) swap(a[now][j],a[k][j]);
for (int j=1;j<=n;j++)
{
if (j!=k)
{
double t = a[j][i]/a[k][i];
for(int p=1;p<=n+1;p++) a[j][p]=a[j][p]-a[k][p]*t;
}
}
++k;
}
for (int i=n;i>=1;i--)
{
for (int j=i+1;j<=n;j++) a[i][n+1]-=a[i][j]*x[j];
x[i]=a[i][n+1]/a[i][i];
}
}
int main()
{
n=read();
for (int i=1;i<=n+1;i++)
{
for (int j=1;j<=n;j++)
scanf("%lf",&b[i][j]),sum[i]+=b[i][j]*b[i][j];
}
for (int i=1;i<=n;i++)
{
for (int j=1;j<=n;j++)
a[i][j]=b[i+1][j]-b[i][j];
a[i][n+1]=(sum[i+1]-sum[i])/2;
}
gauss();
for (int i=1;i<=n;i++) printf("%.3lf ",x[i]);
return 0;
}