题目

实现高斯消去、列主元高斯消去，LU分解分别求解线性方程组

方法一：高斯消去法

消去过程：
$$l_{ik}=a_{ik}^{(k-1)}/a_{kk}^{(k-1)}(i=k+1,k+2,...,n)$$
$$a_{ij}^{(k)}=a_{ij}^{(k-1)}-l_{ik}{a}_{jk}^{(k-1)}(i,j=k+1,k+2,...,n)$$
$$b_i^{(k)}=b_i^{(k-1)}-l_{ik}{b}_k^{(k-1)}(i=k+1,k+2,...,n)$$
回代过程：
$$x_n=b_n^{(n-1)}/a_{nn}^{(n-1)}$$
$$x_k=(b_k^{(k-1)}-\sum _{j=k+1}^n{a}_{kj}^{(k-1)}{x}_j)/a_{kk}^{(k-1)}$$
在消元法第k次迭代时，首先 ，寻找第k列中从k到n处的绝对值最大的一个，然后互换该行和第k行的所有位置。 然后再消元。这样的好处是能控制舍入误差的扩大和传播。 基本算法如下：

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]])
def guass(A,b):
    n = b.shape[0]
    A=A.astype('float')
    b=b.astype('float')
    #消元
    for k in range(n-1):    #k第一层循环，第（0~n-1）行
#         print('第',k+1,'次的初值')
#         print(A)
#         print(b)
        for i in range(k+1,n):    #i第二层循环，第（k+1,n）行
            mik = A[i,k]/A[k,k]
            for j in range(k+1,n):    
                A[i,j] = A[i,j] - mik*A[k,j]
            b[i] = b[i] - mik*b[k]
#         print('第',k+1,'次消元后的值')
#         print(A)
#         print(b)
    #回代
    solution = np.zeros(n)
    solution[n-1] = b[n-1]/A[n-1,n-1]    #改为n-1
    for i in range(n-2,-1,-1):
        for j in range(i+1,n):
            solution[i] = solution[i] + A[i,j]*solution[j]
        solution[i] = (b[i] - solution[i])/A[i,i]
    return solution

guass(A,b)

结果截图

方法二：列主元素高斯消元法

def guass_zhuyuan(A,b):
    n = b.shape[0]
    A=A.astype('float')
    b=b.astype('float')
    #寻找列主元素
    for k in range(n-1):     #k第一层循环，第（0~n-2）行
#         print('第',k+1,'次的初值')
#         print(A)
#         print(b)
        arr = abs(A[k:n,k])  #每一列最大值的数组
        value = max(arr)     #最大值
        position = np.argmax(arr)    #最大值索引
        
        #print(position)
        if position != 0:
#             A[k,k:n-1],A[position+k,k:n-1] = A[position+k,k:n-1],A[k,k:n-1] 
#             b[k],b[position+k] = b[position+k],b[k]
            A[[k,position+k],k:n] = A[[position+k,k],k:n]
            b[[k,position+k],:] = b[[position+k,k],:] 
#         print('第',k+1,'次交换后的值')
#         print(A)
#         print(b)
        #消元
        for i in range(k+1,n):    #i第二层循环，第（k+1,n-1）行
            mik = A[i,k]/A[k,k]
            for j in range(k+1,n):    #j从k开始循环，不是k+1
                A[i,j] -= mik*A[k,j]
            b[i] = b[i] - mik*b[k]
#         print('第',k+1,'次消元后的值')
#         print(A)
#         print(b)

    #回代
    solution = np.zeros(n)
    solution[n-1] = b[n-1]/A[n-1,n-1]    #改为n-1
    for i in range(n-2,-1,-1):
        for j in range(i+1,n):
            solution[i] = solution[i] + A[i,j]*solution[j]
        solution[i] = (b[i] - solution[i])/A[i,i]
    return solution                     
        

结果截图

方法三：LU分解

def LUfenjie(A,b):
    n = b.shape[0]
    A = A.astype('float')
    b = b.astype('float')
    U = np.zeros(A.shape)
    L = np.zeros(A.shape)
    U[0,:] = A[0,:]    #U的第一行元素和A的第一行元素相等
    L[1:n,0] = A[1:n,0]/U[0,0]    #L的每一行的第一列元素和U第一个元素相乘=A的对应元素
    #分解
    #U为上三角元素 k=1~n-1行 j=k~n-1列 r=0~k-2行 i=k+1~n-1行
#     for k in range(1,n):    
#         for j in range(k,n):
#             U[k,j] = A[k,j] - L[k,0:k-1] * U[0:k-1,j]
#         for i in range(k+1,n):
#             L[i,k] = (A[i,k]-L[i,0:k-1] * U[0:k-1,k])/U[k,k]
    #分解
    #U为上三角元素 k=1~n-1行 j=k~n-1列 r=0~k-2行 i=k+1~n-1行
    for k in range(1,n):    
        for j in range(k,n):
            delta = 0
            for r in range(0,k):
                delta += L[k,r] * U[r,j]
            U[k,j] = A[k,j] - delta
        for i in range(k+1,n):
            theta = 0
            for r in range(0,k):
                theta += L[i,r] * U[r,k]
            L[i,k] = (A[i,k] - theta)/U[k,k]
    #求解y
    y = np.zeros(b.shape[0])
    #忘了初始化……
    y[0] = b[0]
    for i in range(1,n):
        alpha = 0
        for j in range(0,i):
            alpha += L[i,j]*y[j].T
        y[i] = b[i] - alpha

        
    #求解x
    x = np.zeros(b.shape[0])
    x[n-1] = y[n-1]/U[n-1,n-1]
    for i in range(n-2,-1,-1):
        gama = 0 
        for j in range(i+1,n):
            gama += U[i,j]*x[j].T
        x[i] = (y[i] - gama )/U[i,i]
    return (x)           

结果截图

结果总结