分析

• 假设球心的坐标为：$(x_1, x_2, …, x_n)$，假设球的半径为R；给定的n+1个点的坐标为：$(a_{i,1}, a_{i,2}, …, a_{i,n})$。根据圆上的点到圆心的距离相等，可得到方程组：

$$\begin{cases} (a_{0,1} - x_1)^2 + (a_{0,2} - x_2)^2 + … + (a_{0,n} - x_n)^2 = R^2 \quad (1) \\\\ (a_{1,1} - x_1)^2 + (a_{1,2} - x_2)^2 + … + (a_{1,n} - x_n)^2 = R^2 \quad (2) \\\\ ...... \\\\ (a_{n,1} - x_1)^2 + (a_{n,2} - x_2)^2 + … + (a_{n,n} - x_n)^2 = R^2 \quad (n+1) \end{cases}$$

• 理论上，上述方程组中有n+1个未知数，有n+1个方程，我们可以求出这些未知数，下面考虑如何求出这些未知数。
• 考虑让上述方程组从第(2)项开始都减去第(1)项，现在考虑(2)-(1)中的第一项相减的结果：

$$(a_{1,1} - x_1)^2 - (a_{0,1} - x_1)^2 = (a_{1,1}^2 - a_{0,1}^2) - 2(a_{1,1} - a_{0,1}) \times x_1$$

所以(2)-(1)的结果如下：
$$2(a_{1,1} - a_{0,1}) \times x_1 + … + 2(a_{1,n} - a_{0,n}) \times x_n = (a_{1,1}^2 - a_{0,1}^2) + .. + (a_{1,n}^2 - a_{0,n}^2)$$
因此(i)-(1)的结果为：
$$2(a_{i,1} - a_{0,1}) \times x_1 + … + 2(a_{i,n} - a_{0,n}) \times x_n = (a_{i,1}^2 - a_{0,1}^2) + .. + (a_{i,n}^2 - a_{0,n}^2)$$

• 因为上述i可以取值范围是1~n，因此我们得到n个线性方程，一共n个未知数，可以使用高斯消元法求解。这里使用数组b替换上述系数，则有：

$$b_{i,1} \times x_1 + … + b_{i,n} \times x_n = b_{i, n+1}$$

代码

• C++

#include <iostream>
#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; c <= n; c++, r++) {
        // 找主元
        int t = r;
        for (int i = r + 1; i <= n; i++)
            if (fabs(b[i][c]) > fabs(b[t][c]))
                t = i;

        // 交换第t行和第r行
        for (int i = c; i <= n + 1; i++) swap(b[t][i], b[r][i]);
        // 归一化第r行
        for (int i = n + 1; i >= c; i--) b[r][i] /= b[r][c];
        // 消
        for (int i = c + 1; i <= n; i++)
            for (int j = n + 1; j >= c; j--)
                b[i][j] -= b[i][c] * b[r][j];
    }

    // 转换成对角矩阵(从系数矩阵最后一列开始消, 使用b[i][i]消上面的数据)
    for (int i = n; i > 1; i--)  // 枚举列
        for (int j = i - 1; j; j--) {  // 当前要把b[j][i]这个元素消成0
            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 + 1; 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();

    for (int i = 1; i <= n; i++) printf("%.3lf ", b[i][n + 1]);

    return 0;
}