hiho一下第56周高斯消元

最新推荐文章于 2021-09-26 11:39:22 发布

fisty

最新推荐文章于 2021-09-26 11:39:22 发布

阅读量499

点赞数

分类专栏：模板

本文链接：https://blog.youkuaiyun.com/u013008291/article/details/47779271

版权

模板同时被 2 个专栏收录

7 篇文章

订阅专栏

hihocoder

2 篇文章

订阅专栏

本文详细介绍了如何使用一元方程组解决实际问题，并深入讲解了高斯消元法的基本原理及具体应用过程。通过实例演示了如何通过一系列行变换将方程组转化为易于求解的形式。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

小Ho：<吧唧><吧唧><吧唧>

小Hi：小Ho，你还吃呢。想好了么?

小Ho：肿抢着呢(正想着呢)...<吞咽>...我记得这个问题上课有提到过，应该是一元一次方程组吧。

我们把每一件商品的价格看作是x[1]..x[n]，第i个组合中第j件商品数量记为a[i][j]，其价格记作y[i]，则可以列出方程式：

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1]
a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2]
                                  ...
a[m][1] * x[1] + a[m][2] * x[2] + ... + a[m][n] * x[n] = y[m]

我们可以对方程组进行3种操作而不改变方程组的解集：

1. 交换两行。

2. 把第i行乘以一个非0系数k。即对于j = 1..n，令a[i][j] = k*a[i][j], y[i]=k*y[i]

3. 把第p行乘以一个非0系数k之后加在第i行上。即对于j=1..n，令a[i][j] = a[i][j]+k*a[p][j]，y[i]=y[i]+k*p[i]

以上三个操作叫做初等行变换。

我们可以使用它们，对这个方程组中的a[i][j]进行加减乘除变换，举个例子：

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y[1]    式子(1)
a[2][1] * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y[2]    式子(2)

我们可以通过式子(1) - 式子(2) * (a[1][1] / a[2][1])，将第1行第1列的a[1][1]变换为0。

对整个方程组进行多次变换之后，可以使得a[i][j]满足：

a[i][j] = 1 (i == j)
a[i][j] = 0 (i != j)

则整个方程组变成了：

x[1] = y'[1]
x[2] = y'[2]
...
x[n] = y'[n]
0 = y'[n + 1]
0 = y'[n + 2]
...
0 = y'[m]

这样的话，y'[1] .. y'[n]就是我们要求的x[1]..x[n]

小Hi：挺不错的嘛，继续？

小Ho：好，关于如何变换，我们可以利用一个叫高斯消元的算法。高斯消元分成了2个步骤：

首先我们要计算出上三角矩阵，也就是将方程组变为：

a[1][1] * x[1] + a[1][2] * x[2] + ... + a[1][n] * x[n] = y'[1]
      0 * x[1] + a[2][2] * x[2] + ... + a[2][n] * x[n] = y'[2]
      0 * x[1] +       0 * x[2] + ... + a[3][n] * x[n] = y'[3]
                                   ...
      0 * x[1] +       0 * x[2] + ... + a[n][n] * x[n] = y'[n]
      0 * x[1] +       0 * x[2] + ... +       0 * x[n] = y'[n + 1]
	                               ...
      0 * x[1] +       0 * x[2] + ... +       0 * x[n] = y'[m]

也就是通过变换，将所有a[i][j](i>j)变换为0。同时要保证对角线上的元素a[i][i]不为0。

方法也很见简单，从第1行开始，我们利用当前行第i列不为0，就可以通过变换将i+1..M行第一列全部变换为0，接着对于第2行，我们用同样的方法将第3..M行第2列也变换为0...不断重复直到第n行为止。

假如计算到第i行时，第i列已经为0，则我们需要在第i+1..M行中找到一行第i列不为0的行k，并交换第i行和第k行，来保证a[i][i] != 0。但这时候还有可能出现一个情况，就是第i..M行中的i列均为0，此时可以判定，该方程组有多解。

当得到上三角矩阵后，就可以从第n行开始逆推，一步一步将a[i][j](i<j)也变换为0.

因为第n行为a[n][n] * x[n] = y'[n]，则x[n] = y'[n] / a[n][n]。

第n-1行为a[n-1][n-1] * x[n - 1] + a[n][n] * x[n] = y'[n - 1]。我们将得到的x[n]代入，即可计算出x[n-1]。

同样的依次类推就可以得到所有的x[1]..x[n]。

而对于多解和无解的判定：

当在求出的上三角矩阵中出现了 a[i][1] = a[i][2] = ... = a[i][n] = 0, 但是y'[i] != 0时，产生了矛盾，即出现了无解的情况。

而多解的证明如下：

假设n=3,m=3，而我们计算出了上三角矩阵为：

a * x[1] + b * x[2] + c * x[3] = d
                      e * x[3] = f
                             0 = 0

当我们在第一个式子中消去x[3]后，有a * x[1] + b * x[2] = g，显然x[1]和x[2]有无穷多种可能的取值。

小Hi：既然小Ho你都已经把整个算法讲了，那么我就只能给出伪代码了：

// 处理出上三角矩阵
For i = 1..N
    Flag ← False
    For j = i..M                // 从第i行开始，找到第i列不等于0的行j
        If a[j][i] != 0 Then
            Swap(j, i)          // 交换第i行和第j行
            Flag ← True
            Break
        End If
    End For
    // 若无法找到，则存在多个解
    If (not Flag) Then
        manySolutionsFlag ← True // 出现了秩小于N的情况
        continue;
    End If
    // 消除第i+1行到第M行的第i列
    For j = i+1 .. M
        a[j][] ← a[j][] - a[i][] * (a[j][i] / a[i][i])
        b[j] ← b[j] - b[i] * (a[j][i] / a[i][i])
    End For
End For 

// 检查是否无解，即存在 0 = x 的情况
For i = 1..M
    If (第i行系数均为0 and (b[i] != 0)) Then
        return "No solutions"
    End If
End For

// 判定多解
If (manySolutionsFlag) Then
	return "Many solutions"
End If

// 此时存在唯一解
// 由于每一行都比前一行少一个系数，所以在M行中只有前N行有系数
// 解析来从第N行开始处理每一行的解
For i = N .. 1
    // 利用已经计算出的结果，将第i行中第i+1列至第N列的系数消除
    For j = i + 1 .. N
        b[i] ← b[i] - a[i][j] * value[j]
        a[i][j] ← 0
    End For
    value[i] ← b[i] / a[i][i]
End For

那最后能够拜托你实现一下这个算法么？

小Ho：没问题，等我吃完这包薯片就去！

/***********************************************
 * Author: fisty
 * Created Time: 2015-08-19 13:22:23
 * File Name   : 56_Gauss.cpp
 *********************************************** */
#include <iostream>
#include <cstring>
#include <deque>
#include <cmath>
#include <queue>
#include <stack>
#include <list>
#include <map>
#include <set>
#include <string>
#include <vector>
#include <cstdio>
#include <bitset>
#include <algorithm>
using namespace std;
#define Debug(x) cout << #x << " " << x <<endl
#define Memset(x, a) memset(x, a, sizeof(x))
const int INF = 0x3f3f3f3f;
typedef long long LL;
typedef pair<int, int> P;
#define FOR(i, a, b) for(int i = a;i < b; i++)
#define lson l, m, k<<1
#define rson m+1, r, k<<1|1
    
const double eps = 1e-8;  
const int Max_M = 1000;       ///m个方程,n个变量  
const int Max_N = 500;  
int m, n;  
double Aug[Max_M][Max_N+1]; ///增广矩阵  
bool free_x[Max_N];         ///判断是否是不确定的变元  
double x[Max_N];            ///解集  

int sign(double x){ return (x>eps) - (x<-eps);}  

/** 
  返回值: 
  -1 无解 
  0 有且仅有一个解 
  >=1 有多个解,根据free_x判断哪些是不确定的解 
  */  
int Gauss()  
{  
    int i,j;  
    int row,col,max_r;  
    for(row=0,col=0; row<m&&col<n; row++,col++)  
    {  
        max_r = row;  
        for(i = row+1; i < m; i++)  ///找到当前列所有行中的最大值(做除法时减小误差)  
        {  
            if(sign(fabs(Aug[i][col])-fabs(Aug[max_r][col]))>0)  
                max_r = i;  
        }  
        if(max_r != row)            ///将该行与当前行交换  
        {  
            for(j = row; j < n+1; j++)  
                swap(Aug[max_r][j],Aug[row][j]);  
        }  
        if(sign(Aug[row][col])==0)  ///当前列row行以下全为0(包括row行)  
        {  
            row--;  
            continue;  
        }  
        for(i = row+1; i < m; i++)  
        {  
            if(sign(Aug[i][col])==0)  
                continue;  
            double ta = Aug[i][col]/Aug[row][col];  
            for(j = col; j < n+1; j++)  
                Aug[i][j] -= Aug[row][j]*ta;  
        }  
    }  
    for(i = row; i < m; i++)    ///col=n存在0...0,a的情况,无解  
    {  
        if(sign(Aug[i][col]))  
            return -1;  
    }  
    if(row < n)     ///存在0...0,0的情况,有多个解,自由变元个数为n-row个  
    {  
        for(i = row-1; i >=0; i--)  
        {  
            int free_num = 0;   ///自由变元的个数  
            int free_index;     ///自由变元的序号  
            for(j = 0; j < n; j++)  
            {  
                if(sign(Aug[i][j])!=0 && free_x[j])  
                    free_num++,free_index=j;  
            }  
            if(free_num > 1) continue;  ///该行中的不确定的变元的个数超过1个,无法求解,它们仍然为不确定的变元  
            ///只有一个不确定的变元free_index,可以求解出该变元,且该变元是确定的  
            double tmp = Aug[i][n];  
            for(j = 0; j < n; j++)  
            {  
                if(sign(Aug[i][j])!=0 && j!=free_index)  
                    tmp -= Aug[i][j]*x[j];  
            }  
            x[free_index] = tmp/Aug[i][free_index];  
            free_x[free_index] = false;  
        }  
        return n-row;  
    }  
    ///有且仅有一个解,严格的上三角矩阵(n==m)  
    for(i = n-1; i >= 0; i--)  
    {  
        double tmp = Aug[i][n];  
        for(j = i+1; j < n; j++)  
            if(sign(Aug[i][j])!=0)  
                tmp -= Aug[i][j]*x[j];  
        x[i] = tmp/Aug[i][i];  
    }  
    return 0;  
}  

void init(){
    memset(Aug,0.0,sizeof(Aug));  
    memset(x,0.0,sizeof(x));  
    memset(free_x,1,sizeof(free_x));    ///都是不确定的变元   
}
int main(){
    scanf("%d%d", &n, &m);
    init();
    int ans;
    for(int i = 0;i < m; i++){
        for(int j = 0;j <= n; j++){
            scanf("%lf", &Aug[i][j]);
        }
    }
    ans = Gauss();
    if(ans == -1){
        puts("No solutions");
    }else if(ans >= 1){
        puts("Many solutions");
    }else{
        for(int i = 0;i < n; i++){
            printf("%d\n", int(x[i] + eps));
        }
    }
}