数值分析——直接法解线性方程组

  数值分析中，直接法解线性方程组，python编写。具体代码见我的下载。

'''
    解线性方程组
    =====================================================
    =   直接法：（A为方阵）
    =       1. 高斯消去法/列主元消去法
    =       2. 三角分解法（A为方阵）
    =             直接分解法: A = L(DU)
    =             平方根法:   A = LDL.T
    =             追赶法:     A = (LD)U
    =====================================================
    =   创建日期： 2018/10/29 20:52
    =   修改日期： 2019/01/05 16:55
'''

# 引入numpy库
import numpy as np

# 测试矩阵

A1 = np.array([[ 5,  2,  1],
               [-1,  4,  2],
               [ 2, -3,  10]])
b1 = np.array([[-12],
               [ 20],
               [ 3]])

A2 = np.array([[ 0.001,  2.000,  3.000],
               [-1.000,  3.712,  4.622],
               [-2.000, 1.072,  5.643]])
b2 = np.array([[1.000],
               [2.000],
               [3.000]])

A3 = np.array([[1,2,3,4],\
               [1,2,3,4],\
               [1,2,3,4],\
               [1,2,3,4]])
b3 = np.array([[1],\
               [1],\
               [1],\
               [1]])

def gauss_eli(A, b):
    '''
    高斯消去法
    :param A:     系数矩阵
    :param b:
    :return:      解向量X
    '''

    # 将A和b组合成增广矩阵
    m,n = A.shape[0], A.shape[1]
    amat = np.zeros((m,n+1))
    amat[:, 0:n] = A
    amat[:, n] = b[:,0]
    # 获取待求解增广矩阵的系数矩阵的行数和列数
    (rmat_row, rmat_col) = np.shape(amat)  # rmat -- 系数矩阵(ratio matrix)
    rmat_col -= 1

    # 消元过程
    for k in range(rmat_row-1):  # n-1次消元
        # 计算Mik
        m_ik = np.array([0.0]*(rmat_row-k-1))  # 存放m_ik
        if amat[k, k] != 0:
            for c in range(rmat_row-k-1):
                m_ik[c] = amat[k+1+c, k] / amat[k, k]  # 公式
            # 计算消元K次后各元素的值
            for i in range(rmat_row-k-1):
                for j in range(rmat_col+1):  
                    amat[i+1+k, j] = amat[i+1+k, j] - m_ik[i]*amat[k, j]  # 公式
        else:
            print('主对角线元素为0！')
            break
        # 输出A
        print('第',k+1,'次：')
        print(amat)

    # 回代过程 -- 求解解向量
    x_arr = np.array([0.0]*rmat_row)  # 存放解向量
    x_arr[rmat_row-1] = amat[rmat_row-1, rmat_col] / amat[rmat_row-1, rmat_col-1]  # Xn
    for r in range(rmat_row-1):
        sum = 0.0
        for c in range(rmat_col-1-(rmat_row-r-2)):
            sum = sum + amat[rmat_row-r-2, (rmat_row-r-2)+c+1]*x_arr[rmat_row-r-1+c]
        x_arr[rmat_row-r-2] = (amat[rmat_row-r-2, rmat_col] - sum) / amat[rmat_row-r-2, rmat_row-r-2]
    # 返回解向量
    return x_arr.reshape(rmat_row, 1)

def gauss_col(A, b):
    '''
    :列主消元法
    :param A:
    :param b:
    :return:
    '''

    # 将A和b组合成增广矩阵
    m, n = A.shape[0], A.shape[1]
    amat = np.zeros((m, n + 1))
    amat[:, 0:n] = A
    amat[:, n] = b[:, 0]
    # 获取待求解增广矩阵的系数矩阵的行数和列数
    (rmat_row, rmat_col) = np.shape(amat)  # rmat -- 系数矩阵(ratio matrix)
    rmat_col -= 1
    # 消元过程
    for k in range(rmat_row-1):  # n-1次消元
        # 是否换行
        _max_index = np.abs(amat[:, k]).argmax()
        temp = np.array([0.0]*(rmat_col+1))
        for i in range((rmat_col+1)):
            temp[i] = amat[_max_index, i]
        amat[_max_index,:] = amat[k,:]
        for i in range((rmat_col+1)):
            amat[k,i] = temp[i]
        print('第',k+1,'换行：')
        print(amat)
        # 计算Mik
        m_ik = np.array([0.0]*(rmat_row-k-1))  # 存放m_ik
        if amat[k, k] != 0:
            for c in range(rmat_row - k - 1):
                m_ik[c] = amat[k+1+c, k] / amat[k, k]  # 公式
            # 计算消元K次后各元素的值
            for i in range(rmat_row-k-1):
                for j in range(rmat_col + 1):
                    amat[i+1+k, j] = amat[i+1+k, j] - m_ik[i]*amat[k,j]  # 公式
        else:
            print('主对角线元素为0！')
            break
        # 输出A
        print('第', k + 1, '次消元：')
        print(amat)

    # 回代过程 -- 求解解向量
    x_arr = np.array([0.0] * rmat_row)  # 存放解向量
    x_arr[rmat_row-1] = amat[rmat_row-1, rmat_col] / amat[rmat_row-1, rmat_col-1]  # Xn
    for r in range(rmat_row-1):
        sum = 0.0
        for c in range(rmat_col-1-(rmat_row-r-2)):
            sum = sum + amat[rmat_row-r-2, (rmat_row-r-2)+c+1] * x_arr[rmat_row-r-1+c]
        x_arr[rmat_row-r-2] = (amat[rmat_row-r-2, rmat_col]-sum) / amat[rmat_row-r-2, rmat_row-r-2]
    # 返回解向量
    return x_arr.reshape(rmat_row, 1)

def LU(A, b):
    '''
    :直接三角形分解法 -- 使用的是杜利特尔分解
    :param A:   系数矩阵
    :param b:
    :return:    解向量X、L、U
    '''

    m = A.shape[0]   # 行数
    n = A.shape[1]   # 列数

    A = np.array(A, dtype=float)  # 如果A的所有元素为整数，必须变换A的类型！！！

    # 矩阵LU分解
    L = np.eye(m)   # 单位下三角矩阵
    for k in range(m-1):
        # 计算Mik
        m_ik = np.array([0.0]*(m-1-k))  # 存放m_ik
        if A[k, k] != 0:
            for c in range(len(m_ik)):
                m_ik[c] = A[k+1+c, k] / A[k, k]  # 公式
                A[c+1+k,:] = A[c+1+k,:] + (-m_ik[c]*A[k,:])
        for i in range(m-k-1):
            L[i+1+k, k] = m_ik[i]    # L矩阵
    U = A  # 上三角矩阵

    # Ly = b
    y = np.dot(np.linalg.inv(L), b)
    # Ux = y 等价于 (DU)x = y
    x = np.dot(np.linalg.inv(U), y)
    # 返回x, L, U
    return x, L, U

def purchase(A, b):
    '''
    :追赶法 -- 使用的是crout分解
    :param A:  系数矩阵
    :param b:
    :return:   返回x、L、D、U
    '''

    m = A.shape[0]  # 行数
    n = A.shape[1]  # 列数

    A = np.array(A, dtype=float)

    # 矩阵LU分解
    L = np.eye(m)  # 单位下三角矩阵
    for k in range(m - 1):
        # 计算Mik
        m_ik = np.array([0.0] * (m - 1 - k))  # 存放m_ik
        if A[k, k] != 0:
            for c in range(len(m_ik)):
                m_ik[c] = A[k + 1 + c, k] / A[k, k]  # 公式
                A[c + 1 + k, :] = A[c + 1 + k, :] + (-m_ik[c] * A[k, :])
        for i in range(m - k - 1):
            L[i + 1 + k, k] = m_ik[i]  # L矩阵
    U = A  # 上三角矩阵
    D = np.zeros((m,n))  # 对角矩阵
    for i in range(m):
        if U[i,i] != 0:
            D[i,i] = U[i,i]
            U[i,:] = U[i,:] / D[i,i]  # 把U变为单位下三角矩阵

    # (LD)y = b
    y = np.dot(np.linalg.inv(np.dot(L,D)), b)
    # Ux = y
    x = np.dot(np.linalg.inv(U), y)
    # 返回x, L, U
    return x, L, D, U

def sqrt2(A, b):
    '''
    :改进的平方根法
    :param A:  A为正定实对称矩阵
    :param b:
    :return:
    '''

    m = A.shape[0]  # 行数
    n = A.shape[1]  # 列数

    A = np.array(A, dtype=float)

    # 矩阵LU分解
    L = np.eye(m)  # 单位下三角矩阵
    for k in range(m - 1):
        # 计算Mik
        m_ik = np.array([0.0] * (m - 1 - k))  # 存放m_ik
        if A[k, k] != 0:
            for c in range(len(m_ik)):
                m_ik[c] = A[k + 1 + c, k] / A[k, k]  # 公式
                A[c + 1 + k, :] = A[c + 1 + k, :] + (-m_ik[c] * A[k, :])
        for i in range(m - k - 1):
            L[i + 1 + k, k] = m_ik[i]  # L矩阵
    # 在该方法中，L.T = U
    U = A  # 上三角矩阵
    D = np.zeros((m, n))  # 对角矩阵
    for i in range(m):
        if U[i, i] != 0:
            D[i, i] = U[i, i]
            U[i, :] = U[i, :] / D[i, i]  # 把U变为单位下三角矩阵

    # Ly = b
    y = np.dot(np.linalg.inv(np.dot(L, D)), b)
    # (DL.T)x = y
    x = np.dot(np.linalg.inv(U), y)
    # 返回x, L, U
    return x, L, D, U

def main():

    print('=' * 25, '高斯消去法', '=' * 25)
    x = gauss_eli(A2, b2)
    print('解：')
    print(x)

    print('=' * 25, '列主消去法', '=' * 25)
    x1 = gauss_col(A2, b2)
    print('解：')
    print(x1)

    xlu,L1,U1 = LU(A2, b2)
    print('='*25,'直接三角形分解法','='*25)
    print('单位下三角矩阵-L：')
    print(L1)
    print('上三角矩阵-U：')
    print(U1)
    print('解：')
    print(xlu)

    xpur, L2, D2, U2 = purchase(A2, b2)
    print('=' * 25, '追赶法', '=' * 25)
    print('单位下三角矩阵-L：')
    print(L2)
    print('对角矩阵-D：')
    print(D2)
    print('单位上三角矩阵-U：')
    print(U2)
    print('解：')
    print(xpur)
    print('=' * 25, '平方根法', '=' * 25)
    print('对于平方根法，由于A为正定实堆成矩阵，所以求的L的转置=U')
    print('只要将追赶法中的U写成L的转置即可，其余的同追赶法！')

    input('')

if __name__ == "__main__":
    main()