第十一课 最陡下降法解方程
梯度法
前面描述的几种方法有时被称为“平稳”方法,因为它们没有根据试解中误差的值来修改收敛过程。
相比之下,在梯度法中,误差被重复评估,并生成新的试解。在我们的第一个例子中,我们将会针对正定系数矩阵的方程。在这种情况下,梯度方法编程非常简单,但它们背后的数学推理有些复杂。感兴趣的同学可以看一些详细介绍梯度法的书籍。
最陡下降法
对于任何试用解 x ,误差或“残差”可以表示为
我们从解向量{x}0开始计算。在最陡下降法中,根据以下的算法中使{r}k中隐含的误差最小化:
按照上面的方法迭代求出误差值r,从而求出解向量x。
本篇的例子和前两篇的例子一致,具体可以看高斯赛德尔迭代
程序如下:
其中有一个主程序,一个检查是否收敛的子程序checkit
主程序:
#线性联立方程的最陡下降法
import numpy as np
import math
import B
n=3
converged=np.array([False])
r=np.zeros((n,1))
u=np.zeros((n,1))
xnew=np.zeros((n,1))
a=np.array([[16,4,8],[4,5,-4],[8,-4,22]],dtype=np.float)
b=np.array([[4],[2],[5]],dtype=np.float)
x=np.array([[1],[1],[1]],dtype=np.float)
tol=1.0e-5
limit=100
print('系数矩阵')
print(a[:])
print('右手边向量',b[:,0])
print('初始猜测值',x[:,0])
print('前几次迭代值')
r[:]=b[:]-np.dot(a,x)
iters=0
while(True):iters=iters+1u[:]=np.dot(a,r)alpha=np.dot(np.transpose(r),r)/np.dot(np.transpose(r),u)xnew[:]=x[:]+alpha*r[:]r[:]=r-alpha*u[:]if iters<5:print(x[:,0])B.checkit(xnew,x,tol,converged)if converged==True or iters==limit:breakx[:,0]=xnew[:,0]
print('到收敛需要迭代次数',iters)
print('解向量',x[:,0])
checkit
def checkit(loads,oldlds,tol,converged):
#检查前后两个量的收敛性neq=loads.shape[0]big=0.0 converged[:]=Truefor i in range(1,neq+1):if abs(loads[i-1,0])>big:big=abs(loads[i-1,0])for i in range(1,neq+1):if abs(loads[i-1,0]-oldlds[i-1,0])/big>tol:converged[:]=False
终端输出结果:
程序结果与计算结果一致。