只要一杯秋天的奶茶，就能学会Python数值分析(2)
上节(/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次的情况下，与前面雅可比迭代的结果有一点点偏差，但这都在正常范围之内。与书上的结果基本吻合。
今日小结
今天可太忙了，从早到晚。时间全靠课间挤，关于迭代法也没有考虑收敛条件的问题。这个系列对我来说还挺好玩，打算有始有终做下去。原来，写东西(还不用负责)这件事情能这么快乐。
今晚毕，明天继续。。。。。
