AcWing 207. 球形空间产生器
\(AcWing\) \(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,...,x_n)\)
对于给定的\(n+1\)个点坐标为\((b_{i1},b_{i2},...,b_{in})\)
根据题意有如下\(n+1\)个方程构成的方程组
\[\large \left\{\begin{matrix}
(b_{01}-x_1)^2 +(b_{02}-x_2)^2 +...+(b_{0n}-x_n)^2=R^2& \ ① \\
(b_{11}-x_1)^2 +(b_{12}-x_2)^2 +...+(b_{1n}-x_n)^2=R^2& \ ② \\
(b_{21}-x_1)^2 +(b_{22}-x_2)^2 +...+(b_{2n}-x_n)^2=R^2& \ ③ \\
... \\
(b_{n1}-x_1)^2 +(b_{n2}-x_2)^2 +...+(b_{nn}-x_n)^2=R^2& \
\end{matrix}\right.
\]
因为我们学过的知识，只有 高斯消元可以解这样的方程组，但本题的方程是二次方的，肯定是不行，需要想办法把二次方消掉，化简为一次的多元线性方程组，才能使用高斯消元。
进行公式变换，把二项全部消掉：
\(\large \displaystyle \left\{\begin{matrix}
②-① \\
③-① \\
... \\
(n+1) -①
&
\end{matrix}\right.\)
有：
\[\large \displaystyle \left\{\begin{matrix}
2(b_{11}-b_{01})x_1+2(b_{12}-b_{02})x_2+...+2(b_{1n}-b_{0n})x_n=b_{11}^2+b_{12}^2+...+b_{1n}^2-b_{01}^2-b_{02}^2-...b_{0n}^2 \\
2(b_{21}-b_{01})x_1+2(b_{22}-b_{02})x_2+...+2(b_{2n}-b_{0n})x_n=b_{11}^2+b_{12}^2+...+b_{1n}^2-b_{01}^2-b_{02}^2-...b_{0n}^2 \\
... \\
2(b_{n1}-b_{01})x_1+2(b_{n2}-b_{02})x_2+...+2(b_{n1}-b_{0n})x_n=b_{n1}^2+b_{n2}^2+...+b_{nn}^2-b_{01}^2-b_{02}^2-...b_{0n}^2 \\
\end{matrix}\right.
\]
令：
\[\large \displaystyle
\left\{\begin{matrix}
2(b_{ij}-b_{0j})=a_{ij} \\
b_{i1}^2+b_{i2}^2+...+b_{in}^2-b_{01}^2-b_{02}^2-...b_{0n}^2=a_{i,{n+1}}
\end{matrix}\right. \ i,j \in [1,n]
\]
则有：
\[\large \left\{\begin{matrix}
a_{11}x_1 +a_{12}x_2+...+a_{1n}x_n=a_{1,{n+1}} \\
a_{21}x_1 +a_{22}x_2+...+a_{2n}x_n=a_{2,{n+1}} \\
... \\
a_{n1}x_1 +a_{n2}x_2+...+a_{nn}x_n=a_{n,{n+1}}
\end{matrix}\right.
\]
如此变化后，就得到一个\(n\)个变量，\(n\)个等式的 线程方程组，可以使用 高斯消元 来处理。
三、实现代码
#include <bits/stdc++.h>
using namespace std;
const int N = 15;
const double eps = 1e-8; // 小数的精度值
int n;
double a[N][N], b[N][N];
int gauss() { // 高斯消元，答案存于a[i][n]中，0 <= i < n
int r = 1; // 先按行后按列进行计算，当前行是第1行
for (int c = 1; c <= n; c++) { // 枚举每一列
int t = r; // 防破坏r,复制出t
for (int i = r; i <= n; i++) // 当前行需要找它的后续行
if (abs(a[i][c]) > abs(a[t][c])) t = i; // t的任务是找出c列中系数最大值是哪一行
if (abs(a[t][c]) < eps) continue; // 如果c列绝对值最大的系数是0, 那么处理下一列
for (int i = c; i <= n + 1; i++) swap(a[t][i], a[r][i]); // 将绝对值最大的行与当前行交换
for (int i = n + 1; i >= c; i--) a[r][i] /= a[r][c]; // a[r][c]:行首系数，将当前行的行首通过除法变为1,倒序
for (int i = r + 1; i <= n; i++) // 用当前行r的c列，通过减法将后续行c列消成0
for (int j = n + 1; j >= c; j--) // 倒序,需要保留行首，逻辑和上面是一样的，行首值是变更系数，如果正序就把系数变成1了，后面就不对了
a[i][j] -= a[r][j] * a[i][c]; // a[i][c]:需要变化的乘法系数，减法：对位相消
r++; // 下一行
}
if (r <= n) { // 如果没有成功执行完所有行，意味着中间存在continue,也就是某一列的系数都是0
for (int i = r; i <= n; i++)
if (abs(a[i][n + 1]) > eps) return 0; // 系数是0，但结果不是0，无解
return 2; // 系数是0，结果也是0，x取啥都对，有无穷多组解
}
// 代回求每个变量值
for (int i = n - 1; i; i--) // 行,倒序
for (int j = i + 1; j <= n; j++) // 列,倒三角,右上角应该都是0,对角线全是1
a[i][n + 1] -= a[i][j] * a[j][n + 1]; // 系数消为0
return 1; // 有唯一解
}
int main() {
cin >> n;
for (int i = 0; i <= n; i++) // n+1行,下标从0开始
for (int j = 1; j <= n; j++) // n列
cin >> b[i][j];
// 转换为线性方程组
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++) {
a[i][j] += 2 * (b[i][j] - b[0][j]);
a[i][n + 1] += b[i][j] * b[i][j] - b[0][j] * b[0][j];
}
// Gauss消元
gauss();
// 输出答案
for (int i = 1; i <= n; i++) printf("%.3lf ", a[i][n + 1]);
return 0;
}
浙公网安备 33010602011771号