高斯消元 - JSOI 2008 球形空间产生器 - 洛谷 P4035
有一个球形空间产生器能够在n维空间中产生一个坚硬的球体。
现在,你被困在了这个n维球体中,你只知道球面上n+1个点的坐标,你需要以最快的速度确定这个n维球体的球心坐标,以便于摧毁这个球形空间产生器。
注意: 数据保证有唯一解。
输入格式
第一行是一个整数n。
接下来的n+1行,每行有n个实数,表示球面上一点的n维坐标。
每一个实数精确到小数点后6位,且其绝对值都不超过20000。
输出格式
有且只有一行,依次给出球心的n维坐标(n个实数),两个实数之间用一个空格隔开。
每个实数精确到小数点后3位。
数据范围
1 ≤ n ≤ 10 1≤n≤10 1≤n≤10
输入样例:
2
0.0 0.0
-1.0 1.0
1.0 0.0
输出样例:
0.500 1.500
分析:
设 n 维 坐 标 下 的 球 心 为 ( x 1 , x 2 , . . . , x n ) , n + 1 个 点 的 坐 标 为 ( a i 1 , a i 2 , . . . , a i n ) 1 ≤ i ≤ n + 1 , 球 半 径 为 R , 则 : 设n维坐标下的球心为(x_1,x_2,...,x_n),n+1个点的坐标为(a_{i1},a_{i2},...,a_{in})_{1≤i≤n+1},球半径为R,则: 设n维坐标下的球心为(x1,x2,...,xn),n+1个点的坐标为(ai1,ai2,...,ain)1≤i≤n+1,球半径为R,则:
{ ( a 11 − x 1 ) 2 + ( a 12 − x 2 ) 2 + . . . + ( a 1 n − x n ) 2 = R 2 ( a 21 − x 1 ) 2 + ( a 22 − x 2 ) 2 + . . . + ( a 2 n − x n ) 2 = R 2 . . . ( a n + 1 , 1 − x 1 ) 2 + ( a n + 1 , 2 − x 2 ) 2 + . . . + ( a n + 1 , n − x n ) 2 = R 2 \begin{cases}(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_{n+1,1}-x_1)^2+(a_{n+1,2}-x_2)^2+...+(a_{n+1,n}-x_n)^2=R^2\end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧(a11−x1)2+(a12−x2)2+...+(a1n−xn)2=R2(a21−x1)2+(a22−x2)2+...+(a2n−xn)2=R2...(an+1,1−x1)2+(an+1,2−x2)2+...+(an+1,n−xn)2=R2
本 题 看 似 是 二 次 方 程 组 , 但 是 我 们 发 现 , 每 个 方 程 的 二 次 项 x i 2 都 对 应 相 等 , 本题看似是二次方程组,但是我们发现,每个方程的二次项x_i^2都对应相等, 本题看似是二次方程组,但是我们发现,每个方程的二次项xi2都对应相等,
故 我 们 对 任 意 两 个 方 程 作 差 , 就 能 够 得 到 关 于 x i 的 线 性 方 程 。 故我们对任意两个方程作差,就能够得到关于x_i的线性方程。 故我们对任意两个方程作差,就能够得到关于xi的线性方程。
本 题 , 将 所 有 方 程 减 去 第 一 个 方 程 , 得 到 新 的 线 性 方 程 组 : 本题,将所有方程减去第一个方程,得到新的线性方程组: 本题,将所有方程减去第一个方程,得到新的线性方程组:
{ a 21 2 − a 11 2 − 2 ( a 21 − a 11 ) x 1 + a 22 2 − a 12 2 − 2 ( a 22 − a 12 ) x 2 + . . . + a 2 n 2 − a 1 n 2 − 2 ( a 2 n − a 1 n ) x n = 0 a 31 2 − a 11 2 − 2 ( a 31 − a 11 ) x 1 + a 32 2 − a 12 2 − 2 ( a 32 − a 12 ) x 2 + . . . + a 3 n 2 − a 1 n 2 − 2 ( a 3 n − a 1 n ) x n = 0 . . . a n + 1 , 1 2 − a 11 2 − 2 ( a n + 1 , 1 − a 11 ) x 1 + a n + 1 , 2 2 − a 12 2 − 2 ( a n + 1 , 2 − a 12 ) x 2 + . . . + a n + 1 , n 2 − a 1 n 2 − 2 ( a n + 1 , n − a 1 n ) x n = 0 \begin{cases}a_{21}^2-a_{11}^2-2(a_{21}-a_{11})x_1+a_{22}^2-a_{12}^2-2(a_{22}-a_{12})x_2+...+a_{2n}^2-a_{1n}^2-2(a_{2n}-a_{1n})x_n=0 \\\\a_{31}^2-a_{11}^2-2(a_{31}-a_{11})x_1+a_{32}^2-a_{12}^2-2(a_{32}-a_{12})x_2+...+a_{3n}^2-a_{1n}^2-2(a_{3n}-a_{1n})x_n=0 \\\\...\\\\ a_{n+1,1}^2-a_{11}^2-2(a_{n+1,1}-a_{11})x_1+a_{n+1,2}^2-a_{12}^2-2(a_{n+1,2}-a_{12})x_2+...+a_{n+1,n}^2-a_{1n}^2-2(a_{n+1,n}-a_{1n})x_n=0 \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a212−a112−2(a21−a11)x1+a222−a122−2(a22−a12)x2+...+a2n2−a1n2−2(a2n−a1n)xn=0a312−a112−2(a31−a11)x1+a322−a122−2(a32−a12)x2+...+a3n2−a1n2−2(a3n−a1n)xn=0...an+1,12−a112−2(an+1,1−a11)x1+an+1,22−a122−2(an+1,2−a12)x2+...+an+1,n2−a1n2−2(an+1,n−a1n)xn=0
等 价 于 : 等价于: 等价于:
{ 2 ( a 21 − a 11 ) x 1 + 2 ( a 22 − a 12 ) x 2 + . . . + 2 ( a 2 n − a 1 n ) x n = a 21 2 − a 11 2 + a 22 2 − a 12 2 + . . . + a 2 n 2 − a 1 n 2 2 ( a 31 − a 11 ) x 1 + 2 ( a 32 − a 12 ) x 2 + . . . + 2 ( a 3 n − a 1 n ) x n = a 31 2 − a 11 2 + a 32 2 − a 12 2 + . . . + a 3 n 2 − a 1 n 2 . . . 2 ( a n + 1 , 1 − a 11 ) x 1 + 2 ( a n + 1 , 2 − a 12 ) x 2 + . . . + 2 ( a n + 1 , n − a 1 n ) x n = a n + 1 , 1 2 − a 11 2 + a n + 1 , 2 2 − a 12 2 + . . . + a n + 1 , n 2 − a 1 n 2 \begin{cases}2(a_{21}-a_{11})x_1+2(a_{22}-a_{12})x_2+...+2(a_{2n}-a_{1n})x_n=a_{21}^2-a_{11}^2+a_{22}^2-a_{12}^2+...+a_{2n}^2-a_{1n}^2 \\\\2(a_{31}-a_{11})x_1+2(a_{32}-a_{12})x_2+...+2(a_{3n}-a_{1n})x_n=a_{31}^2-a_{11}^2+a_{32}^2-a_{12}^2+...+a_{3n}^2-a_{1n}^2 \\\\...\\\\ 2(a_{n+1,1}-a_{11})x_1+2(a_{n+1,2}-a_{12})x_2+...+2(a_{n+1,n}-a_{1n})x_n=a_{n+1,1}^2-a_{11}^2+a_{n+1,2}^2-a_{12}^2+...+a_{n+1,n}^2-a_{1n}^2 \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧2(a21−a11)x1+2(a22−a12)x2+...+2(a2n−a1n)xn=a212−a112+a222−a122+...+a2n2−a1n22(a31−a11)x1+2(a32−a12)x2+...+2(a3n−a1n)xn=a312−a112+a322−a122+...+a3n2−a1n2...2(an+1,1−a11)x1+2(an+1,2−a12)x2+...+2(an+1,n−a1n)xn=an+1,12−a112+an+1,22−a122+...+an+1,n2−a1n2
增广矩阵:
[
2
(
a
21
−
a
11
)
2
(
a
22
−
a
12
)
.
.
.
2
(
a
2
n
−
a
1
n
)
a
21
2
−
a
11
2
+
a
22
2
−
a
12
2
+
.
.
.
+
a
2
n
2
−
a
1
n
2
2
(
a
31
−
a
11
)
2
(
a
32
−
a
12
)
.
.
.
2
(
a
3
n
−
a
1
n
)
a
31
2
−
a
11
2
+
a
32
2
−
a
12
2
+
.
.
.
+
a
3
n
2
−
a
1
n
2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2
(
a
n
+
1
,
1
−
a
11
)
2
(
a
n
+
1
,
2
−
a
12
)
.
.
.
2
(
a
n
+
1
,
n
−
a
1
n
)
a
n
+
1
,
1
2
−
a
11
2
+
a
n
+
1
,
2
2
−
a
12
2
+
.
.
.
+
a
n
+
1
,
n
2
−
a
1
n
2
]
\left[ \begin{array}{cccc|c} 2(a_{21}-a_{11})&2(a_{22}-a_{12})&...&2(a_{2n}-a_{1n})&a_{21}^2-a_{11}^2+a_{22}^2-a_{12}^2+...+a_{2n}^2-a_{1n}^2\\\\ 2(a_{31}-a_{11})&2(a_{32}-a_{12})&...&2(a_{3n}-a_{1n})&a_{31}^2-a_{11}^2+a_{32}^2-a_{12}^2+...+a_{3n}^2-a_{1n}^2\\\\ ...&...&...&...&...\\\\ 2(a_{n+1,1}-a_{11})&2(a_{n+1,2}-a_{12})&...&2(a_{n+1,n}-a_{1n})&a_{n+1,1}^2-a_{11}^2+a_{n+1,2}^2-a_{12}^2+...+a_{n+1,n}^2-a_{1n}^2 \end{array} \right]
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡2(a21−a11)2(a31−a11)...2(an+1,1−a11)2(a22−a12)2(a32−a12)...2(an+1,2−a12)............2(a2n−a1n)2(a3n−a1n)...2(an+1,n−a1n)a212−a112+a222−a122+...+a2n2−a1n2a312−a112+a322−a122+...+a3n2−a1n2...an+1,12−a112+an+1,22−a122+...+an+1,n2−a1n2⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
令 b i j = a i + 1 , j − a 1 j , 则 b i , n + 1 = ∑ j = 1 n a i + 1 , j 2 − a 1 j 2 , 用 高 斯 消 元 解 方 程 组 即 可 。 令b_{ij}=a_{i+1,j}-a_{1j},则b_{i,n+1}=\sum_{j=1}^{n}a_{i+1,j}^2-a_{1j}^2,用高斯消元解方程组即可。 令bij=ai+1,j−a1j,则bi,n+1=∑j=1nai+1,j2−a1j2,用高斯消元解方程组即可。
代码:
#include<cstdio>
#include<cstring>
#include<algorithm>
#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;r++,c++)
{
//找主元
int t=r;
for(int i=r+1;i<=n;i++)
if(fabs(b[i][c])>fabs(b[t][c]))
t=i;
//交换
for(int i=c;i<=n+1;i++) swap(b[r][i],b[t][i]);
//归一化
for(int i=n+1;i>=c;i--) b[r][i]/=b[r][c];
//消成上三角阵
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];
}
//转化成对角阵
for(int i=n;i>1;i--) //i为列,j为行
for(int j=i-1;j;j--)
{
b[j][n+1]-=b[i][n+1]*b[j][i];
b[j][i]=0;
}
}
int main()
{
scanf("%d",&n);
for(int i=1;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+1][j]-a[1][j]);
b[i][n+1]+=a[i+1][j]*a[i+1][j]-a[1][j]*a[1][j];
}
gauss();
for(int i=1;i<=n;i++) printf("%.3lf ",b[i][n+1]);
return 0;
}