欧氏空间 + Gauss消元法 解决方程组问题
设球心坐标为 $(x_1, x_2, \cdots, x_n)^T$
对于给定的$n+1$个点坐标令为 $(a_{i1}, a_{i2}, \cdots, a_{in})^T \quad(i=1,2,\cdots,n)$
根据题意有如下 $n+1$ 个方程构成的方程组
$$
\begin{cases}
(a_{11} - x_1)^2 + (a_{12} - x_2)^2 + \cdots + (a_{1n} - x_n)^2 = R^2 \quad(1)\\\
(a_{21} - x_1)^2 + (a_{22} - x_2)^2 + \cdots + (a_{2n} - x_n)^2 = R^2 \quad(2)\\\
\cdots \\\
(a_{n+1~1} - x_1)^2 + (a_{n+1~2} - x_2)^2 + \cdots + (a_{n+1~n} - x_n)^2 = R^2 \quad(n+1)\\\
\end{cases}
$$
进行公式变换,把二次项全部消掉$\begin{cases}(2)-(1) \\\ (3)-(1) \\\ \cdots \\\ (n+1)-(1) \end{cases}$,有:
$$ \begin{cases} 2(a_{21} - a_{11})x_1 + 2(a_{22} - a_{12})x_2 + \cdots + 2(a_{2n} - a_{1n})x_n = a_{21}^2 + a_{22}^2 + \cdots + a_{2n}^2 - a_{11}^2 - a_{12}^2 - \cdots a_{1n}^2 \\\ 2(a_{31} - a_{11})x_1 + 2(a_{32} - a_{12})x_2 + \cdots + 2(a_{3n} - a_{1n})x_n = a_{31}^2 + a_{32}^2 + \cdots + a_{3n}^2 - a_{11}^2 - a_{12}^2 - \cdots a_{1n}^2 \\\ \cdots \\\ 2(a_{n+1~1} - a_{11})x_1 + 2(a_{n+1~2} - a_{12})x_2 + \cdots + 2(a_{n+1~n} - a_{1n})x_n = a_{n+1~1}^2 + a_{n+1~2}^2 + \cdots + a_{n+1~n}^2 - a_{11}^2 - a_{12}^2 - \cdots a_{1n}^2 \end{cases} $$
令$\begin{cases}2(a_{ij} - a_{1j}) = c_{ij} \\\ a_{i1}^2 + a_{i2}^2 + \cdots + a_{in}^2 - a_{11}^2 - a_{12}^2 - \cdots a_{1n}^2 = d_i\end{cases}$,则有
$$ \begin{cases} c_{21}x_1 + c_{22}x_2 + \cdots + c_{2n}x_n = d_2 \\\ c_{31}x_1 + c_{32}x_2 + \cdots + c_{3n}x_n = d_3 \\\ \cdots \\\ c_{n+1~1}x_1 + c_{n+1~2}x_2 + \cdots + c_{n+1~n}x_n = d_{n+1} \\\ \end{cases} $$
$n$个独立的线性方程 + $n$个未知数 $\Rightarrow$ 唯一解
化成 线性方程组 后,我们就可以使用 Gauss消元法 了
以下是配上详细注释的代码
#include <iostream>
#include <cstdio>
#include <cmath>
using namespace std;
const int N = 15;
int n;
double a[N][N], b[N][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;
//交换第r行和第t行的元素
for (int i = c; i <= n + 1; ++ i) swap(b[r][i], b[t][i]);
//主元归一(第r行除以主元系数)
for (int i = n + 1; i >= c; -- i) b[r][i] /= b[r][c];
//消元(用该行把下面所有行的第c列消为0)
for (int i = r + 1; i <= n; ++ i) {
for (int j = n + 1; j >= c; -- j) {
b[i][j] -= b[r][j] * b[i][c];
}
}
}
//化成行最简阶梯型矩阵(本题唯一解,故同样也是对角矩阵)
for (int i = n; i > 1; -- i) {
for (int j = i - 1; j >= 1; -- j) {
b[j][n + 1] -= b[i][n + 1] * b[j][i];
b[j][i] = 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][j] - a[0][j]);
b[i][n + 1] += a[i][j] * a[i][j] - a[0][j] * a[0][j];
}
}
//Gauss消元
Gauss();
//输出答案
for (int i = 1; i <= n; ++ i)
printf("%.3lf ", b[i][n + 1]);
return 0;
}
彩铅佬牛逼
方程组那里半径 $R$ 是不是少了个平方
是的
最讨厌高精度,其次是高斯消元(悲)