第十四课线性联立方程的预处理共轭梯度(PCG)

系数矩阵病态

百度解释：求解方程组时如果对数据进行较小的扰动，则得出的结果具有很大波动，这样的矩阵称为病态矩阵
在直接求解的诸多方法中，线性方程可能无法准确地求解，尤其是在使用“手动”计算时。在这种情况下，方程组就被说成是“病态的”。现代电子计算器已经能够精确到小数点后许多位，因此方程组的处理并不像“手动”计算，只需要选择精确到小数点后几位。在下表中，四舍五入到小数点后4、5、6、7和8位的数值解与真实解进行比较。可以看出，对于一个适当的工程项目而言，精确到小数点后8位是必要的。当使用32位字的计算机时，建议同学们使用“双精度”，即64位字，达到大约15位小数的精度，以满足大多数实际用途。

本文描述的程序执行的计算精度通常为“双精度”，虽然可以去用数字去度量这个‘病态’，称为“条件数”，但这个过程是很辛苦的，而且并不适用在工程实践。到目前为止，最准确性的解决方案是对结果的物理一致性进行独立检查。
在迭代过程中，这种‘病态’会是的方程变得收敛速度慢或者是根本就不收敛。。因此，要看下面的方程

是否可以转化成一个等效的方程

之后迭代过程的收敛特性会变得更好。这种使用“预处理”矩阵[P]转化的过程称为“预处理”，得到

上面的叫做左处理

这种叫做右处理
如果[P]是[A]的逆，则不需要迭代，但计算工作却很难实现。而得到的另外一种有效方式是，寻找可以方便获得的[A]-1的近似值。其中最简单的是通过取[A]的主对角项的倒数形成的对角矩阵(存储作为一个向量)。这个想法可以通过在逆近似(称为“不完全因子分解”)中包含更多的[A]的主对角线项来扩展，但是在下面的程序中，我们只针对纯主对角线形式。对共轭梯度法的计算过程进行预处理，形成“预处理共轭梯度”或“PCG”算法。算法的步骤如下面所示。
初始阶段：

迭代阶段：

算例采用第一篇中的基础算例

具体求解过程可以看高斯赛德尔迭代
程序如下：
其中有一个主程序和一个检查是否收敛的子程序checkit
主程序：

#线性联立方程的预条件共轭梯度法
import numpy as np
import math
import B
n=3
converged=np.array([False])
d=np.zeros((n,1))
p=np.zeros((n,1))
u=np.zeros((n,1))
precon=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])
for i in range(1,n+1):precon[i-1,0]=1.0/a[i-1,i-1]
b[:]=b[:]-np.dot(a,x)
d[:]=precon[:]*b[:]
p[:]=d[:]
print('前几次迭代值')
iters=0
while(True):iters=iters+1u[:]=np.dot(a,p)up=np.dot(np.transpose(b),d)alpha=up/np.dot(np.transpose(p),u)xnew[:]=x[:]+p[:]*alphab[:]=b[:]-u[:]*alphad[:]=precon[:]*b[:]beta=np.dot(np.transpose(b),d)/upp[:]=d[:]+p[:]*betaif 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

终端输出结果：

程序结果与计算结果一致

预处理共轭梯度(PCG)解线性联立方程(python,数值积分)相关推荐

线性联立方程的雅可比迭代解(jacobi Iteration).(python,数值积分)
第八课线性联立方程的雅可比迭代迭代方法在前面几节中,已经描述了线性代数方程组的"直接"解法.我们所说的"直接"是指求解方法是通过固定次数的运算得到答案.由 ...
python实现共轭梯度算法
python实现共轭梯度优化算法一.共轭梯度算法简介二.实现共轭梯度方法的两块重要积木 1.共轭方向的确定 2.方向优化步长的确定 note 三.共轭梯度算法优化过程四.python实现共轭梯度 ...
线性联立方程的高斯赛德尔迭代(Gauss-Seidel iteration)(python，数值积分)
第九课线性联立方程的高斯赛德尔迭代在雅可比迭代中,{x}k+1(在上一课中被称为xnew)的所有分量都是使用{x}k(在上一课中为x)的分量推出的.因此,新的向量值完全是根据旧的向量值获得的. 然 ...
【优化理论】共轭梯度下降算法实现
实验目的实验步骤此次实验要求采用共轭梯度下降来解决二次优化问题,在解决此问题前,首先分析下什么是共轭梯度法,共轭梯度法能解决什么问题. 我们来看一个线性方程组Ax=b,求解此方程组的过程可 ...
共轭梯度算法之FR算法
共轭梯度算法之FR算法引理 FR算法算法步骤引理若f(x)=12xTGx+δTx+γf(x)=\frac{1}{2}x^TGx+\delta^T x+\gammaf(x)=21xTGx+δT ...
一文详解线性最小二乘与非线性最小二乘
一文详解线性最小二乘与非线性最小二乘一.最小二乘法的引出二.线性最小二乘法 1.线性最小二乘的描述 2.线性最小二乘特殊情况的求解 3.线性最小二乘一般情况的求解三.非线性最小二乘法 1.非线性 ...
matlab ncg,2步非线性共轭梯度(NCG)法
2步非线性共轭梯度(NCG)法提出一种求解非线性方程组的2步非线性共轭梯度法(2步NCG法),并在一定条件下证明 (本文共9页) 阅读全文>> 作者提出的自适应共轭梯度方法是对共轭梯度方 ...
python英语字典程序修改_详解如何修改python中字典的键和值
我们知道python中字典是无序的,它们都是通过hash去对应的.一般的如果我们需要修改字典的值,只需要直接覆盖即可,而修改字典的键,则需要使用字典自带的pop函数,示例如下: t = {} t['a ...
python java混合编程_详解java调用python的几种用法(看这篇就够了)
java调用python的几种用法如下: 在java类中直接执行python语句在java类中直接调用本地python脚本使用Runtime.getRuntime()执行python脚本文件(推荐 ...

预处理共轭梯度(PCG)解线性联立方程(python,数值积分)

第十四课线性联立方程的预处理共轭梯度(PCG)

系数矩阵病态

预处理共轭梯度(PCG)解线性联立方程(python,数值积分)相关推荐

最新文章

热门文章

预处理共轭梯度(PCG)解线性联立方程(python,数值积分)

第十四课 线性联立方程的预处理共轭梯度(PCG)

系数矩阵病态

预处理共轭梯度(PCG)解线性联立方程(python,数值积分)相关推荐

最新文章

热门文章

第十四课线性联立方程的预处理共轭梯度(PCG)