数值分析之高斯消去法解线性方程组

@数值分析之非线性方程AX=B求解

一、高斯消去法（偏序选主元策略）

处理过程中，不可能每次都是三角矩阵。所以我们要采用新的方法来处理一般方程组，将其转化为上三角形式，再采用回代算法求解。这里我们介绍高斯消去法。

高斯消去法之所以可以对方程组进行一定的变换，因为线性代数中定理可知：对线性方程组进行3种初等变换可以将AX=B变换为另一个等价的线性方程组UX=Y（U是上三角矩阵），这段过程就是高斯消去。

同理，对矩阵也可等价变换。而编程中第一步就是让AB合并为一个增广矩阵。让[A|B]一起被变换。

1.1何为高斯消去法？

高斯消去法：对AX=B的增广矩阵[A|B]进行初等行变换，使之变成UX=Y的形式。A变成了U（上三角矩阵）。

其中系数矩阵A中的元素A[r,r]用来消去A[k,r]，k=r+1,r+2…,N，这里称A[r,r]为主元，第r行称为主元行。
下面给出一个例子。高斯消去过程其实就是我们手算线性方程组的过程。

1.2 高斯消去法的问题

但是单纯的高斯消去法，并不普遍适用。而下面一种情况会导致处理过程失败。

若A[k,k]=0，则不能使用第k行来消除第k列的元素，而需要将第k行与对角线下的某行进行交换，来得到一个非0主元。若没有非0主元，则系数矩阵是奇异的，因此不存在方程组唯一解。

这需要采用选主元的策略来规避这种情况。

1.3偏序选主元策略的意义

问题1：选主元策略可以避免A[k,k]=0，但是又有问题出现。==一般情况下，供选择的行很多，应该选哪一行？

答：选择导致误差更小的那一行。
由于计算机使用固定精度计算，每次计算中，都可能出现微小的误差，在指数级的计算后，误差传播，导致计算值偏离实际值。所以在选择主元时，要主动选择导致计算误差更小的那一行。
这就是偏序选主元的作用。

问题2：何为偏序选主元？

答：就是将某一列元素（不算对角线以上）中的最大值与主元行交换。在上面的例子3.16中，最大的误差出现在求解m的过程中，m是某一列中，对角线下的元素除以主元行得到的倍数。m越小，经过指数次计算后，误差传播就越小。

1.4 matlab版算法：

最后一行的backsub就是上篇博客中的上三角回代算法

1.5 病态情况

当存在矩阵B，当矩阵B和A中系数元素微小的变化使得X=A^(-1)B变化很大，则矩阵A是病态矩阵。

当A近似于奇异而且行列式接近0时，就可能发生病态情况。病态情况可能会导致出现错误解。而这类灵敏性分析是高级数值分析的领域，本人并未有了解。

二、题目及实现代码

2.1 题目

为求解一个线性方程组，首先构造增广矩阵[A|B]，
采用偏序选主元策略的高斯消去法变换成上三角矩阵，再执行回代过程得到解。

2.2 输入输出格式

【输入形式】在屏幕上依次输入方阵阶数n，系数矩阵A和常数矩阵B。

【输出形式】首先输出上三角矩阵（变换后的增广矩阵），然后每一行输出一个根

2.3 样例

输入

4
1 2 1 4
2 0 4 3
4 2 2 1
-3 1 3 2
13
28
20
6

输出

[[ 4. 2. 2. 1. 20. ]
[ 0. 2.5 4.5 2.75 21. ]
[ 0. 0. 4.8 3.6 26.4 ]
[ 0. 0. 0. 3.75 7.5 ]]
[[ 3.]
[-1.]
[ 4.]
[ 2.]]

2.4 思路和要点

处理过程：
1、偏序选主元（使得参与运算的列中，最大值所在一行与其中最高行交换，来增加精度）
2、让[A|B]的增广矩阵变成上三角矩阵
3、回代求解上三角矩阵，得X

2.5 代码

import numpy as np
def uptrbk(Aug,n): #偏序选主元，构建上三角矩阵
    for i in range(n):
        C = Aug[i:,i]
        #偏序选主元
        row = np.argmax(abs(C))
        #找出该列最大值，并交换，保证a(i,i)!=0
        Aug[[i, row+i ], :] = Aug[[i+row, i], :]
        for k in range(i+1,n):
            m = Aug[k,i]/Aug[i,i]
            Aug[k,:] = Aug[k,:]-m*Aug[i,:]
    print(Aug)
    return Aug

def backsub(A, B, X): #上三角回代求解
    n = len(X)
    X[n - 1] = B[n - 1] / A[n - 1, n - 1]
    for i in range(n - 2, -1, -1):
        temp = np.reshape(A[i, i + 1:], (1, len(X[i + 1:]))).dot(X[i + 1:])
        # 截断的A虽然看似可以和X想点乘，但是程序会判断为形状不同，所以要重新赋予形状
        X[i] = (B[i] - temp) / A[i, i]
    return X

def main():
    n = int(input())
    A = np.zeros([n, n], dtype=np.double)
    B = np.zeros([n, 1], dtype=np.double)
    X = np.zeros([n, 1], dtype=np.double)
    for i in range(n):
        A[i, :] = np.array(input().split(), dtype=np.double)
    for i in range(n):
        B[i] = np.array(input(), dtype=np.double)
    Aug = np.append(A, B, axis=1)#将A，B合并
    uptrbk(Aug,n)
    backsub(Aug[0:n, 0:n], Aug[0:n, n], X)
    print(X)

if __name__ == '__main__':
    main()

2.6 结果

样例：