第十四课 线性联立方程的预处理共轭梯度(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