Python数值计算（5）

这部分讨论一下矩阵的求逆和线性方程组的求解。以下的操作需要引入numpy模块下的子模块linalg，也就是linear algebra：

from numpy import linalg

矩阵求逆

如果两个矩阵A，B的矩阵乘积为单位矩阵，则互为逆矩阵，且有AB*BA=I。

计算矩阵的逆可以使用inv函数：

A=np.array([[2,1,-1],[1,0,-2],[1,1,3]])
B=linalg.inv(A)
print(B)
'''
[[-1.   2.   1. ]
[ 2.5 -3.5 -1.5]
[-0.5  0.5  0.5]]
'''
print(A@B)
'''
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
'''

并不是所有的矩阵都有逆矩阵，行列式为0的矩阵就没有逆矩阵，计算行列式的函数为det：

A=np.array([[2,1,-1],[1,0,-2],[1,1,3]])
print(linalg.det(A)) # -2
C=np.array([[2,1,-1],[1,0,-2],[1,1,1]])
print(linalg.det(C)) # 0
Cinv=linalg.inv(C) # raise LinAlgError("Singular matrix")

通常， 一个标准的线性方程组具有如下形式:

其中A为矩阵，x和b为列向量。

例如求解线性方程：

，

一个很经典的求解方法是克莱姆法则，实现代码如下：

def cramerSolver(A:np.ndarray,b:np.ndarray):
    '''
    A为n×n方阵，返回值x，b为n×1矩阵
    '''
    n=b.shape[0]
    x=b.copy()
    for i in range(n):
        B=A.copy()
        B[:,i]=b[:,0]
        x[i,0]=np.linalg.det(B)/np.linalg.det(A)
    return x

求解调用该函数如下：

A = np.array([[1,2,1],[3,-1,0],[0,3,-2]])
b = np.array([[7],[4],[4]])
X0=cramerSolver(A,b)
print(X0)
'''
[[2]
[2]
[1]]   
'''

事实上这种方式由于计算量太大，通常不大会使用，但是这个法则确实很优雅。

另一种方式是使用逆矩阵。如果在等式两边同时左乘以A的逆矩阵，则：

还有一种方式是使用linalg下的solve函数，上述两种方式的求解可以编码如下：

A = np.array([[1,2,1],[3,-1,0],[0,3,-2]])
b = np.array([[7],[4],[4]])
X1=np.dot(np.linalg.inv(A),b)
X2 = linalg.solve(A,b)
print('X1=',X1)
print('X2=',X2)

其结果是：

X1= [[2.]
 [2.]
 [1.]]
X2= [[2.]
 [2.]
 [1.]]

所以方程组的解为