列主元消去法和LU矩阵分解以及矩阵迭代求解

最新推荐文章于 2022-10-10 20:00:37 发布

仁，上古

最新推荐文章于 2022-10-10 20:00:37 发布

阅读量1k

点赞数 15

文章标签：线性代数抽象代数

本文链接：https://blog.csdn.net/qq_45458674/article/details/117625919

版权

上三角矩阵求解

为了方便后文求解矩阵，先将上三角矩阵、和下三角矩阵求解打包，这里输入上三角矩阵A和解矩阵b。注意这里的矩阵是A，b都是二维数组。例如：

A =np.array([3.01,6.03,1.99,1.27,4.16,-1.23,0.987,-4.81,9.34],dtype=float).reshape(3,3)
B = np.array([1,1,1], dtype=float).reshape(3, 1)

#上三角求解：
def up_solve(A,b):
    n = len(A[0])
    A[n-1,n-1] = b[n-1][0]/A[n-1,n-1]
    for i in range(1,n):
        k = n-i-1#取逆序
        A[k,k] = (b[k][0] - sum(A[k]*A[k+1]))/A[k,k]
        A[k]=np.hstack((A[k,:k+1],A[k+1,k+1:]))#拼接新数组
    # print('方程解为',A[0])
    return A[0]

下三角矩阵求解

同上文一样，这里下三角矩阵求解：

#下三角求解
def down_solve(A,b):
    n = len(A[0])
    A[0,0] = b[0][0] /A[0,0]
    for i in range(1, n):
        k = i  # 取逆序
        A[k, k] = (b[k][0] - sum(A[k] * A[k-1])) / A[k, k]
        A[k] = np.hstack((A[k-1, :k], A[k, k:]))  # 拼接新数组
    # print('方程解为', A[-1])
    return A[-1]

高斯消去法

这里主要是用到的列主元消去法，为了更直观的观察消去过程，在代码中增加了消去过程输出。消元后得到上三角矩阵，然后调用上面上三角矩阵求解以及输出了系数矩阵的行列式。

#列主元高斯消去
def gaussian(A,b):
    #最大值索引
    n = len(A[0])
    A = np.hstack([A,b])
    change_index=0
    for i in range(n):
        tem_lie = abs(A[:,i])#列主元最大值
        index = np.where(tem_lie==np.max(tem_lie[i:]))[0][0]#注意消去后的新索引
        if index !=i:
            change_index +=1
            print('交换次序%s行与%s行：'%(i,index))
            A[i,:],A[index,:]=A[index,:].copy(),A[i,:].copy()
            print(A)
        #消去部分：
        for  j in range(i+1,n):
            tem = A[i].copy()/A[i,i]
            A[j]=tem*A[j,i]-A[j]
        print('消去结果：\n',A)
    #回带求解
    A,b = A[:,:n],A[:,-1]
    # print(A,b)
    print('高斯消去后行列式|A|', np.linalg.det(A) * ((-1) ** change_index))
    x = up_solve(A,b.reshape(n,1))
    print('求解x为：',x)
    return x

LU矩阵分解及其求解

这里借鉴过其他博主的思路，先求得分解的的LU矩阵，再利用下三角求解获取y,再通过下三角求解获取X

#Lu及其求解
def LU_decomposition(A):
    n = len(A[0])
    L = np.zeros([n, n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i] = 1
        if i == 0:
            U[0][0] = A[0][0]
            for j in range(1, n):
                U[0][j] = A[0][j]
                L[j][0] = A[j][0] / U[0][0]
        else:
            for j in range(i, n):  # U
                temp = 0
                for k in range(0, i):
                    temp = temp + L[i][k] * U[k][j]
                U[i][j] = A[i][j] - temp
            for j in range(i + 1, n):  # L
                temp = 0
                for k in range(0, i):
                    temp = temp + L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - temp) / U[i][i]
    return L, U
# 计算线性方程组的解
def answer(matrix_l,matrix_u,matrix_an):#LU以及B矩阵

    y = down_solve(matrix_l,matrix_an)
    print('y的值为：',y)

    x = up_solve(matrix_u,y.reshape(len(y),1))
    print('x的解为：',x)
    return x

雅克比迭代法

通过对谱半径求解，判断改迭代矩阵是否发散，再利用范数判断是否退出迭代：输出迭代次数，以及迭代的解值，初值全部默认设为0，误差默认为1e-6有需要可以自己设置初值，与上文一样，A,b都是二维数组

#Jacobian雅克比迭代
def jacobian(A,b,x0=None,mistake=1e-6):
    if x0 ==None:
        x0 = [[0] for i in range(len(A[0]))]
    L = (-1)*np.tril(A,-1)
    D =np.diag(np.diag(A))
    U = (-1)*np.triu(A,1)
    B = np.dot(np.linalg.inv(D),(L+U))
    b = np.dot(np.linalg.inv(D),b)
    t_a,t_b = np.linalg.eig(B)
    if max(abs(t_a)) >1:
        print('改矩阵发散')
        return
    k = 1
    while True:
        x1 = np.dot(B,x0)+b
        print('第%s次迭代'%k,x1)
        norm = max(abs(x1-x0))
        if norm <mistake:
            break
        x0 = x1
        k+=1

高斯塞得迭代法

与上文类似，更改迭代矩阵，如下：

#高斯迭代
def gauss_seidel(A,b,x0=None,mistake=1e-6):
    if x0 == None:
        x0 = [[0] for i in range(len(A[0]))]
    L = (-1) * np.tril(A, -1)
    D = np.diag(np.diag(A))
    U = (-1) * np.triu(A, 1)
    B = np.dot(np.linalg.inv(D-L),U)
    b = np.dot(np.linalg.inv(D-L), b)
    t_a, t_b = np.linalg.eig(B)
    if max(abs(t_a)) > 1:
        print('改矩阵发散')
        return
    k = 1
    while True:
        x1 = np.dot(B, x0) + b
        print('第%s次迭代' % k, x1)
        norm = max(abs(x1 - x0))
        # print(norm)#这个莫可以用来判断下降速度
        if norm < mistake:
            break
        x0 = x1
        k += 1

松弛法

这里多出了一个松弛因子w,w=1时就是高斯塞得迭代，这里w为必须参数

def sor(A,b,w,x0=None,mistake=1e-6):
    if x0 == None:
        x0 = [[0] for i in range(len(A[0]))]
    L = (-1) * np.tril(A, -1)
    D = np.diag(np.diag(A))
    U = (-1) * np.triu(A, 1)
    B = np.dot(np.linalg.inv(D-w*L),(1-w)*D+w*U)
    b = np.dot(np.linalg.inv(D-w*L), b)*w
    t_a, t_b = np.linalg.eig(B)
    if max(abs(t_a)) > 1:
        print('改矩阵发散')
        return

    k = 1
    while True:
        x1 = np.dot(B, x0) + b
        print('第%s次迭代' % k, x1)
        norm = max(abs(x1 - x0))
        # print(norm)
        if norm < mistake:
            break
        x0 = x1
        k += 1

源码

这里给出上面全部源码以及使用案列

import numpy as np

#上三角求解：
def up_solve(A,b):
    n = len(A[0])
    A[n-1,n-1] = b[n-1][0]/A[n-1,n-1]
    for i in range(1,n):
        k = n-i-1#取逆序
        A[k,k] = (b[k][0] - sum(A[k]*A[k+1]))/A[k,k]
        A[k]=np.hstack((A[k,:k+1],A[k+1,k+1:]))#拼接新数组
    # print('方程解为',A[0])
    return A[0]

#下三角求解
def down_solve(A,b):
    n = len(A[0])
    A[0,0] = b[0][0] /A[0,0]
    for i in range(1, n):
        k = i  # 取逆序
        A[k, k] = (b[k][0] - sum(A[k] * A[k-1])) / A[k, k]
        A[k] = np.hstack((A[k-1, :k], A[k, k:]))  # 拼接新数组
    # print('方程解为', A[-1])
    return A[-1]

#列主元高斯消去
def gaussian(A,b):
    #最大值索引
    n = len(A[0])
    A = np.hstack([A,b])
    change_index=0
    for i in range(n):
        tem_lie = abs(A[:,i])#列主元最大值
        index = np.where(tem_lie==np.max(tem_lie[i:]))[0][0]#注意消去后的新索引
        if index !=i:
            change_index +=1
            print('交换次序%s行与%s行：'%(i,index))
            A[i,:],A[index,:]=A[index,:].copy(),A[i,:].copy()
            print(A)
        #消去部分：
        for  j in range(i+1,n):
            tem = A[i].copy()/A[i,i]
            A[j]=tem*A[j,i]-A[j]
        print('消去结果：\n',A)
    #回带求解
    A,b = A[:,:n],A[:,-1]
    # print(A,b)
    print('高斯消去后行列式|A|', np.linalg.det(A) * ((-1) ** change_index))
    x = up_solve(A,b.reshape(n,1))
    print('求解x为：',x)
    return x


#Lu及其求解
def LU_decomposition(A):
    n = len(A[0])
    L = np.zeros([n, n])
    U = np.zeros([n, n])
    for i in range(n):
        L[i][i] = 1
        if i == 0:
            U[0][0] = A[0][0]
            for j in range(1, n):
                U[0][j] = A[0][j]
                L[j][0] = A[j][0] / U[0][0]
        else:
            for j in range(i, n):  # U
                temp = 0
                for k in range(0, i):
                    temp = temp + L[i][k] * U[k][j]
                U[i][j] = A[i][j] - temp
            for j in range(i + 1, n):  # L
                temp = 0
                for k in range(0, i):
                    temp = temp + L[j][k] * U[k][i]
                L[j][i] = (A[j][i] - temp) / U[i][i]
    return L, U
# 计算线性方程组的解
def answer(matrix_l,matrix_u,matrix_an):#LU以及B矩阵

    y = down_solve(matrix_l,matrix_an)
    print('y的值为：',y)

    x = up_solve(matrix_u,y.reshape(len(y),1))
    print('x的解为：',x)
    return x

#Jacobian雅克比迭代
def jacobian(A,b,x0=None,mistake=1e-6):
    if x0 ==None:
        x0 = [[0] for i in range(len(A[0]))]
    L = (-1)*np.tril(A,-1)
    D =np.diag(np.diag(A))
    U = (-1)*np.triu(A,1)
    B = np.dot(np.linalg.inv(D),(L+U))
    b = np.dot(np.linalg.inv(D),b)
    t_a,t_b = np.linalg.eig(B)
    if max(abs(t_a)) >1:
        print('改矩阵发散')
        return
    k = 1
    while True:
        x1 = np.dot(B,x0)+b
        print('第%s次迭代'%k,x1)
        norm = max(abs(x1-x0))
        if norm <mistake:
            break
        x0 = x1
        k+=1
#高斯迭代
def gauss_seidel(A,b,x0=None,mistake=1e-6):
    if x0 == None:
        x0 = [[0] for i in range(len(A[0]))]
    L = (-1) * np.tril(A, -1)
    D = np.diag(np.diag(A))
    U = (-1) * np.triu(A, 1)
    B = np.dot(np.linalg.inv(D-L),U)
    b = np.dot(np.linalg.inv(D-L), b)
    t_a, t_b = np.linalg.eig(B)
    if max(abs(t_a)) > 1:
        print('改矩阵发散')
        return
    k = 1
    while True:
        x1 = np.dot(B, x0) + b
        print('第%s次迭代' % k, x1)
        norm = max(abs(x1 - x0))
        # print(norm)#这个莫可以用来判断下降速度
        if norm < mistake:
            break
        x0 = x1
        k += 1
def sor(A,b,w,x0=None,mistake=1e-6):
    if x0 == None:
        x0 = [[0] for i in range(len(A[0]))]
    L = (-1) * np.tril(A, -1)
    D = np.diag(np.diag(A))
    U = (-1) * np.triu(A, 1)
    B = np.dot(np.linalg.inv(D-w*L),(1-w)*D+w*U)
    b = np.dot(np.linalg.inv(D-w*L), b)*w
    t_a, t_b = np.linalg.eig(B)
    if max(abs(t_a)) > 1:
        print('改矩阵发散')
        return

    k = 1
    while True:
        x1 = np.dot(B, x0) + b
        print('第%s次迭代' % k, x1)
        norm = max(abs(x1 - x0))
        # print(norm)
        if norm < mistake:
            break
        x0 = x1
        k += 1
if __name__ == '__main__':
    #系数矩阵
    print('实验一(1)\n')
    matrix_a = np.array([3.01,6.03,1.99,1.27,4.16,-1.23,0.987,-4.81,9.34],dtype=float).reshape(3,3)
    matrix_b = np.array([1,1,1], dtype=float).reshape(3, 1)
    print('原增广矩阵为:')
    print(np.hstack([matrix_a,matrix_b]))
    A=matrix_a
    L, U = LU_decomposition(A)
    print('矩阵L为：\n',L, '\n矩阵U为\n', U)
    x = answer(L,U,matrix_b)
    print('\n高斯求解：')
    gaussian(matrix_a,matrix_b)
    
    #这里为矩阵迭代
    #教材案列
    matrix_a = np.array([10,-2,-1,-2,10,-1,-1,-2,5],dtype=float).reshape(3,3)
    matrix_b = np.array([3,15,10], dtype=float).reshape(3, 1)
    jacobian(matrix_a,matrix_b)#雅克比迭代
    gauss_seidel(matrix_a,matrix_b)#高斯迭代
    sor(matrix_a,matrix_b,w=0.9)#松弛法迭代w=0.9

仁，上古

关注

15
点赞
踩
15

收藏

觉得还不错? 一键收藏
11
评论
列主元消去法和LU矩阵分解以及矩阵迭代求解

上三角求解wei#上三角求解：def up_solve(A,b): n = len(A[0]) A[n-1,n-1] = b[n-1][0]/A[n-1,n-1] for i in range(1,n): k = n-i-1#取逆序 A[k,k] = (b[k][0] - sum(A[k]*A[k+1]))/A[k,k] A[k]=np.hstack((A[k,:k+1],A[k+1,k+1:]))#拼接新数组 #.
复制链接

扫一扫