代数多重网格算法

代数多重网格算法∗
2007年2月14日
1基本思想
Gauss-Seidel算法的特点是,最初几步收敛的很快,但是很快就开始停滞不前.到最后几乎不收敛.从数值试验的图像可以看出,Gauss-Seidel迭代当插值点少的时候,收敛速度极快,但当插值点多的时候,由于上述效应收敛速度极慢.因此,代数多重网格(Algebraic Multi-Grid)算法利用这些特点,将由具体方程离散出来的矩阵,重投到一系列由细到粗的网格上,在每一层网格上只做若干次Gauss-Seidel迭代.与传统的多重网格算法不同,该算法不需要提供任何网格的信息.所有的信息完全只来自方程离散后的矩阵.
假设Possion方程
−∆u=f(1)用某种离散方法(比如,有限元或有限差分),在某个相当细的网格上,最后产生线性问题
A x=b.(2)现在考虑如何将其投影到一个较粗的网格上.
假设φ={φi},i=1,···,N为细网格上的一组分片一次线性有限元基函数.则矩阵A是一个N×N的矩阵,且元素a ij可以看作是对应基函数的一个双线性运算a(φi,φj).我们如果要将A重新投影到一个对应基函数为ψ={ψi},i=1,···,M,M<<N的粗网格上,则根据用φ表出ψ的关系,我们可以得到
ψ=Pφ,(3)这里P是一个M×N的矩阵.相应的,如果令
˜A=P AP T,˜x=P x,˜b=P b,(4)则
˜A˜x=˜b(5)∗该文档为李若教授讲授的《数值分析高等算法》的课堂笔记,由王何宇整理.
就可以看作是(2)投影到粗网格上以后的问题.代数多重网格的做法,就是对第k步的线性问题
A k x=b k(6)先用Gauss-Seidel迭代进行几步迭代,得到一个近似解x k,然后将残问题
A k x=b k−A k x k(7)用投影矩阵P k重投到粗一层的网格上得到第k+1步的问题
,b k+1=P k b k−A k x k,(8)
A k+1x=b k+1,A k+1=P k A k P T
k
如此不断迭代和重投,直到得到一个规模相当小的线性问题后,可以用直接法(Gauss消去法)求得精确解,然后用记录下的一系列P k矩阵,还原出原问题的解.在还原的时候,仍然使用Gauss-Seidel迭代在每一层来改进数值解.如此整个过程为一步AMG迭代.
2算法步骤
现在给出严格的算法步骤.
对过程AMG(A k,x k,b k),
第一步如果A k的阶数小于一个给定的整数,比如20,则用Gauss消去法解出并x k并返回;否则,对问题(6)做3至5步Gauss-Seidel迭代;
第二步产生问题(8),令x k+1=0;
第三步递归调用AMG(A k+1,x k+1,b k+1);
第四步x k=x k+P T
x k+1;
k
第五步做3至5步Gauss-Seidel迭代,返回.
所以对问题(2),执行AMG(A,x,b)完成一次AMG迭代(迭代更新了x).而整个求解过程为
do AMG(A,x,b)while(|b-Ax|<e).
这里e为控制误差.
现在整个算法过程,还剩下的问题就是如何产生P.假设(2)是在网格点x={x i},i= 1,···,N上离散的,那么我们首先要确定一个粗网格,也就是说,要考虑在x移除一些点,保留一些点.我们将保留的点称为核心点(core points).核心点的选取,应该满足如下两点:
1.不能太多:核心点彼此之间,不能相邻;
2.不能太少:所有被移除的点,必须至少与一个核心点相邻.
根据这个原则,和稀疏矩阵存储规则,我们设计算法如下:
第一步产生一维数组c[1:N],并置所有元素初值为零;
第二步选出核心点:for i=1:N,如果c[i]==0,则x i为核心点,同时将所有满足a ij=0的c[j]+=1;
第三步假设x i为第k个选出的核心点,则P的第k行元素为:0如果a ij==0,
1
c[j]
如果a ij=0.
注意该算法仍然只用了矩阵的信息而没有使用网格的信息.
3算法分析
该算法实际上的迭代过程完全基于Gauss-Seidel迭代.所以其收敛的要求和Gauss-Seidel一样,为矩阵对称正定.但是要达到加速收敛的效果,要求矩阵必须有网格和方程的背景,否则这样做是没有意义的.。