高斯消元求解线性方程组

    蒟蒻 Nanjo_Qi 前天考了一次试……第一题就华丽丽地爆零了。

    解一次方程组我会啊，但是解一千个有百来八十个未知数的……弃了弃了orz。

    考完了才知道有高斯消元这个神奇的东西，于是就去简单了解了一下。

    高斯消元法是线性代数规划中的一个算法，可用来为线性方程组求解，还可以求出矩阵的秩，以及求出可逆方阵的逆矩阵。消元法就是将方程组中的一方程的未知数用含有另一未知数的代数式表示，并将其代人到另一方程中，这就消去了一未知数，得到一解；或将方程组中的一方程倍乘某个常数加到另外一方程中去，也可达到消去一未知数的目的（其实就是中小学最常用的解二元，三元一次方程组的思想嘛）。我们把构成n个未知数的，含n个方程的方程组，乱搞一下，一个一个未知数求解，求出一个就再代入回去求其它的，这样就能解出所有的来啦（也有可能会无解或有多组解）！

       那该怎么搞呢？

       我们举个栗子，有这样一个三元一次方程组：

              5x + 4y + 7z = 8

              3x + 9y – 1z = -2

              8x + 3y + 2z = 9

       我们默认每一列上的未知数是相同的，如果没有就补个0，这样抛除未知数，我们就得到了一个n*(n+1)的矩阵，然后以后我们规定未知数用x1(x)，x2(y)，x3(z)等表示：

              5  4  7  8

              3  9  -1  -2

              8  3  2  9

       下面到了算法的核心：

       如果我们正在求解x_i，那么首先我们选取一个方程j>=i，使a[j][i]不为0，然后把这个方程和第i行方程互换一下。随后对于每一行j>i，我们把这个方程通过加减方程i的倍数，将其中的x_i的系数化为0，即进行消元操作。

       这个地方建议手动拿例子模拟一下，代码长这样：

　　int p = i;
　　for(int j=i; j<=n; ++j) if( fabs(a[p][i]) ) { p = j; break; }　　//选取一个方程使a[j][i]不为0
　　for(int j=i; j<=n+m; ++j) swap(a[i][j], a[p][j]);
　　for(int j=i+1; j<=n; ++i) {　　　　//向下消元
    　　if( !a[j][i] )  continue;　　　　//如果系数已经是0了，不用消元
    　　for(int k=i+1; k<=n+m; ++k)
        　　a[j][k] -= a[i][k] / a[i][i] * a[j][i];
    　　a[j][i] = 0;
　　}

       然后外面套一层循环i = 1 ~ n，完成每一个消元的操作，这时对矩阵的操作就完成啦。

       这时这个矩阵被搞成了这样，去掉常数那一边的话会是一个严格上三角矩阵（保留两位小数）：

　　在高斯消元算法中，当且仅当它是严格上三角矩阵时，才是有唯一解的。

　　　　5.00   4.00   7.00    8.00

　　　　0.00   6.60   -5.20   -6.80

　　　　0.00   0.00   -11.88  -7.30   

　　也就是：

              5x + 4y + 7z = 8

              6.6y – 5.2z = -6.8

              -11.88z = -7.3

       此时从后往前看，我们容易解出z！然后就容易解出y！然后代入解出x，就做完了。

       其实还是很简单的，这一步可以这样解决：

　　x[n] = a[n][n+i] / a[n][n];　　//未知数x[n]
　　for(int j=n-1; j>=1; --j) {　　//往前推
    　　x[j] = a[j][n+i];
    　　for(int k=j+1; k<=n; ++k)
        　　x[j] -= a[j][k] * x[k];
    　　x[j] /= a[j][j];
　　}

       这样，高斯消元法求解线性方程组的时间是O(n^3)的。

　　那无解的时候怎么办呢？

　　只需要判断是否出现某时对xi进行消元时，能否找到找到一个xi的系数不为零就可以了，就是：

　　if( fabs(a[p][i])<eps ) return 0;

　　下面是考试题目：

 

Pear 上数学课（equation.pas/c/cpp） 

题目描述：

　　Pear 被一群长得很像的一次方程组围住了，快去拯救他！

　　Pear 最近的作业之一，就是练习解一次方程组，老师给了他 m 个方程组，每个方程组未知数个数一样（n 个），未知数的系数一样，唯一不同的是每个方程的常数项。

　　比如 n=2，m=3，pear 手上有这样的 3 个方程组：

　　　　𝑥 + 𝑦 = 3     𝑥 + 𝑦 = 20     𝑥 + 𝑦 = 4

　　　　𝑥 − 𝑦 = 5     𝑥 − 𝑦 = 10     𝑥 − 𝑦 = 0

　　那他会分别解得：

　　　　𝑥 = 4      𝑥 = 15     𝑥 = 2

　　　　𝑦 = −1     𝑦 = 5      𝑦 = 2

　　现在这个任务就交给你了！

　　保证方程组有解，允许误差范围在 1e-3 以内。

输入格式：

　　第一行两个正整数 n，m；

　　之后 n 行每行 n 个整数，第 i 行的第 j 个整数 a_i,j 表示第 i 个方程中 x_i 的系数；

　　之后 m 行每行 n 个整数，第 i 行的第 j 个整数 b_i,j 表示第 m 个方程组的第 j 个方程的常数项（约定方程的左边写未知数，右边写常数项）。

输出格式：

　　输出 m 行 n 个实数，第 i 行第 j 个实数表示方程组 i 中的 x_j 的值。

输入样例：

　　2 3

　　1 1

　　1 -1

　　3 5

　　20 10

　　4 0

输出样例：

　　4 -1

　　15 5

　　2 2

数据规模：

　　对于 10%的数据 n=2；

　　对于 30%的数据 n,m<=20；

　　对于 50%的数据 n<=50,m<=100；

　　对于 100%的数据 n<=100,m<=1000。

HINT

　　本题包含 special judge。

 

　　有所不同的是题目中含有m个方程组，一个个做的时间是O(n^3*m)，会T。

　　很方便的是每一组数据所有方程组系数是一样的，只有常数不一样，而可以发现在各种对矩阵的操作中实际上常数项并不会影响什么，只是随之加加减减而已，所以不妨将这些常数全放到矩阵对应位置，我们一起操作，这样的时间是O(n^2*(n+m))的，大概是可过的范围。

　　最后贴一个代码：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const double eps = 1e-10;
const int maxn = 100 + 10;
const int maxm = 1000 + 10;
int n, m;  double a[maxn][maxn+maxm], x[maxn];

inline double Abs(double x) { return x>=0 ? x : -x; }

inline int Gaussian() {
    for(int i=1; i<=n; ++i) {
        int p = i;
        for(int j=i+1; j<=n; ++j) 
            if( Abs(a[j][i]) > Abs(a[p][i]) ) p = j;
                                //这样选可以有所避免精度问题
        if( Abs(a[p][i])<eps ) return 0;
        for(int j=i; j<=n+m; ++j) swap(a[i][j], a[p][j]);
        for(int j=i+1;j<=n;j++) {
            if( a[j][i]==0 )  continue;
            for(int k=i+1; k<=n+m; ++k)
                a[j][k] -= a[i][k] / a[i][i] * a[j][i];
            a[j][i] = 0;
        }
    }
    return 1;
}

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

    int flag = Gaussian();
    if( !flag ) { printf("-1\n"); return 0; }

    for(int i=1; i<=m; ++i) {        //对于m组方程组分别求解
        x[n] = a[n][n+i] / a[n][n];
        for(int j=n-1; j>=1; --j) {
            x[j] = a[j][n+i];
            for(int k=j+1; k<=n; ++k)
                x[j] -= a[j][k] * x[k];
            x[j] /= a[j][j];
        }
        for(int j=1; j<=n; ++j)
            printf("%.9lf ", x[j]);
        printf("\n");
    }
    return 0;
}