高斯消元

高斯消元

1. 高斯消元原理

原理

  • 高斯消元是用于求解线性方程组的方法。给定如下方程组:

{ a 11 × x 1 + a 12 × x 2 + . . . . . . + a 1 n × x n = b 1 a 21 × x 1 + a 22 × x 2 + . . . . . . + a 2 n × x n = b 2 . . . . . . a n 1 × x 1 + a n 2 × x 2 + . . . . . . + a n n × x n = b n \begin{cases} a_{11} \times x_1 + a_{12} \times x_2 + ...... + a_{1n} \times x_n = b_1 \\ a_{21} \times x_1 + a_{22} \times x_2 + ...... + a_{2n} \times x_n = b_2 \\ ...... \\ a_{n1} \times x_1 + a_{n2} \times x_2 + ...... + a_{nn} \times x_n = b_n \\ \end{cases} a11×x1+a12×x2+......+a1n×xn=b1a21×x1+a22×x2+......+a2n×xn=b2......an1×x1+an2×x2+......+ann×xn=bn

  • 我们是通过初等变换对上述方程组进行求解的,初等变换存在三种:

    (1)把某一行乘以一个非零的数;

    (2)交换某两行;

    (3)把某行的若干倍加到另一行上去;

  • 总的来说高斯消元分为两大步骤:第一步将系数矩阵通过初等变化变为上三角矩阵;回代求解。

  • 第一步:转为上三角矩阵:枚举每一列c

    (1)找到绝对值最大的一行;

    (2)将该行换到最上面;

    (3)将该行第1个数变为1;

    (4)将下面所有行的第c列消成0;

    类似于下图,但是缺少(2)(3)步。

在这里插入图片描述

  • 第二步:回代求解。

    从最后一个方程进行回代可以求出所有的答案。

  • 结果存在三种情况(做完第一步之后):

    (1)唯一解:满秩,即系数矩阵没有全0行;

    (2)无解:系数矩阵存在全0行,且该行对应的b值不为0;

    (3)无穷多组解:系数矩阵存在全0行,且0行对应的b值都为0;

2. AcWing上的高斯消元题目

AcWing 883. 高斯消元解线性方程组

问题描述

分析

  • 根据上述分析实现代码即可。

代码

  • C++
#include <iostream>
#include <cmath>

using namespace std;

const int N = 110;
const double eps = 1e-6;  // 浮点数小于eps认为是0

int n;
double a[N][N];  // 存储增广矩阵

void print() {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n + 1; j++)
            printf("%10.2lf", a[i][j]);
        puts("");
    }
    puts("");
}

// 返回0:存在唯一解; 返回1:存在无穷多组解; 返回2:无解
int gauss() {
    
    int r, c;  // 行、列
    for (r = 0, c = 0; c < n; c++) {
        // 第一步:(1) 找到绝对值最大的一行
        int t = r;  // 当前考察的为第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;  // 说明后面的行都为0
        
        // 第一步:(2) 将该行换到最上面
        for (int i = c; i < n + 1; i++) swap(a[t][i], a[r][i]);
        // 第一步:(3) 将第r行第c列对应的数变为1
        for (int i = n; i >= c; i--) a[r][i] /= a[r][c];
        // 第一步:(4) 将下面所有行的第c列消成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++;
    }
    
    if (r < n) {  // 说明系数矩阵存在0行
        for (int i = r; i < n; i++)
            if (fabs(a[i][n]) > eps)
                return 2;
        return 1;
    }
    
    // 第二步:回代
    for (int i = n - 1; i >= 0; i--)
        for (int j = i + 1; j < n; j++)  // j+1行到最后一行
            a[i][n] -= a[j][n] * a[i][j];
    
    return 0;
}

int main() {
    
    scanf("%d", &n);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n + 1; j++)
            scanf("%lf", &a[i][j]);
    
    int t = gauss();
    
    if (t == 0) {
        for (int i = 0; i < n; i++) printf("%.2lf\n", a[i][n]);
    } else if (t == 1) puts("Infinite group solutions");
    else puts("No solution");
    
    return 0;
}

AcWing 884. 高斯消元解异或线性方程组

问题描述

分析

  • 异或又可以被看成不进位的加法。我们将常规高斯消元中的加减法换成异或运算,乘法换成与运算即可。

  • 这里的步骤是和高斯消元法一样的,只不过需要将运算换成异或运算。

  • 第一步:转为上三角矩阵:枚举每一列c

    (1)找到非零行;

    (2)将该行换到最上面;

    (3)将下面所有行的第c列消成0;

  • 第二步:回代求解。

    从最后一个方程进行回代可以求出所有的答案。

  • 结果存在三种情况:

    (1)唯一解:做完第一步之后,满秩,即系数矩阵没有全0行;

    (2)无解:系数矩阵存在全0行,且该行对应的b值不为0;

    (3)无穷多组解:系数矩阵存在全0行,且0行对应的b值都为0;

代码

  • C++
#include <iostream>

using namespace std;

const int N = 110;

int n;
int a[N][N];

// 返回0:存在唯一解; 返回1:存在无穷多组解; 返回2:无解
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;
        
        for (int i = c; i < n + 1; i++) swap(a[t][i], a[r][i]);
        
        for (int i = r + 1; i < n; i++)
            if (a[i][c]) 
                for (int j = c; j < n + 1; j++)
                    a[i][j] ^= a[r][j];
        
        r++;
    }
    
    if (r < n) {
        for (int i = r; i < n; i++)
            if (a[i][n])
                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];
    
    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("%d\n", a[i][n]);
    } else if (t == 1) {
        puts("Multiple sets of solutions");
    } else puts("No solution");
    
    return 0;
}

AcWing 207. 球形空间产生器

问题描述

分析

  • 假设球心的坐标为: ( x 1 , x 2 , . . . , x n ) (x_1, x_2, ..., x_n) (x1,x2,...,xn),假设球的半径为R;给定的n+1个点的坐标为: ( a i , 1 , a i , 2 , . . . , a i , n ) (a_{i,1}, a_{i,2}, ..., a_{i,n}) (ai,1,ai,2,...,ai,n)。根据圆上的点到圆心的距离相等,可得到方程组:

{ ( a 0 , 1 − x 1 ) 2 + ( a 0 , 2 − x 2 ) 2 + . . . + ( a 0 , n − x n ) 2 = R 2 ( 1 ) ( a 1 , 1 − x 1 ) 2 + ( a 1 , 2 − x 2 ) 2 + . . . + ( a 1 , n − x n ) 2 = R 2 ( 2 ) . . . . . . ( a n , 1 − x 1 ) 2 + ( a n , 2 − x 2 ) 2 + . . . + ( a n , n − x n ) 2 = R 2 ( n + 1 ) \begin{cases} (a_{0,1} - x_1)^2 + (a_{0,2} - x_2)^2 + ... + (a_{0,n} - x_n)^2 = R^2 \quad (1) \\ (a_{1,1} - x_1)^2 + (a_{1,2} - x_2)^2 + ... + (a_{1,n} - x_n)^2 = R^2 \quad (2) \\ ...... \\ (a_{n,1} - x_1)^2 + (a_{n,2} - x_2)^2 + ... + (a_{n,n} - x_n)^2 = R^2 \quad (n+1) \end{cases} (a0,1x1)2+(a0,2x2)2+...+(a0,nxn)2=R2(1)(a1,1x1)2+(a1,2x2)2+...+(a1,nxn)2=R2(2)......(an,1x1)2+(an,2x2)2+...+(an,nxn)2=R2(n+1)

  • 理论上,上述方程组中有n+1个未知数,有n+1个方程,我们可以求出这些未知数,下面考虑如何求出这些未知数。

  • 考虑让上述方程组从第(2)项开始都减去第(1)项,现在考虑(2)-(1)中的第一项相减的结果:

( a 1 , 1 − x 1 ) 2 − ( a 0 , 1 − x 1 ) 2 = ( a 1 , 1 2 − a 0 , 1 2 ) − 2 ( a 1 , 1 − a 0 , 1 ) × x 1 (a_{1,1} - x_1)^2 - (a_{0,1} - x_1)^2 = (a_{1,1}^2 - a_{0,1}^2) - 2(a_{1,1} - a_{0,1}) \times x_1 (a1,1x1)2(a0,1x1)2=(a1,12a0,12)2(a1,1a0,1)×x1

所以(2)-(1)的结果如下:
2 ( a 1 , 1 − a 0 , 1 ) × x 1 + . . . + 2 ( a 1 , n − a 0 , n ) × x n = ( a 1 , 1 2 − a 0 , 1 2 ) + . . + ( a 1 , n 2 − a 0 , n 2 ) 2(a_{1,1} - a_{0,1}) \times x_1 + ... + 2(a_{1,n} - a_{0,n}) \times x_n = (a_{1,1}^2 - a_{0,1}^2) + .. + (a_{1,n}^2 - a_{0,n}^2) 2(a1,1a0,1)×x1+...+2(a1,na0,n)×xn=(a1,12a0,12)+..+(a1,n2a0,n2)
因此(i)-(1)的结果为:
2 ( a i , 1 − a 0 , 1 ) × x 1 + . . . + 2 ( a i , n − a 0 , n ) × x n = ( a i , 1 2 − a 0 , 1 2 ) + . . + ( a i , n 2 − a 0 , n 2 ) 2(a_{i,1} - a_{0,1}) \times x_1 + ... + 2(a_{i,n} - a_{0,n}) \times x_n = (a_{i,1}^2 - a_{0,1}^2) + .. + (a_{i,n}^2 - a_{0,n}^2) 2(ai,1a0,1)×x1+...+2(ai,na0,n)×xn=(ai,12a0,12)+..+(ai,n2a0,n2)

  • 因为上述i可以取值范围是1~n,因此我们得到n个线性方程,一共n个未知数,可以使用高斯消元法求解。这里使用数组b替换上述系数,则有:

b i , 1 × x 1 + . . . + b i , n × x n = b i , n + 1 b_{i,1} \times x_1 + ... + b_{i,n} \times x_n = b_{i, n+1} bi,1×x1+...+bi,n×xn=bi,n+1

代码

  • C++
#include <iostream>
#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; c++, r++) {
        // 找主元
        int t = r;
        for (int i = r + 1; i <= n; i++)
            if (fabs(b[i][c]) > fabs(b[t][c]))
                t = i;

        // 交换第t行和第r行
        for (int i = c; i <= n + 1; i++) swap(b[t][i], b[r][i]);
        // 归一化第r行
        for (int i = n + 1; i >= c; i--) b[r][i] /= b[r][c];
        // 消
        for (int i = c + 1; i <= n; i++)
            for (int j = n + 1; j >= c; j--)
                b[i][j] -= b[i][c] * b[r][j];
    }
    
    // 转换成对角矩阵(从系数矩阵最后一列开始消, 使用b[i][i]消上面的数据)
    for (int i = n; i > 1; i--)  // 枚举列
        for (int j = i - 1; j; j--) {  // 当前要把b[j][i]这个元素消成0
            b[j][n + 1] -= b[i][n + 1] * b[j][i];
            b[j][i] = 0;
        }
}

int main() {
    
    scanf("%d", &n);
    for (int i = 0; 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][j] - a[0][j]);
            b[i][n + 1] += a[i][j] * a[i][j] - a[0][j] * a[0][j];
        }
    
    gauss();
    
    for (int i = 1; i <= n; i++) printf("%.3lf ", b[i][n + 1]);
    
    return 0;
}

AcWing 208. 开关问题

问题描述

分析

  • 这里以一个例子讲解这个题目,有三个开关,分别为1、2、3,规则是:如果按1,则1、2的状态会改变;如果按2,则1、2、3的状态会改变;如果按3,则3的状态会改变。初始状态是0、0、0,转变为1、1、1存在多少种方案?

  • 我们使用 x 1 、 x 2 、 x 3 x_1、x_2、x_3 x1x2x3 分别表示每个开关有没有发生变化,没变为0,变了为1。根据上述规则,按1、2会影响1的状态,因此x1 ^ x2 = 1,同理我们可以得到如下异或方程组:

{ x 1 ⊕ x 2 = 1 第 一 个 开 关 x 1 ⊕ x 2 = 1 第 二 个 开 关 x 2 ⊕ x 3 = 1 第 三 个 开 关 \begin{cases} x_1 \oplus x_2 = 1 \quad 第一个开关 \\ x_1 \oplus x_2 = 1 \quad 第二个开关 \\ x_2 \oplus x_3 = 1 \quad 第三个开关 \end{cases} x1x2=1x1x2=1x2x3=1

  • 因此此时我们就将问题转化成为了AcWing 884. 高斯消元解异或线性方程组。对增广矩阵进行变换,变为上三角矩阵,最后如果有k个自由元(即k个行全0),最终的结果就为 2 k 2^k 2k

  • 上述例题存在一个自由元,因此答案是2,两个答案分别是 ( x 1 = 0 , x 2 = 1 , x 3 = 0 ) 、 ( x 1 = 1 , x 2 = 0 , x 3 = 1 ) (x_1=0, x_2=1, x_3=0)、(x_1=1, x_2=0, x_3=1) (x1=0,x2=1,x3=0)(x1=1,x2=0,x3=1)

代码

  • C++
#include <iostream>
#include <cstring>

using namespace std;

const int N = 35;

int n;
int a[N][N];

int gauss() {
    
    int r, c;
    for (r = 1, c = 1; c <= n; c++) {
        // 找主元
        int t = r;
        for (int i = r; i <= n; i++)
            if (a[i][c])
                t = i;
        if (a[t][c] == 0) continue;
        // 交换
        for (int i = c; i <= n + 1; i++) swap(a[r][i], a[t][i]);
        // 消
        for (int i = r + 1; i <= n; i++)
            for (int j = n + 1; j >= c; j--)
                a[i][j] ^= a[i][c] & a[r][j];
        r++;
    }
    
    int res = 1;
    if (r < n + 1) {
        for (int i = r; i <= n; i++) {
            if (a[i][n + 1]) return -1;  // 出现了 0 == !0,无解
            res *= 2;
        }
    }
    return res;
}

int main() {
    
    int T;
    scanf("%d", &T);
    while (T--) {
        scanf("%d", &n);
        memset(a, 0, sizeof a);
        for (int i = 1; i <= n; i++) scanf("%d", &a[i][n + 1]);
        for (int i = 1; i <= n; i++) {
            int x;
            scanf("%d", &x);
            a[i][n + 1] ^= x;  // a[i][n + 1]存储的是第i个灯的状态是否变化
            a[i][i] = 1;  // 默认开关可以控制自己
        }
        
        int x, y;
        while (scanf("%d%d", &x, &y), x || y) a[y][x] = 1;  // 影响的是第y个方程
        
        int t = gauss();
        if (t == -1) puts("Oh,it's impossible~!!");
        else printf("%d\n", t);
    }
    
    return 0;
}
  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值