为顺序消元法和列主元消元法我分别写了两个函数,即gauss_se (gauss sequntial elimination) 与 gauss_cpe (gauss column-pivoting elimination) 。它们的区别只在于,在 gauss_se 中,我每次循环开始检查一下 a[i, i] 是否为 0 ,如果为 0 ,则将其下不为 0 的行与当前行交换,如果下面的对应元素也都为 0 ,那就跳出。而在 gauss_cpe 里,每次循环开始不会检查,而是直接将 a[i, i] 下绝对值最大的那行与当前行交换,当然,如果这个绝对值最大也才是 0 ,那么就跳出。

import numpy as np
from numpy import dot

def gauss_se(a,b):
    a = np.copy(a)
    b = np.copy(b)    
    n = len(b)
    for i in range(n-1):
        if a[i, i] == 0:
            ind = aux_func(a[i:, i+1])
            if ind is None:
                break
            a[i,:], a[i+ind+1, :] = a[i+ind+1, :], a[i,:]
            b[i], b[i+ind+1] = b[i+ind+1], b[i]
        for j in range(i+1,n):
            if a[j,i] != 0:
                lam = a[j,i] / a[i,i]
                a[j,:] = a[j,:] - lam*a[i,:]
                b[j] = b[j] - lam*b[i]
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - dot(a[k,(k+1):n],b[(k+1):n]))/a[k,k]
    result = b
    return result

def gauss_cpe(a, b):
    a = np.copy(a)
    b = np.copy(b)    
    n = len(b)
    for i in range(n-1):
        ind = np.argmax(np.abs(a[i:, i]))
        if a[i, i+ind]==0:
            break
        a[[i, i+ind],:]= a[[i+ind, i],:]
        b[[i, i+ind]] = b[[i+ind, i]]
        for j in range(i+1,n):
            if a[j,i] != 0:
                lam = a[j,i] / a[i,i]
                a[j,:] = a[j,:] - lam*a[i,:]
                b[j] = b[j] - lam*b[i]
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - dot(a[k,(k+1):n],b[(k+1):n]))/a[k,k]
    result = b
    return result

def aux_func_1(x):
    for i,t in enumerate(x):
        if t!=0:
            return i
    return None


def main():
    A = [[1.1348, 3.8326, 1.1651, 3.4017], [0.5301, 1.7875, 2.5330, 1.5435],
         [3.4129, 4.9317, 8.7643, 1.3142], [1.2371, 4.9998, 10.6721, 0.0147]]
    b = [9.5342, 6.3941, 18.4231, 16.9237]
    A = np.array(A)
    b = np.array(b).T
    print(gauss_se(A, b))
    print(gauss_cpe(A, b))

if __name__ == '__main__':
    main()

输出如下:

[ 1.  1.  1.  1.]
[ 1.  1.  1.  1.]