博客

考虑n阶线性方程组：

高斯消去法python

import numpy as np

def gaussin(a, b):
    m, n = a.shape
    c = np.zeros(n)
    for i in range(n):
        #限制条件
        if (a[i][i] == 0):
            print("no answer")
    #k表示第一层循环，(0，n-1)行
    #i表示第二层循环,(k+1,n)行,计算该行消元的系数
    #j表示列
    for k in range(n - 1):
        for i in range(k + 1, n):
            c[i] = a[i][k] / a[k][k]
            for j in range(m):
                a[i][j] = a[i][j] - c[i] * a[k][j]
            b[i] = b[i] - c[i] * b[k]
    x = np.zeros(n)
    x[n - 1] = b[n - 1] / a[n - 1][n - 1]
    #回代求出方程解
    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])


if __name__ == '__main__':
    a = np.array([[2, 4, -6], [1, 5, 3], [1, 3, 2]])
    b = np.array([-4.0, 10.0, 5.0])
    gaussin(a, b)
    print('*' * 20)

但是，在消元过程中，无法使主元素a(ii)≠0，但是很小时，用其做除数，会导致其他元素数量级的严重增长，舍入误差的扩展，最后导致计算结果不可靠。所以这次采用列主元素消去法来进行，思想就是将有小数的那行与该列中数最大的那行进行交换。

python代码如下：

import numpy as np
#高斯列主元消去法
#找到主元并交换
def swap(a,b,k,n):
    ans = -1.
    for i in range(k,n):
        if ans < np.fabs(a[i][k]):
            ans = a[i][k]
            maxn = i
    a[[k,maxn],:] = a[[maxn,k],:]
    b[k],b[maxn] = b[maxn],b[k]
def gaussin(a, b):
    m, n = a.shape
    c = np.zeros(n)
    for i in range(n):
        #限制条件
        if (a[i][i] == 0):
            print("no answer")
    #k表示第一层循环，(0，n-1)行
    #i表示第二层循环,(k+1,n)行,计算该行消元的系数
    #j表示列
    for k in range(n - 1):
        swap(a,b,k,n)
        for i in range(k + 1, n):
            c[i] = a[i][k] / a[k][k]
            for j in range(m):
                a[i][j] = a[i][j] - c[i] * a[k][j]
            b[i] = b[i] - c[i] * b[k]
    x = np.zeros(n)
    x[n - 1] = b[n - 1] / a[n - 1][n - 1]

    #回代求出方程解
    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])
if __name__ == '__main__':
    a = np.array([[2, 4, -6], [1, 5, 3], [1, 3, 2]])
    b = np.array([-4.0, 10.0, 5.0])
    gaussin(a, b)
    print('*' * 20)
    a = np.array([[0.01, 2, -0.5], [-1, -0.5, 2], [5, -4, 0.5]])
    b = np.array([-5.0, 5.0, 9.0])
    gaussin(a, b)
    print('*' * 20)