zoukankan      html  css  js  c++  java
  • 高斯消元解线性方程组

    Gauss

    #include <iostream>
    #include <algorithm>
    #include <cmath>
    
    using namespace std;
    
    const double eps = 1e-6;
    const int N =110;
    
    int n;
    double a[N][N];
    void out()
    {
        for(int i = 0;i < n;i++){
            for(int j = 0;j <= n; j++)
                printf("%10.2lf ",a[i][j]);
            cout<<endl;
        }
        cout<<endl;
    }
    int gauss()
    {
        int c, r;//c表示列,r表示行
        for(c = 0, r = 0; c < n; c++)
        {
            int t = r;
            for(int i = r; i < n; i++)
                if(fabs(a[i][c]) > fabs(a[t][c])) t = i;//找到当前列绝对值最大的行
            
            if(fabs(a[t][c]) < eps) continue;
            
            for(int i = c; i <= n; i++) swap(a[t][i], a[r][i]);
            for(int i = n; i >= c; i--) a[r][i] /= a[r][c];//把该行第一个数的系数变成1
            for(int i = r + 1; i < n; i++) 
                if(fabs(a[i][c]) > eps)
                    for(int j = n; j >= c; j--)
                        a[i][j] -= a[r][j] * a[i][c];
            out();
            r++;
        }
        if(r < n)
        {
            for(int i = r; i < n; i++)
                if(fabs(a[i][n]) > eps)
                    return 2;//无解 一个未知数,两个方程 后面的a[i][n]应该全都是0的,如果不是0,矛盾。
            return 1;//无数个解 <n个方程 n个未知数(就如两个未知数,一个方程)
        }
        
        
        for(int i = n - 1; i >= 0; i--)
            for(int j = i + 1; j < n; j++)
                a[i][n] -= a[i][j] * a[j][n];
                
        
        return 0;
    }
    
    int main()
    {
        cin>>n;
        for(int i = 0;i < n;i++)
            for(int j = 0;j < n + 1;j++)
                cin>>a[i][j];
        
        int t = gauss();
        
        if(t == 0)
        {
            for(int i = 0; i < n; i++) printf("%.2lf
    ",a[i][n]);
        }
        else if(t == 1) puts("Infinite group solutions");
        else puts("No solution");
    }
                

    输入的是系数:

    3
    1.00 2.00 -1.00 -6.00
    2.00 1.00 -3.00 -9.00
    -1.00 -1.00 2.00 7.00

    输出的过程变量:1.00 0.50 -1.50 -4.50       0.00 1.50 0.50 -1.50       0.00 -0.50 0.50 2.50 

    
          1.00       0.50      -1.50      -4.50 
          0.00       1.00       0.33      -1.00 
          0.00       0.00       0.67       2.00 
    
          1.00       0.50      -1.50      -4.50 
          0.00       1.00       0.33      -1.00 
          0.00       0.00       1.00       3.00 
    
    -4.50
    -1.00
    3.00
    矩阵初等行列变换的三定律:
    1)把某一行乘一个非0的数。
    2)交换某两行
    3)把某行的若干倍加到另一行上面去。(目的是另一行消去变成0)
    高斯消元的步骤:
    1)枚举每一列column,找到当前绝对值最大的那一行。
    2)把这行换到上面去。
    3)将下面所有行的当前列column消成0。

    最后变成的应该是一个上三角。
    消去的这个步骤要对矩阵二维数组有着较为敏锐的观察:
    for(int i = r + 1; i < n; i++) 
                if(fabs(a[i][c]) > eps)
                    for(int j = n; j >= c; j--)
                        a[i][j] -= a[r][j] * a[i][c];
    比如这里代表的是:消去除r行外的后面其他所有行的当前列c为0.表示的是后面的所有数a[i][j]都减去本行第一个数a[i][c](因为第r行的c列数已经被置为1了)乘以上一行比较行a[r][j]的每一列的数字。
    for(int i = n - 1; i >= 0; i--)
            for(int j = i + 1; j < n; j++)
                a[i][n] -= a[i][j] * a[j][n];
    而这里得到的结果是每个式子的最后一个,即右边的值:a[i][n],因为左边全都只剩一个了代表各自未知数的结果。每个都是依据下一个直到最后一个式子来消元,打到左边只剩一个数字1右边的值即为结果的目的。
    a[n-1][n]即为最后一个式子未知数的解。倒数第二个未知数的值:a[n-2][n] = a[n-2][n] - a[n-2][n-1] * a[n-1][n]
                                   a[i][n] = a[i][n] - a[i][j] * a[j][n](i: n-1 -> 0, j: i + 1 ->> n - 1)(因为下一行a[j][j]已经是为1了,就是减去本行的a[i][j] * a[j][n]下一行右边的解得值)

     高斯消元解异或方程组:
    #include <iostream>
    #include <algorithm>
    using namespace std;
    
    const int N = 110;
    int n, a[N][N];
    int gauss()
    {
        int r, c;
        for(r = c = 0; c < n; c++)
        {
            int t = r;
            for(int i = r; i < n;i++)
                if(a[i][c]) {
                    t = i;
                    break;
                }
            if(!a[t][c]) continue;//如果找一圈都没有找到非0的就跳到下一列
            
            for(int i = c; i <= n;i++) swap(a[t][i], a[r][i]);//交换
            
            //用当前这一行的c列,把下面所有行的这列消成0
            for(int i = r + 1;i < n;i++)
                if(a[i][c])//已经是0的不需要消0
                    for(int j = c;j <= n; j++)
                        a[i][j] ^= a[r][j];//普通高斯消元是:a[i][j] -= a[r][j] * a[i][c];
            r++;
        }
        
        if(r < n)
        {
            for(int i = r; i < n;i++)
                if(a[i][n])//说的是如果r没有到达n的话,到了n - 2然后变为n-1就结束了,那么最后一行它就没有消因为都是0,假如还剩最后一行n-1没消,r = n - 1,但是如果右边a[n-1][n]有值,但左边为0,矛盾
                    return 2;//无解
            return 1;//无穷多组解
        }
        
        for(int i = n - 1;i >= 0;i--)
            for(int j = i + 1;j < n;j++)
                a[i][n] ^= a[i][j] & a[j][n];//普通高斯消元结果是:a[i][n] -= a[j][n] * a[i][j];
        return 0;
    }
    int main()
    {
        cin>>n;
        for(int i = 0;i < n;i++)
            for(int j = 0;j < n + 1;j++)
                cin>>a[i][j];
    
        int res = gauss();
        if(res == 0)
        {
            for(int i = 0;i<n;i++) cout<<a[i][n]<<endl;
        }
        else if(res == 1) cout<<"Multiple sets of solutions"<<endl;
        else cout<<"No solution"<<endl;
        
    }

     异或高斯校园的过程如下:

    输出

             1         1         0         1
             0         1         1         0
             0         1         0         0
    
             1         1         0         1
             0         1         1         0
             0         0         1         0
    
             1         1         0         1
             0         1         1         0
             0         0         1         0
    
    1
    0
    0
     也可以用我自己的方法,这里:

    if(a[i][j]){
    a[i][n] ^= a[j][n];
    // a[i][n] ^= a[i][j] & a[j][n];//根据模拟的结果确实是要&, 普通高斯消元结果是:a[i][n] -= a[j][n] * a[i][j];
    }

    如果a[i][j] = 1, 那么就用异或 a[i][n] ^= a[j][n]

  • 相关阅读:
    django 之 用户忘记密码的解决办法
    Django 富文本ckeditor 在模板中的实现
    MySQL密码的恢复方法
    sublime 快捷键
    linux 修改用户密码
    ubuntu 下重启 mysql
    python 控制浏览器模块
    读书笔记:从小工到专家(一)
    urlparse 模块
    python 标准内建函数
  • 原文地址:https://www.cnblogs.com/longxue1991/p/12731053.html
Copyright © 2011-2022 走看看