数值分析：使用python实现列主元消去法、原始高斯消去法

1 实验目的

       正确理解原始高斯消去法，并清楚其优缺点，正确理解列主元消去法，并且可以在程序中体现出其优势。

2 实验任务

3 实验相关知识

原始高斯消去原理：
       将方程组换成增广矩阵，利用初等行变化把增高矩阵转换成行阶梯阵，然后回代求出方程的解。
列主元加入的原因：
       列主元消去法是为控制舍入误差而提出来的，在高斯消去法的消元过程中,若出现 a=0,则消元无法进行,即使其不为 0,但很小,把它作为除数,就会导致其他元素量级的巨大增长和舍入误差的扩散,最后使计算结果不可靠.使用列主元消去法计算,基本上能控制舍入误差的影响。
列主元加入的方法：
       对给定的 n 阶线性方程组 Ax=b，首先进行列主元消去过程，然后进行回代过程，最后得到解或确定该线性方程组是奇异的。
       其具体方法是：在进行第 k(k=1,2,…,n-1)步消元时，从第 k 列的 akk 及其以下的各元素中选取绝对值最大的元素，然后通过行变换将它交换到主元素 akk 的位置上，再进行消元。

4 实验流程图

原始高斯消去法：

列主元消去法：

5 实验代码

原始高斯消去法：

import numpy as np
def gauss(a, b):
    cout = 0           #定义计算次数
    m, n = a.shape   #矩阵a的行数和列数
    if ( m < n ):
        print("There is a 解空间。")#保证方程个数大于未知数个数
    else:
        l = np.zeros((n,n))
        for i in range(n):
            # 限制条件
            if (a[i][i] == 0):
                print("no answer")


        # j表示列
        for k in range(n - 1):          # k表示第一层循环，(0，n-1)行
            for i in range(k + 1, n):   # i表示第二层循环,(k+1,n)行,计算该行消元的系数
                l[i][k] = a[i][k] / a[k][k]     #计算l
                cout += 1               #计算次数加一
                for j in range(m):      # j表示列，对每一列进行运算
                    a[i][j] = a[i][j] - l[i][k] * a[k][j]
                    cout += 1
                b[i] = b[i] - l[i][k] * b[k]
        # 回代求出方程解
        x = np.zeros(n)                       #先将解赋值为零，再一一计算
        x[n - 1] = b[n - 1] / a[n - 1][n - 1] #先算最后一位的x解

        for i in range(n - 2, -1, -1):      #依次回代倒着算每一个解
            for j in range(i + 1, n):
                b[i] -= a[i][j] * x[j]       #自增自减
            x[i] = b[i] / a[i][i]
        for i in range(n):
            print("x" + str(i + 1) + " = ", x[i])
        print("x" " = ", x)
        print("计算次数", "=", cout)


if __name__ == '__main__':
    a = np.array([[0.9428, 0.3475, -0.8468], [0.3475, 1.8423, 0.4759], [-0.8468, 0.4759, 1.2147]])  # 输入的系数矩阵
    b = np.array([0.4127, 1.7321, -0.8621])                                     #增广的一列矩阵
    gauss(a, b)                                                          #进行函数guassin运算

列主元消去法：

import numpy as np
def swap(a, b, k, n):           # 找到主元并交换，这仅是一个仅用来交换的函数
    ans = 0
    for i in range(k, n):
        if ans < np.fabs(a[i][k]):    #fabs是绝对值，将a中绝对值最大的找出来
            ans = a[i][k]
            maxn = i
    a[[k, maxn], :] = a[[maxn, k], :]     #交换
    b[k], b[maxn] = b[maxn], b[k]

#主算法
def gaussin(a, b):
    cout = 0                         #定义计算次数
    m, n = a.shape                  #矩阵a的行数和列数
    if ( m < n ):
        print("There is a 解空间。")#保证方程个数大于未知数个数
    else:
        l = np.zeros((n,n))
        for i in range(n):
            # 限制条件
            if (a[i][i] == 0):
                print("no answer")
        # j表示列
        for k in range(n - 1):          # k表示第一层循环，(0，n-1)行
            swap(a, b, k, n)            #在每次计算前，找到最大主元，进行换行
            for i in range(k + 1, n):   # i表示第二层循环,(k+1,n)行,计算该行消元的系数
                l[i][k] = a[i][k] / a[k][k]     #计算
                cout += 1
                for j in range(m):      # j表示列，对每一列进行运算
                    a[i][j] = a[i][j] - l[i][k] * a[k][j]
                    cout += 1
                b[i] = b[i] - l[i][k] * b[k]
        # 回代求出方程解
        x = np.zeros(n)
        x[n - 1] = b[n - 1] / a[n - 1][n - 1] #先算最后一位的x解

        for i in range(n - 2, -1, -1):      #依次回代倒着算每一个解
            for j in range(i + 1, n):
                b[i] -= a[i][j] * x[j]
                               #自增自减
            x[i] = b[i] / a[i][i]
        for i in range(n):
            print("x" + str(i + 1) + " = ", x[i])
        print("x" " = ", x)
        print("计算次数","=",cout)


if __name__ == '__main__':      #当模块被直接运行时，以下代码块将被运行，当模块是被导入时，代码块不被运行。
    a = np.array([[0.9428, 0.3475, -0.8468], [0.3475, 1.8423, 0.4759], [-0.8468, 0.4759, 1.2147]])  # 输入的系数矩阵
    b = np.array([0.4127, 1.7321, -0.8621])
    gaussin(a, b)

6 实验结果

原始高斯消去法：

列主元消去法：

7 总结

       线性方程组有很多种解法，可以最简单的直接代入消元计算，但是运算量较大，且过程复杂不直观。而高斯消元法，不仅可用来求解线性方程组，并可以求出矩阵的秩，以及求出可逆方阵的逆矩阵；高斯消元法目的是预处理方程组的系数矩阵，将系数矩阵变换为上三角矩阵，这样整个方程就变得清晰直观很多。
       而列主元高斯消元法，是为了克服高斯消元法的”除法 bug“的改进版本，因为消元过程中用到了除法，对于计算机来说，如果一个很小的值作为分母，则其在迭代中产生的误差，在做除法后会被放大很多倍。如 0.0001 与 0.0002 的误差小，但是1/0.0001 与 1/0.0002 的误差巨大。列主元思想，就是顺应而生，在消元前，把列中最大的元素所在行放到矩阵的最上方，那么就可以让较大值做分母，从而避免了误差被除法放大的问题。所以高斯列主元消去法在求解方程组时，使得求解得到的结果更加的准确。