09

2024-07-15 10:14| 来源: 网络整理| 查看: 265

旁听了今天的上机课，收获良多。

方阵A求逆，先做LU分解。A的逆等于U的逆乘于L的逆，L的逆就利用下三角矩阵求逆算法进行求解，U的逆可以这样求：先将U转置成下三角矩阵，再像对L求逆一样对U的转置求逆，再将得到的结果转置过来，得到的就是U的逆。因此，关键是下三角矩阵的求逆。

1.下三角矩阵求逆算法

我利用的公式计算公式如下：

对角元素.png 对角元素以下的元素.png

我的代码如下：

def triInverse(matA): ''' @author：zengwei 输入： matA:一个等待求逆的下三角矩阵,大小为n*n，并且希望里面的元素值是浮点数输出： matInv:matA的逆矩阵 ''' numRows = matA.shape[1] matL = matA.copy() matInv = np.zeros((numRows,numRows)) for row in np.arange(0,numRows): matInv[row,row] = 1/matL[row,row] for k in np.arange(row-1,-1,-1): matInv[row,k] = -(np.dot(matInv[row,k+1:row+1],matL[k+1:row+1,k]))/matL[k,k] return matInv

同样，生成一个矩阵来测试一下上面的函数，并得到输出。

import numpy as np test_mat = np.array([[1,0,0,0],[2,3,0,0],[4,5,6,0],[7,8,9,10]],dtype=float) ''' test_mat = array([[ 1., 0., 0., 0.], [ 2., 3., 0., 0.], [ 4., 5., 6., 0.], [ 7., 8., 9., 10.]]) ''' Inv = triInverse(test_mat) ''' Inv = array([[ 1. , 0. , 0. , 0. ], [-0.66666667, 0.33333333, 0. , 0. ], [-0.11111111, -0.27777778, 0.16666667, 0. ], [-0.06666667, -0.01666667, -0.15 , 0.1 ]]) '''

显然，我不知道解出来对不对，但是我可以用这个逆矩阵Inv和测试矩阵test_mat相乘看看是不是单位矩阵来判断。

np.dot(test_mat,Inv) ''' array([[1., 0., 0., 0.], [0., 1., 0., 0.], [0., 0., 1., 0.], [0., 0., 0., 1.]]) '''

是单位矩阵，说明下三角矩阵求逆成功。接下来，利用上面的函数来进行矩阵的求逆。

2.矩阵求逆

首先，先贴出我LU分解的函数：

def getLandU(A): ''' @author：zengwei 输入： A:系数矩阵；array格式，且数值类型必须为浮点数，否则会出现截断误差。输出： LU分解中的L和U矩阵 ''' matA = A.copy() if matA[0,0]==0: print('不符合要求！需要换行操作。') else: numRows = matA.shape[0] matL = np.identity(numRows) # 准备给L矩阵用 for row in np.arange(0,numRows-1): # 前n-1行，row代表第几行 for k in range(row+1,numRows): # 在第row行的情况下，遍历剩下的行 matL[k,row] = matA[k,row]/matA[row,row] # 获得L矩阵 if matA[k,row] != 0: matA[k,:]=matA[k,:] - (matA[k,row]/matA[row,row])*matA[row,:] return matL,matA # matL为L矩阵，matA变为U矩阵

然后我们利用getLandU函数和triInverse函数来写矩阵求解函数。如下：

def matInverse(A): ''' @author：zengwei 输入： A:想要求逆的系数矩阵，n*n，希望里面的数值是浮点数输出： A的逆矩阵 ''' L,U = getLandU(A) #LU分解 L_inv = triInverse(L) #L求逆 U_inv = triInverse(U.T).T #U求逆 return np.dot(U_inv,L_inv)

下面，我们生成一个随机矩阵来测试，并验证求得的逆矩阵和原矩阵相乘是否为单位矩阵。

TestMat = np.float64(np.random.randint(0,10,(4,4))) ''' TestMat = array([[3., 7., 4., 0.], [8., 3., 9., 3.], [5., 4., 2., 4.], [7., 0., 7., 6.]]) ''' TestMatInv = matInverse(TestMat) isI = np.int64(np.dot(TestMatInv,TestMat)) # 验证 ''' isI = array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]], dtype=int64) '''

可见，逆矩阵乘与原矩阵是单位矩阵，这里用了一下np.int64()函数，是为了让结果好看。否则原本是0的地方就会是非常非常非常小的数，这是由于计算机用浮点数运算得到的效果。

放假。

【本文地址】

09

09

今日新闻

推荐新闻