列主元高斯消去法解线性方程组——C语言实现

原理

高斯消去法

的基本原理就是用初等变换将用行的，逐次消去未知数的方法，把原来的方程组，化为与其等价的上三角方程组。
设有线性方程$\mathbf{A}\mathbf{x}=\mathbf{B}$,写成矩阵如下
$$\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]$$

该线性方程组的初始增广矩阵为

$$\left[\begin{array}{ccccc}{a}_{11}^{(0)}& a_{12}^{(0)}& \cdots & a_{1n}^{(0)}& b_1^{(0)}\\ a_{21}^{(0)}& a_{22}^{(0)}& \cdots & a_{2n}^{(0)}& b_2^{(0)}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n1}^{(0)}& a_{n2}^{(0)}& \cdots & a_{mn}^{(0)}& b_n^{(0)}\end{array}\right]$$
利用初等行变换的方法，即消元，将其化为如下形式

$$\left[\begin{array}{ccccc}1& a_{12}^{\left(n\right)}& \cdots & a_{1n}^{\left(n\right)}& b_1^{\left(n\right)}\\ 0& 1& \cdots & a_{2n}^{\left(n\right)}& b_2^{\left(n\right)}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& b_n^{\left(n\right)}\end{array}\right]$$
其上标$n$表示进行第$n$消元后的值

列主元高斯消去法

列主元高斯消去法较高斯消去法多了一步选取该列主元，即选取该列上绝对值最大的元素作为主元的步骤，交换该行与主元所在的行。目的是使得高斯消去法更稳定。举例说明，比如说现在有一个线性方程组$\mathbf{A}\mathbf{x}=\mathbf{B}$如下
$$\begin{array}{l}2x_1-x_2+3x_3=1\\ 4x_1+2x_2+5x_3=4\\ x_1+2x_2=7\end{array}$$
其增广矩阵为
$$\left[\begin{array}{cccc}2& -1& 3& 1\\ 4& 2& 5& 4\\ 1& 2& 0& 7\end{array}\right]$$
在进行第一步消元时，先选取第一列中绝对值最大的元素4作为主元，交换第一行和第二行，得到
$$\left[\begin{array}{cccc}4& 2& 5& 4\\ 2& -1& 3& 1\\ 1& 2& 0& 7\end{array}\right]$$
再进行消元，得到
$$\left[\begin{array}{cccc}4& 2& 5& 4\\ 0& -2& 0.5& -3\\ 0& 1.5& -1.25& 6\end{array}\right]$$

算法

1. ∃ l ∈ { k , k + 1 , … , n } a l k = max ⁡ k ≤ i ≤ n a i k , a l j ⇔ a k j , j = k , k + 1 , … , n \exists l \in \{k,k+1,\ldots,n\}\\a_{lk}=\max\limits_{k\le i \le n}a_{ik},a_{lj}\Leftrightarrow a_{kj},j=k,k+1,\ldots,n ∃l∈{k,k+1,…,n}alk​=k≤i≤nmax​aik​,alj​⇔akj​,j=k,k+1,…,n

2. $$\begin{gathered} a_{kj}^{(k)}=a_{kj}^{(k-1)}/a_{kk}^{(k-1)} \\ {a}_{ij}^{(k)}=a_{ij}^{(k-1)}-a_{ik}^{(k-1)}{a}_{kj}^{(k)},i,j=k,k+1,\dots ,n \end{gathered}$$

3. $\begin{array}{l}{x}_n=a_{n,n+1}^{(n)},i=n\\ x_i=a_{i,n+1}^{(n)}-{\sum }_{j=i+1}^n{a}_{ij}^{(n)}{x}_j,i=n-1,n-2,\cdots ,1\end{array}$

c程序

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

void gauss(double[], double[], double[], int);

int main(int argc, char* argv[])
{
    int n           = 4;
    double a[4 * 4] = { 1, 1, 0, 3,
                        2, 1, -1, 1,
                        3, -1, -1, 3,
                        -1, 2, 3, -1 };
    double b[4]     = { 4, 1, -3, 4 };
    double x[4]     = {};

    printf("增广矩阵\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%f\t", a[i * n + j]);
        }
        printf("%f\n", b[i]);
    }
    printf("解向量\n");
    gauss(a, b, x, n);
    for (int i = 0; i < n; i++)
        printf("%f\t", x[i]);
    printf("\n");
    return 0;
}

void gauss(double a[], double b[], double x[], int n)//列主元高斯消去法，a[]系数矩阵，b[]结果向量，x[]解向量，n阶数
{
    int i, j, k, exchangeline, exchangeflag = 0;
    double temp, max;

    for (k = 0; k < n - 1; k++) { //k迭代次数
        max = a[k * n + k];
        for (i = k + 1; i < n; i++) { //寻找主元，i行号
            if (fabs(max) < fabs(a[n * i + k])) {
                max          = a[n * i + k];
                exchangeflag = 1;     //交换标志
                exchangeline = n * i; //记录需要交换的行号
            }
        }
        if (exchangeflag) { //换行，j列号
            for (j = 0; j < n; j++) {
                temp                = a[exchangeline + j]; //对系数矩阵操作
                a[exchangeline + j] = a[k * n + j];
                a[k * n + j]        = temp;
            }
            temp                = b[exchangeline / n]; //对结果向量操作
            b[exchangeline / n] = b[k];
            b[k]                = temp;
            exchangeflag        = 0; //清除交换标志
        }
        for (i = k + 1; i < n; i++) { //消元
            temp = a[i * n + k] / a[k * n + k];
            b[i] = b[i] - b[k] * temp; //对结果向量操作
            for (j = k; j < n; j++)
                a[i * n + j] = a[i * n + j] - a[k * n + j] * temp; //对系数矩阵操作
        }
    }
    for (i = n - 1; i > -1; i--) { //回代
        temp = b[i];
        for (j = n - 1; j > i; j--)
            temp = temp - a[n * i + j] * x[j];
        x[i] = temp / a[i * n + i];
    }
}

运行结果