python求解线性方程组（高斯、列主元高斯、直接三角）

求解上述线性方程组

高斯消去法

消元

回带

#高斯消去法
import numpy as np
A = np.array([[10,-7,0,1],[-3,2.099999,6,2],[5,-1,5,-1],[2,1,0,2]])
b = np.array([8,5.900001,5,1])
for k in range(3) :
    for i in range(k+1,4) :
        mik = A[i,k]/A[k,k]
        for j in range(k+1,4) :
            A[i,j] -= mik * A[k,j]
        b[i] = b[i] - mik * b[k]
solution = [0,0,0,0]
solution[3] = b[3] / A[3,3]
for i in range(2,-1,-1) :
    for j in range(i+1,4) :
        solution[i] = solution[i] + A[i,j] * solution[j]
    solution[i] = (b[i] - solution[i]) / A[i,i]
print('高斯消去法，最终的解X_i(角标从小到大):{}'.format(solution))

列主元高斯消去法

#列主元高斯消去法
import numpy as np
A = np.array([[10,-7,0,1],[-3,2.099999,6,2],[5,-1,5,-1],[2,1,0,2]])
b = np.array([8,5.900001,5,1])
for k in range(3) :
    value = max(abs(A[k:4,k]))  #主元的值
    position = np.argmax(A[k:4,k])  #主元所在位置（相对位置）
    if position != k :  #若a(k,k)不是绝对值最大的，换位置
        A[k,k:4],A[position+k,k:4] = A[position+k,k:4],A[k,k:4]
        b[k],b[position+k] = b[position+k],b[k]
    for i in range(k+1,4) :
        mik = A[i,k]/A[k,k]
        for j in range(k+1,4) :
            A[i,j] -= mik * A[k,j]
        b[i] = b[i] - mik * b[k]
        
#回带过程不变
solution = [0,0,0,0]
solution[3] = b[3] / A[3,3]
for i in range(2,-1,-1) :
    for j in range(i+1,4) :
        solution[i] = solution[i] + A[i,j] * solution[j]
    solution[i] = (b[i] - solution[i]) / A[i,i]
print('列主元高斯消去法，最终的解X_i(角标从小到大):{}'.format(solution))

解析

position = np.argmax(A[k:4,k])  #主元所在位置（相对位置）

当k=2时，仅比较第三行和第四行中第三列的元素绝对值大小，若第四行第三列元素绝对值较大则为主元，此时position为1，是在比较序列中的相对位置，而在矩阵中的绝对位置（下标从0开始）为position+k=3

直接三角分解法

#直接三角分解法
import numpy as np
A = np.array([[10,-7,0,1],[-3,2.099999,6,2],[5,-1,5,-1],[2,1,0,2]])
b = np.array([8,5.900001,5,1])
y = np.zeros(4)
x = np.zeros(4)
L = np.zeros((4,4))
U = np.zeros((4,4))

# 分解A = LU
U[0,:] = A[0,:]
L[1:4,0] = A[1:4,0] / U[0,0]
for k in range(1,4) :
    for j in range(k,4) :
        U[k,j] = A[k,j] - L[k,0] * U[0,j]
        for m in range(1,k) :
            U[k,j] -= L[k,m] * U[m,j]
    for i in range(k+1,4) :
        L[i,k] = A[i,k] - L[i,0] * U[0,k]
        for m in range(1,k) :
            L[i,k] -= L[i,m] * U[m,k]
        L[i,k] = L[i,k] / U[k,k]
            
# 解Ly = b
y[0] = b[0]
for i in range(1,4) :
    y[i] = b[i] - L[i,0] * y[0].T
    for n in range(1,i) :
        y[i] -= L[i,n] * y[n].T
    
# 解Ux = y
x[3] = y[3] / U[3,3]
for i in range(2,-1,-1) :
    x[i] = (y[i] - U[i,i+1] * x[i+1].T) 
    for n in range(i+2,4) :
        x[i] -= U[i,n] * x[n].T
    x[i] = x[i] / U[i,i]
print('直接三角分解法，最终的解X_i(角标从小到大):{}'.format(x))

结果列表

方法/解	高斯消去法	列主元高斯消去法	直接三角分解法
x1	-6.033918253933735e-10	-6.033918253933735e-10	-6.0339177e-10
x2	-1.0000000008881784	-1.0000000008881784	-1
x3	1.000000000070277	1.000000000070277	1
x4	0.9999999998166688	0.9999999998166688	1

结果截图