【JSOI2008】球形空间产生器题解

最新推荐文章于 2024-07-25 21:13:04 发布

Duawieh

最新推荐文章于 2024-07-25 21:13:04 发布

阅读量313

点赞数 2

文章标签：线性代数算法 c++

本文链接：https://blog.csdn.net/qq_42097628/article/details/120450021

版权

洛谷传送门

在开始之前，先放一发目录，已经会高斯消元的朋友可以通过目录跳过第一部分。

前置知识：高斯消元法

首先我们先来了解一下高斯消元。

高斯消元，是使用矩阵的初等行变换和初等列变换求解多元一次方程组的通用算法。对于一组方程，我们可以把它们的常数（系数）部分写成矩阵的形式。

对于如下左侧方程组，我们可以将它转化为右侧的矩阵：

x + 3y + 4z = 5	1 3 4 5
x + 4y + 7z = 3	1 4 7 3
9x + 3y + 2z = 2	9 3 2 2

这样，矩阵的最后一列就表示了等号右侧的常数，前n * n的方阵则表示了各方程的系数，我们称这个方阵为系数矩阵，称加上最右侧一列的这整个矩阵为增广矩阵。

接下来我们通过一些操作在增广矩阵上完成消元。

初等行、列变换

由于是方程组转化来的矩阵，所以我们当然可以将方程组中的消元方法引入到矩阵中。小学我们就知道加减消元法和代入消元法，在增广矩阵中我们同样可以使用这两种方法。

由于这个矩阵满足的如下性质：

矩阵的任意一行同乘一个非零的数，方程组仍然成立
矩阵的任意一行同时加减其他任意一行的非零倍数，方程组仍然成立
矩阵的任意两行都可以交换位置

我们可以对矩阵的任意行进行上述操作，仍然不会影响答案的正确性，这些操作就被称为初等行变换。同样的，我们也可以在增广矩阵上定义合理的初等列变换。

高斯 - 约旦消元的思想

下面来简单介绍一下高斯消元法的流程：

选择方程组中的一个元作为“主元”
在增广矩阵中找到没有被操作过的一行，满足主元不为零（若不存在满足条件的行则跳过）
将这一行记为已经操作过的行，通过整行同时除以主元的系数，使得主元的系数化为1
通过加减消元将其他行的该主元列系数消为零
重复2~4步，直到所有主元都遍历过

上面的过程可以图解为如下形式（过程从左至右，下划线表示本行已经找出主元并进行过操作，黑体表示当前阶段进行操作的主元）：

1 3 4 5	1 3 4 5	1 3 4 5	1 3 4 5	1 0 -5 11	1 0 -5 11	1 0 0 $-\tfrac{37}{38}$	x1 = $-\tfrac{37}{38}$
1 4 7 3	1 4 7 3	0 1 3 -2	0 1 3 -2	0 1 3 -2	0 1 3 -2	0 1 0 $\tfrac{197}{38}$	x2 = $\tfrac{197}{38}$
9 3 2 2	9 3 2 2	0 -24 -34 -43	0 -24 -34 -43	0 0 38 -91	0 0 1 $-\tfrac{91}{38}$	0 0 1 $-\tfrac{91}{38}$	x3 = $-\tfrac{91}{38}$

经过上面的一系列操作，矩阵的最后一列就表示了方程组各个元的根。我们称最后得到的矩阵为对角矩阵，也叫简化阶梯矩阵。简化阶梯矩阵可以直接描述方程的根。

无解或有无数个实数解的判定

在消元的过程中，我们有可能会遇到如下情形：

1 2 -1 3	1 2 -1 3
2 4 -8 0	0 0 -6 -6
-1 -2 6 2	0 0 5 5

此时我们发现，这个增广矩阵在从左侧向右侧变化后出现了某个元找不到系数非零且未被操作过的行的情况，此时我们按照上面描述的过程将这个元跳过，继续尝试将增广矩阵化简为简化阶梯矩阵，就会得到这样的结果：

1 2 0 4	x1 + 2x2 = 4
0 0 0 0
0 0 1 1	x3 = 1

现在我们就会发现，不管x2取何值，方程都可以成立，所以我们称x2为自由元，而像x1和x3这样能够通过某式确定其值的元则称为主元。此时方程有无数多个实数解。

当然，方程之间也可能产生冲突，方程组内发生这样的矛盾会使得方程组无解。当且仅当化简过程中出现形如 0 = d (d ≠ 0) 的行时，方程组无解。上面有无数解的情况，也可以这样类似地表示为，当且仅当化简过程中出现形如 0 = 0 的行时，方程组有无数多个解。

题解：球形空间产生器

这道题的难点不在求解，而在于找出方程组。如题，题目中给出了 n+1 个点的坐标，要求在这 n+1 个点之外再求出一个点，使得其与任何其他点的距离相等。

根据题目中给出的距离计算公式:

$dist=\sqrt{(a_{1}-b_{1})_{}^{2}+(a_{2}-b_{2})_{}^{2}+\cdots+(a_{n}-b_{n})_{}^{2}}$

我们可以列出如下方程：

设点 x 为符合要求的点，其坐标为 (x1, x2, x3, ..., xn)，已知给出的 n+1 个点的坐标矩阵，矩阵第 i 行第 j 列的数字表示第 i 个点的 j 轴坐标，则对于每个 i ( $1\leqslant i\leqslant n$ ) 有：

$\large dist_{i} = dist_{i+1}$

$\sum_{j=1}^n(a_{i,j}-x_j)^2=\sum_{j=1}^n(a_{i+1,j}-x_j)^2$

对其进行化简，可得：

$\sum_{j=1}^n(a_{i,j}^2+x_j^2-2a_{i,j}x_j)=\sum_{j=1}^n(a_{i+1,j}^2+x_{j}^2-2a_{i+1,j}x_j)$

$\sum_{j=1}^n(a_{i,j}^2-a_{i+1,j}^2-2a_{i,j}x_j+2a_{i+1,j}x_j)=0$

$\sum_{j=1}^{n}[(a_{i,j}+a_{i+1,j})(a_{i,j}-a_{i+1,j})-2x_{j}(a_{i,j}-a_{i+1,j})]=0$

将其移项，我们就得到了关于 x1, x2, ..., xn 的 n元一次方程组：

$\sum_{j=1}^n2x_j(a_{i,j}-a_{i+1,j})=\sum_{j=1}^n(a_{i,j}+a_{i+1,j})(a_{i,j}-a_{i+1,j})$ $(1\leq i\leq n)$

列出方程，将其填入增广矩阵，我们就可以进行消元求解了，由于本题数据保证有解，所以代码中不需要对根的数量进行特殊判定。

#include<bits/stdc++.h>
using namespace std;
//球形空间产生器
const int maxn = 12;
double p[maxn][maxn];    //输入的点坐标矩阵
double a[maxn][maxn];    //计算得出的方程组增广矩阵
int n;                   //方程阶数

inline bool zr(double x) {    //判断一个浮点数是否非零
    if (abs(x) >= 0.0000001) return true;
    else return false;
}

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n + 1; i++)
        for (int j = 1; j <= n; j++)
            scanf("%lf", &p[i][j]);

    for (int i = 1; i <= n; i++) {
        for (int j = 1; j <= n; j++) {                        //求出一个 n*n+1 的增广矩阵
            a[i][j] = 2 * (p[i][j] - p[i + 1][j]);            //求出系数矩阵
            a[i][n + 1] += (p[i][j] - p[i + 1][j]) * (p[i][j] + p[i + 1][j]);//求出等号右侧常量
        }
    }



    for (int j = 1; j <= n; j++) {        //枚举每个元作为主元进行操作
        for (int i = j; i <= n; i++) {
            if (zr(a[i][j])) {            //在未被操作过的行中找出一行，满足主元的系数非零
                double rate = a[i][j];
                for (int k = 1; k <= n + 1; k++) {
                    a[i][k] /= rate;      //将操作行的主元系数化为“1”
                    if (i != j) swap(a[i][k], a[j][k]);    //交换第j行和操作行
                }
                break;
            }
        }

        for (int i = 1; i <= n; i++) {    //将其他行中主元的系数消为“0”
            if (i == j) continue;
            double rate = a[i][j] / a[j][j];
            for (int k = 1; k <= n + 1; k++) a[i][k] -= a[j][k] * rate;
        }
    }

    for (int i = 1; i <= n; i++) printf("%.3lf ", a[i][n + 1]);    //输出答案
    return 0;
}

Duawieh

关注

2
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
【JSOI2008】球形空间产生器题解

洛谷传送门在开始之前，先放一发目录，已经会高斯消元的朋友可以通过目录跳过第一部分。目录前置知识：高斯消元法初等行、列变换高斯 - 约旦消元的思想无解或有无数个实数解的判定题解：球形空间产生器前置知识：高斯消元法首先我们先来了解一下高斯消元。高斯消元，是使用矩阵的初等行变换和初等列变换求解多元一次方程组的通用算法。对于一组方程，我们可以把它们的常数部分写成矩阵的形式。对于如下左侧方程组，我们可以将它转化为右侧的矩阵...
复制链接

扫一扫