1. AI技术聚合首页
  2. IT

直接法求解线性方程 （高斯Gauss消去法 矩阵三角分解法） [数值计算方法](c语言描述版)

社会演员多 • • IT • 阅读 28

本系列内容主要针对数值计算方法
第一部分 引言（误差、有效数字、减小误差的方法）

第二部分 非线性方程解法（二分法、牛顿迭代法等）

第三部分 线性方程解法（直接法和迭代法）

第四部分 差值与曲线拟合（拉格朗日、牛顿插值法、最小二乘法等）

第五部分 数值积分与数值微分（牛顿-柯特斯公式，复化积分法等）

第六部分 常微分方程的初值问题的解法

        线性方程解法一般分为直接法和迭代法,这次内容我们主要讲述线性方程的直接法，本文主要是简述多种直接法求解，内容包含理论、代码、以及例题分析.

一、高斯Gauss消去法

        在线性代数中我们曾经学习过克莱姆Cramer法则，这是一种不实用的直接法，本章介绍几种实用的直接法。高斯Gauss消去法是一种规则化的加减消元法，其基本思想是通过逐次消元计算，把一般线性方程组的求解转化为等价的上三角形方程组的求解。

1.顺序高斯消去法

        先来一道简单的三元一次线性方程组：

        那么我们应该如何用代码实现呢？对于这个题目呢，我们首先要录入线性方程的增广矩阵

                                                                [2    4    -2         |         2]

                                                                [1   -3    -3         |        -1]

                                                                [4    2      2         |         3]

        首先我们需要将列出的矩阵导入到代码当中

double A[N][N+1] = {{2, 4, -2, 2}, {1, -3, -3, -1}, {4, 2, 2, 3}};

         然后构造函数计算方法来对矩阵进行处理，具体的完整代码如下：

#include <stdio.h>

#define N 3 // 线性方程组的阶数

void gaussian_elimination(double A[N][N+1]) {
    int i, j, k;
    double factor, temp;

    // 列主元消去
    for (i = 0; i < N - 1; i++) {
        for (j = i + 1; j < N; j++) {
            factor = A[j][i] / A[i][i];
            for (k = i; k <= N; k++) {
                A[j][k] -= factor * A[i][k];
            }
        }
    }

    // 回代求解
    for (i = N - 1; i >= 0; i--) {
        temp = A[i][N];
        for (j = i + 1; j < N; j++) {
            temp -= A[i][j] * A[j][N];
        }
        A[i][N] = temp / A[i][i];
    }
}

int main() {
    double A[N][N+1] = {{2, 4, -2, 2}, {1, -3, -3, -1}, {4, 2, 2, 3}};

    gaussian_elimination(A);

    printf("Solution:\n");
    for (int i = 0; i < N; i++) {
        printf("x%d = %.2f\n", i+1, A[i][N]);
    }

    return 0;
}

        运行上面代码，你会得到Solution:x1 = 0.50 x2 = 0.33 x3 = 0.17

        即（x1 = 1/2         x2 = 1/3         x3 = 1/6）

2.列主元素高斯消元法

        在高斯消去法中，消元过程中可能出现某个元素接近无穷小的情况，这时消去法将无法进行;即使主元素≠0但很小，用其作除数，也会导致其他元素数量级的严重增长和舍入误差的扩散（[数值计算方法] 零基础的学习资料和总结 （ c语言描述版本））关于误差我在这篇文章中有详细的描述。为避免此种情况的发生，可通过交换方程的次序，选取绝对值大的元素作为主元。基于这种想法导出了主元素法。

步骤：

  1. 构造增广矩阵，即系数矩阵A和常数向量b组合在一起。
  2. 对于矩阵的第k列，找到绝对值最大的元素ar,k，并将该行与第k行交换。
  3. 使用这个主元进行消元，使得第k列的第k+1行到底部的元素都变为零。
  4. 重复步骤2和3，直到所有的主对角线上的元素都成为主元，且其下方的元素都是零，此时矩阵成为上三角形式。
  5. 使用回代算法求解上三角形式，得到方程组的解.

例题: 现在让我们看这么一道题  

2x + y - z = 8
-3x - y + 2z = -11
-2x + y + 2z = -3

        使用列主元素高斯消元法，我们首先找到第一列绝对值最大的元素，即-3，然后交换第一行和第二行。接下来，我们对第二列进行消元，使得第二列的第二行和第三行元素变为零。重复这个过程，直到所有列都被消元完毕。最后，通过回代求解出每个变量的值。如下就是代码的实现方式,根据大家的需求进行修改就可以解不同的题目啦!

#include <stdio.h>
#include <math.h>
#define N 3 // 定义矩阵的维度

void pivotGaussianElimination(double A[N][N+1], double x[N]) {
    int i, j, k, r;
    double max, temp, sum;

    // 遍历每一列，进行列主元素高斯消元法
    for (k = 0; k < N - 1; k++) {
        max = fabs(A[k][k]); // 初始化主元素的绝对值为当前列的第一个元素的绝对值
        r = k; // 记录当前主元素所在的行
        // 在当前列中找到绝对值最大的元素，确定主元素
        for (i = k + 1; i < N; i++) {
            if (fabs(A[i][k]) > max) {
                max = fabs(A[i][k]);
                r = i;
            }
        }
        // 将主元素所在的行与当前行交换
        for (j = k; j <= N; j++) {
            temp = A[k][j];
            A[k][j] = A[r][j];
            A[r][j] = temp;
        }

        // 将主元素下方的元素消为零
        for (i = k + 1; i < N; i++) {
            A[i][k] /= A[k][k];
            for (j = k + 1; j <= N; j++) {
                A[i][j] -= A[k][j] * A[i][k];
            }
        }
    }

    // 回代求解未知变量
    x[N - 1] = A[N - 1][N] / A[N - 1][N - 1];
    for (i = N - 2; i >= 0; i--) {
        sum = 0.0;
        for (j = i + 1; j < N; j++) {
            sum += A[i][j] * x[j];
        }
        x[i] = (A[i][N] - sum) / A[i][i];
    }
}

int main() {
    double A[N][N+1] = {
        {2, 1, -1, 8},
        {-3, -1, 2, -11},
        {-2, 1, 2, -3}
    };
    double x[N];

    pivotGaussianElimination(A, x);

    printf("Solution:\n");
    for (int i = 0; i < N; i++) {
        printf("x%d = %f\n", i + 1, x[i]);
    }

    return 0;
}

        列主元Gauss消去法是在每一步消元前，在主元所在的一列选取绝对值最大的元素作为主元素。而全主元Gauss消去法是在每一步消元前，在所有元素中选取绝对值最大的元素作为主元素。但由于运算量大增，实际应用中并不经常使用。

3.高斯-若尔当消去法

        高斯-若尔当消去法是高斯消去法的扩展，它不仅将线性方程组转换为上三角形式，而且进一步通过行操作将其转换为简化行阶梯形（Reduced Row Echelon Form, RREF），即对角线上的元素都是1，而对角线以上的元素都是0的矩阵。这种方法可以精确地解决线性方程组，而不受舍入误差的影响。

步骤：

  1. 构造增广矩阵，将系数矩阵A和常数向量b组合在一起。
  2. 对于矩阵的第k列，找到绝对值最大的元素ar,k，并将该行与第k行交换。
  3. 使用这个主元进行消元，不仅使第k列的第k+1行到底部的元素都变为零，而且还要使第k行的主元变为1。
  4. 对于第k行，除以主元使它变为1，然后使用这一行去消去第k+1行到第n行的第k列元素。
  5. 重复步骤2到4，直到所有的对角线上的元素都成为1，且对角线以上的元素都是零。

例题: 

2x - 3y + z = 4
-6x + 9y - 3z = -12
x - y + 0z = 1

首先，我们构造增广矩阵：

|  2  -3   1 |  4 |
| -6   9  -3 | -12|
|  1  -1   0 |  1 |

 通过行交换和消元操作，我们逐步将矩阵转换为简化行阶梯形：

|  1  -1   0 |  1 |
|  0   3  -1 |  2 |
|  0   0   1 |  4 |

 观察上面矩阵的第三行可以容易的得出Z=4,带入第二行得Y=2,最后计算第一行得出X=1。

下面我们给出具体代码(我已经在关键的步骤上写了中文注释):

#include <stdio.h>
#include <math.h>

#define MAX_SIZE 100

void gaussianJordanElimination(double A[][MAX_SIZE], int n, double x[]) {
    int i, j, k, max_row;
    double max_val, temp;

    // 遍历每一列进行消元
    for (k = 0; k < n - 1; k++) {
        // 找到主元素所在的行并交换到当前行
        max_val = fabs(A[k][k]);
        max_row = k;
        for (i = k + 1; i < n; i++) {
            if (fabs(A[i][k]) > max_val) {
                max_val = fabs(A[i][k]);
                max_row = i;
            }
        }
        // 交换行 k 和 max_row
        for (j = 0; j <= n; j++) {
            temp = A[k][j];
            A[k][j] = A[max_row][j];
            A[max_row][j] = temp;
        }

        // 将对角线元素化为1，将对角线上方和下方的元素化为0
        for (i = 0; i < n; i++) {
            if (i != k) {
                A[i][k] -= A[k][k] * A[i][k];
            }
        }
        A[k][k] = 1.0 / A[k][k]; // 对角线元素归一化为1
        for (j = k + 1; j <= n; j--) {
            A[k][j] *= A[k][k];
        }
    }

    // 最后一行自动成为单位行，直接将解复制过来
    for (i = 0; i < n; i++) {
        x[i] = A[n - 1][i];
    }
}

int main() {
    int n;
    printf("输入线性方程组的阶数: ");
    scanf("%d", &n);
    double A[MAX_SIZE][MAX_SIZE], x[n];

    // 输入增广矩阵
    printf("输入增广矩阵:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j <= n; j++) {
            scanf("%lf", &A[i][j]);
        }
    }

    gaussianJordanElimination(A, n, x);

    // 输出解向量
    printf("解向量为:\n");
    for (int i = 0; i < n; i++) {
        printf("x%d = %f\n", i + 1, x[i]);
    }

    return 0;
}

二、矩阵三角分解法

1.Gauss消去法的矩阵表示

步骤：

  1. 初始化一个单位矩阵作为下三角矩阵的初始状态。
  2. 对于每一步骤的高斯消去，确定需要的行变换（例如，交换行、行乘以标量、行相加）。
  3. 构建一个初等矩阵来表示该步的行变换。
  4. 将这个初等矩阵乘以当前的下三角矩阵，得到新的下三角矩阵。
  5. 应用相同的行变换到原始矩阵上，得到上三角矩阵。
  6. 不断循环步骤2到5，直到原始矩阵完全转换为上三角形式。

 例题: 接下来我们增加一下例题难度，来看下面例题：

3x + 4y - 5z + 2w = 1
x - 2y + 3z - 4w = 2
-6x + 8y - 3z + 7w = 3
9x - 4y + z - 2w = 4

        让我们分步看看这道题：首先依旧是把系数提取出来构造增广矩阵 

int main() {
    double A[N][N+1] = {
        {3, 4, -5, 2, 1},
        {1, -2, 3, -4, 2},
        {-6, 8, -3, 7, 3},
        {9, -4, 1, -2, 4}
    };

         接下来构造函数，在每一步中，我们寻找当前列绝对值最大的元素作为主元，并将其所在的行与该步的初始行交换，以确保每一步的消元过程都是有效的。利用主元，我们通过行变换将当前列的其他行中的对应元素消为零。当矩阵转换为行阶梯形后，我们从最后一行开始，逐步回代求解每个变量的值。最终我们可以构造出下面完整代码。

#include <stdio.h>

#define N 4

// 打印矩阵
void printMatrix(double A[N][N+1]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N+1; j++) {
            printf("%f ", A[i][j]);
        }
        printf("\n");
    }
}

// 高斯消去法求解线性方程组
void gaussianElimination(double A[N][N+1]) {
    int i, j, k, max_row;
    double max_val, multiplier;

    // 进行消元操作，将矩阵转化为上三角形式
    for (k = 0; k < N - 1; k++) {
        max_val = fabs(A[k][k]);
        max_row = k;

        // 找到最大元素所在的行，进行主元素交换
        for (i = k + 1; i < N; i++) {
            if (fabs(A[i][k]) > max_val) {
                max_val = fabs(A[i][k]);
                max_row = i;
            }
        }

        // 交换最大元素所在的行和当前行
        for (j = k; j <= N; j++) {
            double temp = A[max_row][j];
            A[max_row][j] = A[k][j];
            A[k][j] = temp;
        }

        // 使用当前行将下方元素消为0
        for (i = k + 1; i < N; i++) {
            multiplier = A[i][k] / A[k][k];

            // 执行行操作：A[i] = A[i] - multiplier * A[k]
            for (j = k; j <= N; j++) {
                A[i][j] -= multiplier * A[k][j];
            }
        }
    }

    // 矩阵现在已经处于行阶梯形式。进行回代以找到解
    for (i = N - 1; i >= 0; i--) {
        for (j = i + 1; j < N; j++) {
            A[i][N] -= A[i][j] * A[j][N];
        }
        A[i][N] /= A[i][i];
    }
}

int main() {
    double A[N][N+1] = {
        {3, 4, -5, 2, 1},
        {1, -2, 3, -4, 2},
        {-6, 8, -3, 7, 3},
        {9, -4, 1, -2, 4}
    };

    printf("原始增广矩阵:\n");
    printMatrix(A);

    gaussianElimination(A);

    printf("行阶梯形式矩阵:\n");
    printMatrix(A);

    return 0;
}

2.杜里特尔 Doolittle分解法

步骤：

  1. 首先，我们定义并输入系数矩阵A。
  2. 通过Doolittle分解法，我们将矩阵A分解为下三角矩阵L和上三角矩阵U
  3. 将得到的下三角矩阵L和上三角矩阵U回代求解

例题：依旧是一小道

x + 3y + 7z = 6
4x + 2y + 6z = 13
5x + 9y + z = 14

        根据我给出的步骤也不难得出下面代码：

#include <stdio.h>

#define N 3

// Doolittle分解法进行LU分解
void DoolittleDecomposition(double A[N][N], double L[N][N], double U[N][N]) {
    // 初始化L矩阵为对角线为1，其余为0
    for (int i = 0; i < N; i++) {
        // 对角线元素为1
        for (int j = 0; j < (i + 1); j++) {
            L[i][j] = 0;
        }
        L[i][i] = 1; // 对角线元素为1

        // 初始化U矩阵为0
        for (int j = i; j < N; j++) {
            U[i][j] = 0;
        }
    }

    // 执行分解
    for (int j = 0; j < N; j++) {
        for (int i = j; i < N; i++) {
            // 计算U矩阵的元素
            if (j == i) {
                U[i][j] = A[i][j]; // 主对角线元素直接赋值
            } else {
                double sum = 0;
                // 计算累加和
                for (int k = 0; k < j; k++) {
                    sum += L[i][k] * U[k][j];
                }
                U[i][j] = A[i][j] - sum; // 计算U矩阵的非对角线元素
            }
        }
    }
}

// 使用Doolittle分解得到的L和U矩阵进行回代求解方程组
void solveEquations(double L[N][N], double U[N][N], double b[N], double x[N]) {
    double y[N];

    // 前向替换求解Ly = b
    for (int i = 0; i < N; i++) {
        double sum = 0;
        for (int j = 0; j < i; j++) {
            sum += L[i][j] * y[j];
        }
        y[i] = b[i] - sum;
    }

    // 回代求解Ux = y
    for (int i = N - 1; i >= 0; i--) {
        double sum = 0;
        for (int j = i + 1; j < N; j++) {
            sum += U[i][j] * x[j];
        }
        x[i] = (y[i] - sum) / U[i][i];
    }
}

int main() {
    double A[N][N] = {
        {1, 3, 7},
        {4, 2, 6},
        {5, 9, 1}
    };
    double L[N][N], U[N][N], b[N], x[N];

    // 设置右端向量b
    b[0] = 1;
    b[1] = 2;
    b[2] = 3;

    // 进行Doolittle分解
    DoolittleDecomposition(A, L, U);

    // 求解方程组Ax = b
    solveEquations(L, U, b, x);

    // 打印解向量x
    printf("解向量x:\n");
    for (int i = 0; i < N; i++) {
        printf("%f\n", x[i]);
    }

    return 0;
}

        以上就是本次的全部内容了，希望给大家带来指导和帮助。

版权声明：本文为博主作者：OnlyOswald原创文章，版权归属原作者，如果侵权，请联系我们删除！

原文链接：https://blog.csdn.net/m0_74431077/article/details/138195989

共计人评分，平均分

到目前为止还没有投票！成为第一位评论此文章。

赞 (0)
社会演员多的头像社会演员多普通用户
0
上一篇 2024年5月6日
下一篇 2024年5月6日

相关推荐