只要一杯秋天的奶茶，就能学会Python数值分析(2)

上节(https://www.jianshu.com/p/671a94ce586b)说到高斯消元法，今天从高斯列主元消元法开始，拓展到线性方程组的两种迭代方法：雅可比迭代法和高斯-赛德尔迭代法。同样，能力有限，希望读者指出我的问题，用代码和公式和我深入交流。

2.列主元消去法

参考的教材是《数值分析》(李庆扬等)的第148页到151页。

列主元消去法是使用于主元为0或很小的情况，思路是选取这一列中绝对值最大的值所在的行与当前行进行替换，具体可以参考教材。再写代码方面，有了第一节的东西，写列主元消去法就很简单了，只要在高斯消去法之前加上换行就行。代码如下：

def pivotTriangularMatrix(matA):

'''

先找主元，找到后跟高斯顺序消元操作一样

'''

numRows = matA.shape[0]

for row in np.arange(0,numRows-1):

Index = np.argmax(abs(matA[:,row]))

matA[[row,Index],:] = matA[[Index,row],:] # 交换。也可以用np.copy来做交换

# 交换行结束，开始消元

for k in range(row+1,numRows): # 在第row行的情况下，遍历剩下的行

if matA[k,row] != 0:

#第k(k小于row)行的新值等于该行减去第row行的值乘于一个系数。这个系数存在目的就是消元

matA[k,:]=matA[k,:]-(matA[k,row]/matA[row,row])*matA[row,:]

return matA # 此时输入的矩阵也被改变

用教材148页的例题4来测试一下这个代码。

test_mat3 = np.array([[0.001,2.0,3.0,1.0],[-1.0,3.712,4.623,2.0],[-2.0,1.072,5.643,3.0]],dtype=float)

result = Gauss_solve(pivotTriangularMatrix(test_mat3)) # Gauss_solve函数是上一节的函数

# 得出结果

result = array([[-0.49039646],

[-0.05103518],

[ 0.36752025]])

得出的结果和书中给出的答案是吻合的。暂时没有找其它算例测试。

3.今日加餐

今天舍友问我下面这个代码。刚好今天没去旁听。主要问这里面提到的独特的索引方式是什么。其实就是一组布尔值。这里举个例子，具体解释在代码里面。

i,n = 1,3

idx1 = np.arange(i,n+i)

'''

假设i=1,n=3

idx1 = array([1, 2, 3])

'''

idx1>n-1

'''

输出为 array([False, False, True])，这样就拿到了True所在的索引2

那么idx1[idx1>n-1]在这里就相当于idx1[2]，这个值等于3

'''

# 举例

test_mat1 = np.array([[2,0,6],[1,4,3],[3,2,1]],dtype=float)

test_mat[idx0,idx1]

'''

测试矩阵为：

test_mat1 = array([[2., 0., 6.],

[1., 4., 3.],

[3., 2., 1.]])

输出值为 array([0., 3., 3.]) 。实现了上图中取斜对角的操作

'''

二、线性方程组的迭代求解法

1.雅可比迭代法

该处参考的是《数值分析》(李庆扬等)第6章的6.1.1和6.2.1。书上关于雅可比迭代的内容有更多细节。这里我按照我对Gradient Descent的理解来编写该处代码。

这是一个迭代的计算方法，从我的理解来看，这类方法与线性方程组直接求解法有2个较大区别：1.对于解，会有一个初始值，会从这个初始值开始按照迭代公式运算下去。2.迭代法都要有迭代结束条件，要么是迭代解与真实解之间的差小于人为设定的一个阈值，要么是达到人为设定的最大迭代次数。

以下是我的代码，参考的是《数值分析》(李庆扬等)188页2.3这个公式，这个公式是一个展开的公式，我代码里写成矩阵运算的形式，这样可以减少使用for循环。

def Jocobi(A,b,initial,delta,maxTimes):

'''

输入：A是系数矩阵，N阶方阵

b是N*1列向量

initial是解的初始值，N*1大小

delta是人为设置的一个阈值，指到后期两次迭代之间解之间的差距已经很小了，用作迭代结束条件

maxTimes指人为设定的最大迭代次数

教材6.2.1将矩阵A分解成D,L,U，这里与之对应

输出：迭代后的解析解

'''

D = np.diag(np.diag(A)) # 获得D矩阵

L = np.tril(A,-1) # 获得L矩阵

U = np.triu(A,1) # 获得U矩阵

d = np.linalg.inv(D) # 对D矩阵求逆

G = -np.dot(d,L+U) # D的逆矩阵乘于(L+U)矩阵

f = np.dot(d,b) # D的逆矩阵b向量

X = np.dot(G,initial)+f # 初次的解

for i in range(maxTimes):

if np.linalg.norm(X - initial) > delta:

initial = X

X = np.dot(G,initial) + f

return X

利用第180页的例题1来对这个函数进行验证。解的初始值设为0解。

import numpy as np

A = np.array([[8,-3,2],[4,11,-1],[6,3,12]],dtype=float)

b = np.array([[20],[33],[36]],dtype=float)

initial = np.zeros((3,1))

'''

A = array([[ 8., -3., 2.],

[ 4., 11., -1.],

[ 6., 3., 12.]])

b = array([[20.],

[33.],

[36.]])

initial = array([[0.],

[0.],

[0.]])

'''

# 迭代求解

X = Jocobi(A,b,initial,0.001,10)

'''

X = array([[3.00028157],

[1.99991182],

[0.99974048]])

'''

进过10次迭代后解出的结果与书上181页的结果一致。

2.高斯-赛德尔迭代法

有了雅可比迭代的代码，写高斯赛德尔迭代代码就非常简单，你可能会惊讶，这两个代码怎么是一样的。看下图可以看到这两的区别。参考教材的6.2.2节。

以下是我的代码：

def Seidel(A,b,initial,delta,maxTimes):

'''

输入：A是系数矩阵，N阶方阵

b是N*1列向量

initial是解的初始值，N*1大小

delta是人为设置的一个阈值，指到后期两次迭代之间解之间的差距已经很小了，用作迭代结束条件

maxTimes指人为设定的最大迭代次数

教材6.2.1将矩阵A分解成D,L,U，这里与之对应

输出：迭代后的解析解

'''

# 先分解

D = np.diag(np.diag(A))

L = np.tril(A,-1)

U = np.triu(A,1)

d = np.linalg.inv( D + L ) # 这里是和雅克比不同的地方

G = -np.dot(d,U)

f = np.dot(d,b)

X = np.dot( G ,initial ) + f

for i in range(maxTimes):

if np.linalg.norm( X - initial ) > delta:

initial = X

X = np.dot( G,initial )+f

return X

同样用前面雅可比的算例来测试一下。

X = Seidel(A,b,initial,0.0001,10)

'''

X = array([[3.00000201],

[1.9999987 ],

[0.99999932]])

'''

同样的例题，在同样迭代10次的情况下，与前面雅可比迭代的结果有一点点偏差，但这都在正常范围之内。与书上的结果基本吻合。

今日小结

今天可太忙了，从早到晚。时间全靠课间挤，关于迭代法也没有考虑收敛条件的问题。这个系列对我来说还挺好玩，打算有始有终做下去。原来，写东西(还不用负责)这件事情能这么快乐。

今晚毕，明天继续。。。。。