C++解线性方程--高斯消元法/列主元素消元法/克莱姆法则

一、介绍

解线性方程组是我在大一上线性代数中学习的知识

线性方程组求解项目：

1. 采用VSCode作为IDE，熟悉VSCode下的编码，调试和开发环境配置。
2. 采用C++作为开发语言，并学习使用标准C++库STL。熟练使用vector容器作为动态数组(可变长度数组)。
3. 实现线性方程组求解。至少采用3种求解方式，明白原理，并实践。对比几种不同求解方法的优缺点并编写报告。

1、vector容器介绍

vector是C++标准模板库中的部分内容，中文偶尔译作“容器”，但并不准确。它是一个多功能的，能够操作多种数据结构和算法的模板类和函数库。vector之所以被认为是一个容器，是因为它能够像容器一样存放各种类型的对象，简单地说，vector是一个能够存放任意类型的动态数组，能够增加和压缩数据。

下面将由高斯消元法、列主元素消元法、克莱姆法则求解线性方程组。

二、高斯消元法

1、高斯消元的主要步骤：

1).部分选主元：在每一列中选择绝对值最大的元素所在的行，将其与当前行交换，目的是为了避免除数为0的情况，使得矩阵更容易变成上三角矩阵。

2).消元：从第一行开始，将第i列的第i+1行到第n行的元素通过消元操作变为0，操作包括将当前行的元素减去比例系数乘以上一行对应元素的值。

3).回代求解：从最后一行开始，通过回代的方式逐步求解出方程组的解，即先求出最后一个变量的值，然后带入前面的方程逐步求解出其他变量的值。

2、代码实现

vector<double> Gauss(vector<vector<double>> A, vector<double> b) {//尖括号之间一定要有空格！ 
//vector<double> solveLinearEquations(vector<vector<double> > A, vector<double> b) {
    // 在这里实现解线性方程组的算法
    // 返回解向量

    int n = A.size();
    for (int i = 0; i < n; i++) {
        // 部分选主元：使其形成倒三角
        int maxRow = i;
        for (int j = i + 1; j < n; j++) {
            if (abs(A[j][i]) > abs(A[maxRow][i])) {
                maxRow = j;
            }
        }
        swap(A[i], A[maxRow]);
        swap(b[i], b[maxRow]);

        // 将第i列的第i+1行到第n行消元为0
        for (int j = i + 1; j < n; j++) {
            double ratio = A[j][i] / A[i][i];
            for (int k = i; k < n; k++) {
                A[j][k] -= ratio * A[i][k];
            }
            b[j] -= ratio * b[i];
        }
    }

    // 回代求解
    vector<double> x(n);
    for (int i = n - 1; i >= 0; i--) {
        x[i] = b[i];
        for (int j = i + 1; j < n; j++) {
            x[i] -= A[i][j] * x[j];
        }
        x[i] /= A[i][i];
    }

    return x;
}

三、列主元素消元法

1、列主元素消元法的主要步骤

1).首先，对系数矩阵A进行列主元素选取，即选取每一列中绝对最大的元素所在的列作为主元素列。

2).将选取的主元素列交换到对角线上，并将主元素归一化为1。

3).使用主元素所在列对下方的元素进行消元，将其变为0。

4).最后，进行回代求解，求解出线性方程组的解向量x。

2、代码实现

//列主元素
vector<double> solveLinearEquationslie(vector<vector<double>> A, vector<double> b) 
{
    int n = A.size();
    for (int i = 0; i < n; i++) {
        // 选取绝对值最大的系数所在的列作为主元
        int maxCol = i;
        double maxVal = abs(A[i][i]);
        for (int j = i + 1; j < n; j++) {
            if (abs(A[j][i]) > maxVal) {
                maxCol = j;
                maxVal = abs(A[j][i]);
            }
        }

        // 将主元换到对角线上
        swap(A[i], A[maxCol]);
        swap(b[i], b[maxCol]);

        // 将主元归一化
        double pivot = A[i][i];
        for (int k = i; k < n; k++) {
            A[i][k] /= pivot;
        }
        b[i] /= pivot;

        // 将第i列的第i+1行到第n行消元为0
        for (int j = i + 1; j < n; j++) {
            double ratio = A[j][i];
            for (int k = i; k < n; k++) {
                A[j][k] -= ratio * A[i][k];
            }
            b[j] -= ratio * b[i];
        }
    }
    // 回代求解
    vector<double> x(n);
    for (int i = n - 1; i >= 0; i--) {
        x[i] = b[i];
        for (int j = i + 1; j < n; j++) {
            x[i] -= A[i][j] * x[j];
        }
        x[i] /= A[i][i];
    }

    return x;
}

四、克莱姆法则

1、克莱姆法则的主要步骤

1、检查输入的方程数量是否与未知数数量匹配。
2、创建一个增广矩阵，将系数矩阵A和常数向量b合并成一个矩阵。
3、使用高斯消元法解决方程组：
-对每一列找到主元素（即绝对值最大的元素）并将其作为主元素。
-将主元素所在行与当前行交换，确保主元素在对角线位置。
-将主元素所在行的元素除以主元素，使得主元素变为1。
-将其他行的元素变为0，确保主元素所在列的其他元素变为0。
4、提取解，将解存储在向量x中并返回

2、代码实现

// 克莱姆法则
vector<double> Gramer(vector<vector<double>> A, vector<double> b) {
    int n = A.size();
    vector<double> x;

    // 检查方程数量是否与未知数数量匹配
    if (n != b.size()) {
        cout << "方程数量与未知数数量不匹配" << endl;
        return x;
    }

    // 创建矩阵存储增广矩阵 [A | b]
    vector<vector<double>> augmentedMatrix(n, vector<double>(n + 1));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            augmentedMatrix[i][j] = A[i][j];
        }
        augmentedMatrix[i][n] = b[i];
    }

    // 执行高斯消元法解方程组
    for (int i = 0; i < n; i++) {
        // 找到本列的主元素
        int pivot = i;
        for (int j = i + 1; j < n; j++) {
            if (abs(augmentedMatrix[j][i]) > abs(augmentedMatrix[pivot][i])) {
                pivot = j;
            }
        }

        // 交换行
        swap(augmentedMatrix[i], augmentedMatrix[pivot]);

        // 将对角元素变为1
        double divisor = augmentedMatrix[i][i];
        for (int j = i; j < n + 1; j++) {
            augmentedMatrix[i][j] /= divisor;
        }

        // 将本列其他元素变为0
        for (int j = 0; j < n; j++) {
            if (j != i) {
                double multiplier = augmentedMatrix[j][i];
                for (int k = i; k < n + 1; k++) {
                    augmentedMatrix[j][k] -= multiplier * augmentedMatrix[i][k];
                }
            }
        }
    }

    // 提取解
    for (int i = 0; i < n; i++) {
        x.push_back(augmentedMatrix[i][n]);
    }

    return x;
}

五、小结

经过本次实践训练，我复习了大一所学的线性代数问题，尝试使用vector容器（因为可以动态存储，所以帮我解决了不知道数据容量大小的问题）。刚开始时我用的dev编写程序，但是dev正在使用C++98标准，而在C++98标准中，对于类类型的初始化需要使用构造函数，而不是使用列表初始化。（本方法是升级到支持C++11标准的编译器）所以当时用的：

 vector<vector<double> > A;
double row1[] = {2, 1, -1};
double row2[] = {-3, -1, 2};
double row3[] = {-2, 1, 2};
A.push_back(vector<double>(row1, row1 + 3));
A.push_back(vector<double>(row2, row2 + 3));
A.push_back(vector<double>(row3, row3 + 3));
/*遇到了另一个问题，即no matching function for call to 'std::vector<std::vector<double> >::push_back(<brace-enclosed initializer list>)'。这是因为push_back函数期望接收一个vector<double>类型的参数，而不是一个初始化列表。

要解决这个问题，您可以使用临时变量来存储每一行的系数，并将其添加到A中*/