欧氏空间 + Gauss消元法 解决方程组问题
设球心坐标为\((x_1,x_2,...,x_n)\)
对于给定的\(n+1\)个点坐标为\((a_{i1},a_{i2},...,a_{in})\)
根据题意有如下\(n+1\)个方程构成的方程组
\[\large \left\{\begin{matrix}
(a_{01}-x_1)^2 +(a_{02}-x_2)^2 +...+(a_{0n}-x_n)^2=R^2& \ ① \\
(a_{11}-x_1)^2 +(a_{12}-x_2)^2 +...+(a_{1n}-x_n)^2=R^2& \ ② \\
(a_{21}-x_1)^2 +(a_{22}-x_2)^2 +...+(a_{2n}-x_n)^2=R^2& \ ③ \\
... \\
(a_{n1}-x_1)^2 +(a_{n2}-x_2)^2 +...+(a_{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(a_{11}-a_{01})x_1+2(a_{12}-a_{02})x_2+...+2(a_{1n}-a_{0n})x_n=a_{11}^2+a_{12}^2+...+a_{1n}^2-a_{01}^2-a_{02}^2-...a_{0n}^2 \\
2(a_{21}-a_{01})x_1+2(a_{22}-a_{02})x_2+...+2(a_{2n}-a_{0n})x_n=a_{11}^2+a_{12}^2+...+a_{1n}^2-a_{01}^2-a_{02}^2-...a_{0n}^2 \\
... \\
2(a_{n1}-a_{01})x_1+2(a_{n2}-a_{02})x_2+...+2(a_{n1}-a_{0n})x_n=a_{n1}^2+a_{n2}^2+...+a_{nn}^2-a_{01}^2-a_{02}^2-...a_{0n}^2 \\
\end{matrix}\right.
\]
令:
\[\large \displaystyle
\left\{\begin{matrix}
2(a_{ij}-a_{0j})=b_{ij} \\
a_{i1}^2+a_{i2}^2+...+a_{in}^2-a_{11}^2-a_{12}^2-...a_{1n}^2=d_i
\end{matrix}\right. \ i \in [1,n],j \in [1,n]
\]
则有:
\[\large \left\{\begin{matrix}
b_{11}x_1 +b_{12}x_2+...+b_{1n}x_n=d_1 \\
b_{21}x_1 +b_{22}x_2+...+b_{2n}x_n=d_2 \\
... \\
b_{n1}x_1 +b_{n2}x_2+...+b_{nn}x_n=d_{n}
\end{matrix}\right.
\]
如此变化后,就得到一个\(n\)个变量,\(n\)个等式的线程方程组,可以使用高斯消元来处理。
实现代码
#include <bits/stdc++.h>
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++) {
// 找主元
int t = r;
for (int i = r + 1; i <= n; i++)
if (abs(b[i][c]) > abs(b[t][c])) t = i;
//交换第r行和第t行的元素
if (r != t) swap(b[t], b[r]);
//主元归一(第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[i][c] * b[r][j];
//把r++放在这里,防止记忆原始高斯消元时,找不到最大值时只跳列不跳行
r++;
}
//化成行最简阶梯型矩阵(本题唯一解,故同样也是对角矩阵)
for (int i = n; i > 1; i--)
for (int j = i - 1; j; j--)
b[j][n + 1] -= b[i][n + 1] * b[j][i];
}
int main() {
cin >> n;
for (int i = 0; i <= n; i++) // n+1行,下标从0开始
for (int j = 1; j <= n; j++) // n列
cin >> 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;
}