Description
Solution
以二维的时候为例
球面上三个点分别为((a_1,a_2))，((b_1,b_2))，((c_1,c_2))
设球心的点为((x_1,x_2))，球的半径为(T)
列式得
(left{
egin{aligned}
(a_1-x_1)^2+(a_2-x_2)^2&=(a_1)^2-2a_1x_1+(x_1)^2+(a_2)^2-2a_2x_2+(x_2)^2&=T\
(b_1-x_1)^2+(b_2-x_2)^2&=(b_1)^2-2b_1x_1+(x_1)^2+(b_2)^2-2b_2x_2+(x_2)^2&=T\
(c_1-x_1)^2+(c_2-x_2)^2&=(c_1)^2-2c_1x_1+(x_1)^2+(c_2)^2-2c_2x_2+(x_2)^2&=T
end{aligned}
ight.)
联立方程
((a_1-x_1)^2+(a_2-x_2)^2=(b_1-x_1)^2+(b_2-x_2)^2=(c_1-x_1)^2+(c_2-x_2)^2)
暴力展开
((a_1)^2-2a_1x_1+(x_1)^2+(a_2)^2-2a_2x_2+(x_2)^2=(b_1)^2-2b_1x_1+(x_1)^2+(b_2)^2-2b_2x_2+(x_2)^2=(c_1)^2-2c_1x_1+(x_1)^2+(c_2)^2-2c_2x_2+(x_2)^2)
整理得
((a_1)^2-2a_1x_1+(a_2)^2-2a_2x_2=(b_1)^2-2b_1x_1+(b_2)^2-2b_2x_2=(c_1)^2-2c_1x_1+(c_2)^2-2c_2x_2)
也可以写成这个样子
(left{
egin{aligned}
(a_1)^2-2a_1x_1+(a_2)^2-2a_2x_2&=(b_1)^2-2b_1x_1+(b_2)^2-2b_2x_2\
(b_1)^2-2b_1x_1+(b_2)^2-2b_2x_2&=(c_1)^2-2c_1x_1+(c_2)^2-2c_2x_2
end{aligned}
ight.)
移来移去
(left{
egin{aligned}
-2a_1x_1-2a_2x_2+2b_1x_1+2b_2x_2&=-(a_1)^2-(a_2)^2+(b_1)^2+(b_2)^2\
-2b_1x_1-2b_2x_2+2c_1x_1+2c_2x_2&=-(b_1)^2-(b_2)^2+(c_1)^2+(c_2)^2
end{aligned}
ight.)
移来移去
(left{
egin{aligned}
(2b_1-2a_1)x_1+(2b_2-2a_2)x_2&=-(a_1)^2-(a_2)^2+(b_1)^2+(b_2)^2\
(2c_1-2b_1)x_1+(2c_2-2b_2)x_2&=-(b_1)^2-(b_2)^2+(c_1)^2+(c_2)^2
end{aligned}
ight.)
然后我们发现它成了一个(n)元一次方程组，于是我们来高斯消元就行了
Code
#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define eps 1e-7
#define N 20
int n;
double c[N][N], a[N][N], ans[N];
inline int read() {
int s = 0, w = 1;
char c = getchar();
for (; !isdigit(c); c = getchar()) if (c == '-') w = -1;
for (; isdigit(c); c = getchar()) s = (s << 1) + (s << 3) + (c ^ 48);
return s * w;
}
int main() {
n = read();
for (register int i = 1; i <= n + 1; i++)
for (register int j = 1; j <= n; j++)
scanf("%lf", &c[i][j]);
for (register int i = 1; i <= n; i++)
for (register int j = 1; j <= n; j++)
a[i][j] = - 2 * c[i][j] + 2 * c[i + 1][j],
a[i][n + 1] += - c[i][j] * c[i][j] + c[i + 1][j] * c[i + 1][j];
for (register int i = 1; i <= n; i++) {
int r = i;
for (register int j = i + 1; j <= n; j++)
if (fabs(a[r][i]) < fabs(a[j][i]))
r = j;
if (fabs(a[r][i]) < eps) break;
if (i != r) swap(a[i], a[r]);
double div = a[i][i];
for (register int j = i; j <= n + 1; j++)
a[i][j] /= div;
for (register int j = i + 1; j <= n; j++) {
div = a[j][i];
for (register int k = i; k <= n + 1; k++)
a[j][k] -= a[i][k] * div;
}
}
ans[n] = a[n][n + 1];
for (register int i = n - 1; i; i--) {
ans[i] = a[i][n + 1];
for (register int j = i + 1; j <= n; j++)
ans[i] -= a[i][j] * ans[j];
}
for (register int i = 1; i <= n; i++)
printf("%.3lf ", ans[i]);
return 0;
}