线性方程组的定常迭代法简析——by Tangwei
在采用有限元或有限体积进行力学的数值计算时,会涉及到大型的线性方程组,因为我们将比如流体离散为一个个小的有限体积,通过这种离散的方式把大自然中原本连续的物质给分割成一小块一小块,进而可以利用计算机进行计算,避免了对现实世界中物体的解析解有时求解太过复杂,甚至解不出来。
那么在计算离散值时,会涉及到大量的数据,因为模型会被分割成成千上万、甚至百万千万个元素。它们将组成一套线性方程组通过计算机来进行计算。我们将现实中的数据抽象成一个线性方程组,如下:
A=()∈R^(n×n)是非奇异矩阵,x∈R^n是n维向量,b∈R^n是n维向量。若b不全为0向量时,这是一个典型的非齐次方组。
那么当未知数的系数矩阵A是低阶稠密矩阵的时候(n比较小),我们可以用LU分解、高斯消元法等,一个一个去算把它解出来。但是当A是大型的稀疏矩阵的时候,就是说n很大,而0比较多的时候,一个个地去求解就不现实了,这就要用迭代法去求解。
我们把含矩阵A的(1)式变化成
的形式,这种方式是怎么得到的呢,B和f 又是啥东西?来看:
将A分裂成
M是供你选择的一个非奇异矩阵,而且对于Mx=d这种形式的方程组是容易求解的(这个d只是一种形式,就是说式子容易求解就行),我们就把M叫做分裂矩阵。那么这个时候
所以上面我们讲到的B和f 是啥东东,就是这个:
所以我们叫B为迭代法的迭代矩阵。那么不同的M选取,就产生了不同的迭代矩阵,导致迭代方法也不一样。以后会讲到具体的迭代方法:雅克比迭代、高斯-赛德尔迭代和逐次超松弛迭代SOR。本次仅对基础概念进行初步分析。
好了,话又回来了,变成x=Bx+f干嘛?
干嘛?擎好了您嘞!我们在用计算机算的时候,计算机用的都是一个个离散的数据,它也不会像我们人那样有思想地去思考怎么解题,只能根据代码的指令去机械的计算。那怎么办,用迭代法,这种体力活计算机最适合。
我们先给定方程组x=Bx+f的解,设为唯一的解x*,那么
x*=Bx*+f (6)
这个“=”是两边相等时的等号,因为x*是真解.
我们先给定x一个初始向量,通常为x^(0)=(0,0…,0)^T,x的上标(0)是迭代的初始值,带入到方程组x=Bx+f就解得一个x^(1)解向量。这样就得到:
k表示迭代的次数
注意到,这里不管迭代多少次,B和f 始终不变的。
那么问题来了,这样迭代下去,x^(k)会向真值x*靠拢,而收敛于它吗?这个问题的答案是:不一定!
我们把
这个极限若存在,就称此迭代法收敛,收敛的值就是x*,否则就是发散。
那么进一步地,就要研究{ x^(k)}的收敛性了,我们这是要引出一个量:误差向量。
它的意思就是第k+1步时,x向量与真值向量的差值。
由(9)式减去(6)式,得到
我们要确定{ x^(k)}是否收敛,B满足什么样的条件,可以使得
也就是需要:
上面我们讲到,B满足(13)式时迭代收敛,那么什么样的条件下能得出这样的B矩阵,以及有这样的B矩阵会得到什么样的结果?(3)式讲到把矩阵A分裂为A=M-N的形式,就得到(2)式来求解线性方程组。这样一来我们就可以构造一个一阶定常的迭代法了:
下面我们有这样几个结论(定理):
结论(1):对于任意选取的初始向量x^(0),迭代方程组(14)收敛的充要条件为矩阵B的谱半径:
在给出其证明之前,我们先给出另外几个结论,
结论(2):
结论(3):
对于矩阵序列,
结论(4):
下面我们稍证明一下结论(1):
1)先证明收敛的充分性。
已知谱半径,由(2)式得到(I-B)x=f,即Ax=f,
∴ I-B为非奇异矩阵,因此(I-B)x=f,即Ax=f有唯一解,设为x*,得到:
那么误差向量
由(B)<1可知,根据上述结论(2)得出
那么对于任意的x^(0),
2)证明必要性
已知迭代法(14)收敛,所以可以设对任意x^(0)都有:
所以x*是方程组的解,并且对任意x^(0)都有:
那么对于任意的x^(0)其实就是任意的^(0),都有上式成立,那么根据结论(3)可知:
再由结论(2)知道:
以上结论(1)给出了线性方程组迭代法的收敛性的充要条件以及证明,那么不同的B,虽然都收敛,但是收敛的速度却是不一样的,对于计算机计算来说,收敛的快,就可以在较短的时间步得到结果,也节省计算资源。
我们来看迭代法的收敛速度:
假设(14)是收敛的,那么
得到:
根据算子范数的定义:
我们来看一下平均每次迭代后的误差向量的压缩率:
令第一次迭代,前后步向量误差的比率:
那么第二次为:
直至第k次为:
令P为平均压缩率,将,
,…,
相乘:
这里的每个P1,...Pk,不是把每个项的分子分母相消得到总的压缩率的,这里因为取max,所以不能相消,此处只是为了展示平均压缩率的数值得到的路径。平均值P的k次方就是等于每个压缩率的总积,进而等于总压缩率的,从而得到了下面(15)式。
所以平均每次迭代的压缩率P为:
我们希望迭代k次后,误差向量^(k)小于初始误差向量的
倍,即
根据算子范数和上式,就可以得到:
其中设σ<<1,我们可以用10的幂来表示,比较容易理解,取=10^(-s).
又因为矩阵B的谱半径(B)<1,所以根据上述结论(2)得:
由此可以看出,迭代次数k与如下参数:
有关。分母就可以看成是迭代速度,在s一定的情况下,k取决于分母的迭代速度。
我们再来看看两个收敛速度:平均收敛速度和渐进收敛速度。
平均收敛速度的定义为:
因为这个收敛速度和k及范数有关,为了方便计算和分析,我们把(B)进行变换,有结论(4)可知
而对数是连续函数,所以
它的意思就是在迭代次数趋向于∞时的迭代速度,这样就与k无关了。
渐进收敛速度的定义为:
这个收敛速度就与k及范数无关了,它反应的是迭代次数趋向于∞时迭代法的渐进性质的速度。由(18)式可知,当矩阵B的谱半径(B)越小时R(B)越大,收敛的就越快。此时我们可以用(19)式来估计所需的迭代次数:
所以,在s选定后,k的次数,也就是收敛的快慢由渐进收敛速度确定,即由B的谱半径决定,那么选取什么样的B对迭代的收敛速度来说,就是一个很重要的问题了。