高斯消元法 (模板)

/* 用于求整数解得方程组. */
#include <iostream>
#include <string>
#include <cmath>
using namespace std;
const int maxn = 105;  //矩阵大小
int equ, var; // 有equ个方程,var个变元。增广阵行数为equ, 分别为0到equ - 1,列数为var + 1,分别为0到var.
int a[maxn][maxn];
int x[maxn]; // 解集.
bool free_x[maxn]; // 判断是否是不确定的变元.
int free_num;
void Debug(void)
{
    int i, j;
    for (i = 0; i < equ; i++)
    {
        for (j = 0; j < var + 1; j++)
        {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
    cout << endl;
}
inline int gcd(int a, int b)
{
    int t;
    while (b != 0)
    {
        t = b;
        b = a % b;
        a = t;
    }
    return a;
}
inline int lcm(int a, int b)
{
    return a * b / gcd(a, b);
}
// 高斯消元法解方程组(Gauss-Jordan elimination).(-2表示有浮点数解,但无整数解,-1表示无解,0表示唯一解,大于0表示无穷解,并返回自由变元的个数)
int Gauss(void)
{
    int i, j, k;
    int max_r; // 当前这列绝对值最大的行.
    int col; // 当前处理的列.
    int ta, tb;
    int LCM;
    int temp;
    int free_x_num;
    int free_index;
    // 转换为阶梯阵.
    col = 0; // 当前处理的列.
    for (k = 0; k < equ && col < var; k++, col++) // 枚举当前处理的行. // 找到该col列元素绝对值最大的那行与第k行交换.(为了在除法时减小误差)
    {
        max_r = k;
        for (i = k + 1; i < equ; i++)
        {
            if (abs(a[i][col]) > abs(a[max_r][col])) max_r = i;
        }
        if (max_r != k)//与第K行交换
        { 
            for (j = k; j < var + 1; j++) swap(a[k][j], a[max_r][j]);
        }
        if (a[k][col] == 0) //说明该col列第k行以下全是0了,则处理当前行的下一列。
        { 
            k--; continue;
        }
        for (i = k + 1; i < equ; i++)//枚举要删去的行。
        { 
            if (a[i][col] != 0)
            {
                LCM = lcm(abs(a[i][col]), abs(a[k][col]));
                ta = LCM / abs(a[i][col]), tb = LCM / abs(a[k][col]);
                if (a[i][col] * a[k][col] < 0) tb = -tb; // 异号的情况是两个数相加.
                for (j = col; j < var + 1; j++)
                {
                    a[i][j] = a[i][j] * ta - a[k][j] * tb;
                }
            }
        }
    }
    Debug();
    // 1. 无解的情况: 化简的增广阵中存在(0, 0, ..., a)这样的行(a != 0).
    for (i = k; i < equ; i++) // 对于无穷解来说,如果要判断哪些是自由变元,那么初等行变换中的交换就会影响,则要记录交换
    { 
        if (a[i][col] != 0) return -1;
    }
    // 2. 无穷解的情况: 在var * (var + 1)的增广阵中出现(0, 0, ..., 0)这样的行,即说明没有形成严格的上三角阵.且出现的行数即为自由变元的个数.
    if (k < var)
    {
        // 首先,自由变元有var - k个,即不确定的变元至少有var - k个.
        for (i = k - 1; i >= 0; i--)
        {
            // 第i行一定不会是(0, 0, ..., 0)的情况,因为这样的行是在第k行到第equ行.
            // 同样,第i行一定不会是(0, 0, ..., a), a != 0的情况,这样的无解的.
            free_x_num = 0; // 用于判断该行中的不确定的变元的个数,如果超过1个,则无法求解,它们仍然为不确定的变元.
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && free_x[j]) free_x_num++, free_index = j;
            }
            if (free_x_num > 1) continue; // 无法求解出确定的变元.
            // 说明就只有一个不确定的变元free_index,那么可以求解出该变元,且该变元是确定的.
            temp = a[i][var];
            for (j = 0; j < var; j++)
            {
                if (a[i][j] != 0 && j != free_index) temp -= a[i][j] * x[j];
            }
            x[free_index] = temp / a[i][free_index]; // 求出该变元.
            free_x[free_index] = 0; // 该变元是确定的.
        }
        return var - k; // 自由变元有var - k个.
    }
    // 3. 唯一解的情况: 在var * (var + 1)的增广阵中形成严格的上三角阵.
    // 计算出Xn-1, Xn-2 ... X0.
    for (i = var - 1; i >= 0; i--)
    {
        temp = a[i][var];
        for (j = i + 1; j < var; j++)
        {
            if (a[i][j] != 0) temp -= a[i][j] * x[j];
        }
        if (temp % a[i][i] != 0) return -2; // 说明有浮点数解,但无整数解.
        x[i] = temp / a[i][i];
    }
    return 0;
}
int main(void)
{
    freopen("Input.txt", "r", stdin);
    int i, j;
    while (scanf("%d %d", &equ, &var) != EOF)
    {
        memset(a, 0, sizeof(a));
        memset(x, 0, sizeof(x));
        memset(free_x, 1, sizeof(free_x)); // 一开始全是不确定的变元.
        for (i = 0; i < equ; i++)
        {
            for (j = 0; j < var + 1; j++)
            {
                scanf("%d", &a[i][j]);
            }
        }
        //        Debug();
        free_num = Gauss();
        if (free_num == -1) printf("无解!\n");
        else if (free_num == -2) printf("有浮点数解,无整数解!\n");
        else if (free_num > 0)
        {
            printf("无穷多解! 自由变元个数为%d\n", free_num);
            for (i = 0; i < var; i++)
            {
                if (free_x[i]) printf("x%d 是不确定的\n", i + 1);
                else printf("x%d: %d\n", i + 1, x[i]);
            }
        }
        else
        {
            for (i = 0; i < var; i++)
            {
                printf("x%d: %d\n", i + 1, x[i]);
            }
        }
        printf("\n");
    }
    return 0;
}
/* czyuan原创,转载请注明出处。*/
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
高斯消元法是一种线性代数中求解线性方程组的常用方法,它通过一系列的行变换将增广矩阵转化为上三角矩阵,从而求出方程组的解。 下面我们来详细介绍高斯消元法的步骤: 1. 将增广矩阵写成一个矩阵形式:$$\left[\begin{array}{cccc|c}a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\\vdots & \vdots & \ddots & \vdots & \vdots \\a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \\\end{array}\right]$$ 2. 选择第一列系数最大的行作为第一步的主元素所在的行,将该行移到第一行。 3. 通过消元操作将第一列的其他元素变成零。具体地,对于第 $i$ 行,我们将其乘以一个系数 $m_i$ 加到第一行上,使得第一行的第 $i$ 个元素变成零。需要注意的是,系数 $m_i$ 的取值为:$$m_i = -\frac{a_{i1}}{a_{11}}$$ 4. 将第二列系数最大的行作为第二步的主元素所在的行,将该行移到第二行。 5. 通过消元操作将第二列的其他元素变成零。具体地,对于第 $i$ 行,我们将其乘以一个系数 $m_i$ 加到第二行上,使得第二行的第 $i$ 个元素变成零。需要注意的是,系数 $m_i$ 的取值为:$$m_i = -\frac{a_{i2}}{a_{22}}$$ 6. 重复上述步骤,直到将增广矩阵转化为上三角矩阵。此时,方程组的解可以通过回代得到。 7. 回代过程:从最后一行开始,依次求解每个未知量。具体地,对于第 $i$ 个未知量,我们先将第 $i$ 行的解代入第 $i$ 个方程中,然后依次代入已知的第 $i+1$ 到第 $n$ 个未知量的解,得到第 $i$ 个未知量的解。 下面是高斯消元法的代码实现: ```c++ const double eps = 1e-8; int gauss(vector<vector<double>>& a, vector<double>& b) { int n = a.size(); int m = a[0].size() - 1; vector<int> p(n); for (int i = 0; i < n; i++) { p[i] = i; } for (int k = 0; k < m; k++) { int pivot = k; for (int i = k; i < n; i++) { if (abs(a[i][k]) > abs(a[pivot][k])) { pivot = i; } } swap(a[pivot], a[k]); swap(b[pivot], b[k]); if (abs(a[k][k]) < eps) { return -1; } for (int i = k + 1; i < n; i++) { double f = a[i][k] / a[k][k]; b[i] -= f * b[k]; for (int j = k; j < m; j++) { a[i][j] -= f * a[k][j]; } } } vector<double> x(m); for (int k = m - 1; k >= 0; k--) { x[k] = b[k]; for (int i = k + 1; i < m; i++) { x[k] -= a[k][i] * x[i]; } x[k] /= a[k][k]; } return 0; } ``` 其中,输入参数为一个 $n \times (m+1)$ 的增广矩阵 $A$ 和一个长度为 $n$ 的向量 $b$,输出为 $0$ 或者 $-1$,表示方程组有唯一解或者无解,解存储在长度为 $m$ 的向量 $x$ 中。 需要注意的是,为了防止精度误差,我们在进行消元操作时,如果某个数的绝对值小于一个极小值 $\epsilon$,则将其视为零。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值