多极展开广义极小残值算法

1、第4章 基于多极展开法的广义极小残值算法线性方程组的求解的诸多方法中，传统的Gauss消去法或是由Gauss消去法派生的其他方法解方程所需的计算次数与方程秩N的3次方成比例，而在Arnoldi算法67的基础上提出的GMRES迭代算法用于求解稠密系数矩阵线性方程组时，其计算次数与N成比例。GMRES方法因其存储量、计算量较少，每步迭代具有最优性并且迭代不会中断等优点，成为较为可靠的Krylov子空间投影法。但是当问题规模较大时，使用GMRES法求解依然很困难。故将其与多极展开法(FMM)相结合更新计算结构，提高求解效率。其要点是求解时可只用行矩阵求解，而无须生成系数总矩阵，故可减少运算量及计算所

2、需内存，提高计算速度。4.1 广义极小残值(GMRES)算法 Krylov子空间方法设待求线性代数方程组为 (4-1)式中，为阶非对称非奇异稠密实矩阵，为实数组成的向量。设是维线性空间中的一组线性无关向量，由其张成的线性子空间记为，。Krylov子空间投影方法的基本思想是在子空间中构造满足(4-2)的，即式(4-1)近似解。 (4-2)式中，是另一子空间中的任一向量。令，其中为一维向量，则式(4-2)可写为 (4-3)故若令，则近似解可表示为 (4-4)其中，若选取，这样选取的Krylov子空间方法称之为Arnoldi算法。近似解的构造方法如下：首先，取为任一向量，令，则可将(4-1) 式化为

3、 (4-5)之后，特别取，假设其中各向量线性无关，再由其构造一组标准正交基，正交化过程的计算流程为：(1)任取，计算和(2)迭代 记迭代中形成矩阵为，且有，上述正交化过程也叫Arnoldi过程。下一步，构造近似解，其推导过程如下：令，因为又因为由子空间的一组标准正交基组成，故有，由第2步迭代可推出3个重要关系式 (4-6) (4-7) (4-8)(4-6)式和(4-7)式在算法中已直接给出，(4-8)式可由(4-7)式导出，在 (4-7) 式两边乘得 (4-9) (4-10)即 (4-11)两边取模 (4-12)即 (4-13)(4-6)式可写成在(4-7)式中令则，所以 (4-14)推导完毕

4、。最后，求解最小二乘问题，判断并使求得的近似解满足精度要求。 GMRES算法但是上述Arnoldi算法在理论上难以分析其收敛性，并且存在中断问题，所以人们转而探索其它的Krylov子空间投影方法，如果选取，这样的方法称之为GMRES算法。此时，在中极小化就等价于在中极小化。原问题归结为最小二乘问题。GMRES算法的计算步骤可以归结为：(1)初始化 任取，计算；(2)迭代 For （直到满足条件）doStep1 正交化Step2 标准化 Step3 更新与 为上Hessenberg矩阵，当时第一列省略，并且有；(3)求解最小二乘问题，得到；(4)构造近似解。 GMRES算法的实用化处理GMRES

5、算法应用于求解边界元方程组目前来说是一种比较理想的算法，但是经典的算法存在着两大不足:第一，就是随着迭代次数的增加，由Arnoldi过程生成的正交基也不断增加，从而造成存储量和计算量的不断增加；第二，就是随着迭代次数的增加，由于误差积累造成正交基正交性的丧失，从而使迭代的收敛速度放缓甚至发散。因此GMRES实用化的关键就是解决这两个不足。对于第一个不足，传统方法采用重启技术，其基本思想就是计算过程中最多存储m个基向量，每隔m次迭代就重启一次，即用最后一次的余量作为新一轮迭代的初始余量，最后一次的迭代解作为新一轮迭代的初始解，即GMRES(m)算法。但是这种方法不仅要根据经验确定m的值，而且也没

6、有收敛性保证，经验表明该法经常出现停滞甚至发散的现象，为此一些学者做了大量的工作68。但实践表明这些方法应用于三维弹性静力问题的边界元法，效果依然不够理想，反而非重启型GMRES取得了比较好的效果69。而为了能使非重启型GMRES实用化，必须尽量减少迭代的次数，以减少存储量和计算量。本文采用的方法就是在经典算法的基础上引入预条件技术。对于非对称矩阵，预条件的目的就是使预条件后的矩阵的特征值的分布尽量邻近。目前最常用的预条件技术有两类，一类是左预条件，另一类是右预条件，它们的变换式为：左预条件 右预条件 上两式中矩阵称为预条件阵，它的选取应符合两条原则，第一要使预条件处理后的矩阵的特征值分布尽量

7、靠近，第二就是便于求解。在边界元方程组的求解中，经常取为系数矩阵的对角元素组成的对角矩阵，称为对角预条件阵，或者取由对角线元素和次对角线元素组成的带状矩阵，称为带状预条件阵。变换后的方程组记为。对于左预条件，；对于右预条件，。具体算法的实施就是在计算过程中用预条件处理后的方程组代替原方程组，再进行一定的整理。对于第二个不足，可采用双精度存储和运算或修正的Gramm-Achmidt正交化过程。为提高所得向量组的正交性，本文采用前者使正交程度大大提高，同时加快收敛速度。为避免存储量的增加，可以只引入少量双精度计算，对及相应计算采用双精度存储，精度混合运算时先将单精度转化为双精度，再用双精度计算。4

8、.2 快速多极展开法(FMM)4.2.1 FMM基本思想 FMM源于静电场计算，用于计算大量粒子间两两相互作用的电场解析。基本点在于将粒子按空间位置分为不同的集合，分属于距离较近的集合及集合内的粒子间的相互作用两两直接计算；当集合相距足够远时，集合之间多个粒子的相互作用根据多极展开公式近似计算。空间中多个粒子对某一粒子作用叠加的直接计算式为 (4-15)式中，是一组影响系数。(4-15)式在数学、物理、化学及生物等不同学科中表达形式相近，普遍用于化学中分子动力学和量子力学模拟、天文学中大规模引力系统的计算、电子工程中的电容和电感的计算和不可压缩液体的动力学分析70-73。可见，在使用上述公式求

9、解电荷间的相互作用问题时，当需要同时计算的电荷数增加时，计算量将以阶上升。由Greengard等提出的快速多极展开法，能够基于球谐函数在空间中的多极展开，采用递归算法结构，将多个粒子对某处的作用层层近似转化为单一粒子的作用，将计算量降为阶。直接计算如图4-1所示，当距离足够远时，使用FMM方法相当于把Q区域内的所有粒子对位于P点的粒子的作用近似归结为域内某一点对P点粒子的作用，如图4-2所示。图4-1 直接计算Fig. 4-1 Direct Calculation图4-2 采用多极展开法计算Fig. 4-2 Calculation with FMM在多层的递归算法中，对给定的区域，首先定义包含

10、所有源点的最小立方体(或平面)作为计算区域。先将计算区域按某种方式层层细化为越来越小的多个子区域。在0层，有一个立方体包含所有的计算主域。第层是把第层的每个立方体分成若干份形成的。一个立方体是否需要再细分，是根据其中包含源点的个数是否少于某一设定值，这个值是根据解题规模的大小而确定。之后再从最后一层开始，将相互作用力的计算分成两个部分，第一部分是与该立方体相邻的立方体，即有公共的边界点，它们之间的作用力需要直接计算。第二部分是与该立方体不相邻的立方体，即没有公共的边界点，它们之间的作用力用FMM计算，将其上的每一单元格中所有粒子的作用归结为某一粒子的作用，之后再进行其上一层的FMM近似计算，如

11、此递归下去直到最上层。其递归计算结构如图4-3所示。图4-3 多层递归FMMFig. 4-3 Multi-level FMM4.2.2 FMM相关定理 FMM的所有公式，建立在球坐标系上。为适用FMM，定义阶次的球面调和函数：，式中，是连带勒让德函数，可以用罗巨格(Rodrigues)公式定义式中，是次勒让德多项式。下面给出的一系列FMM相关定理适于域积分项的计算。其中，多极展开和局部展开定理，描述了电荷集在远处的电场只与电荷集等效的总电荷(球心处)大小有关，而与其分布形式无关，这与力学的圣维南原理相当。FMM主要依赖多极展开和局部展开的平移。相关定理叙述如下：(1)多极展开定理假设有个强度为

12、的电荷分别位于球坐标分别为的点且这些点均位于球心在原点，半径为的球内, 如图4-4所示。则在任意点，当时，由电荷产生的位势为 (4-16)式中。此外，对任意的有图4-4 多极展开Fig. 4-4 Multipole Expansion(2)局部展开定理 假设有个强度为的电荷分别位于点上，其坐标分别为，如图4-5所示。且点位于球心在原点，半径为的球外。则在任意点，当时，由电荷产生的位势为：，式中，。此外，对任意的有 (4-17)图4-5 局部展开Fig. 4-5 Local Expansion(3)多极展开的平移假设有个强度为的电荷位于球心在，半径为的球内。则由这些电荷在任意点产生的位势为 (4

13、-18)式中，是向量的球坐标。由此，对任意位于球心在原点，半径为的球的外部的点有式中式中 (4-19)对任意的有，如图4-6所示。图4-6 多极展开平移Fig. 4-6 Transiton of the Multipole Expansion(4)局部展开的平移假设是中的一对点，是向量的球坐标，为自然数，为阶局部展开的中心。如图4-7所示。在点处的位势表示为 那么，在中的任意点有 式中。图4-7 局部展开的平移Fig. 4-7 Transiton of the Local Expansion(5)多极展开到局部展开的变换假设个强度为的电荷位于球心在，半径为的球内，并且有。则有对应的多极展开式(

14、4-18)在球心在原点，半径为的球内收敛。因此，对在任意的点，由电荷产生的位势用局部展开式写成 (4-20)式中，由(4-19)式定义，此外，对任意的有 (4-21)如图4-8所示。图4-8 多极展开到局部展开的变换Fig. 4-8 Transformation from the Multipole Expansion to the Local Expansion4.2.3 FMM的使用条件用FMM级数展开式计算粒子集合之间相互作用时，必须保证参与运算的集合相距足够远，否则将严重影响其计算精度。如图4-9所示，假设有两组点电荷分别位于点集和，在下面条件下，定义点集的相互独立性。.x1r.x3.

15、xm.x2.x0.y1.y2.y3.yny0rr图 4-9 点集示意图Fig. 4-9 Diagrammatic Sketch of the Set of Points点集独立定义 若存在和实数，且 则称，点集与是相互独立的，或者称点集与相距足够远。用表示中心在，半径为的球，表示中心在，半径为的球，根据三角不等式，若，则对任意有点处由处电荷产生的为实可表示为，取截断阶数，截断误差。是表征源粒子和场粒子之间间距的量。设两集合中心的距离为，则截断误差为由上式得4.3 基于FMM的GMRES算法随着对边界元等问题研究的进一步深入，由于研究对象结构的复杂和尺寸的巨大，为获得准确的分析结果，需要用较多的

16、单元数去离散该结构，所获得的影响系数矩阵将为一个阶数巨大系数稠密的矩阵，若用Gauss消去法一类的直接解法求解，所需的内存容量需以几十GB计算，运行时间需以天计算，这是不现实的。同Gauss消去法相比，GMRES 方法所需计算时间由O() 减少至O()，但所需内存容量相同，皆为O()。建立在FMM基础上的GMRES 方法将计算所需内存减少到了O()，并提高了计算速度，为这一难题的解决提供了可能。经分析可以发现GMRES算法的运算量主要来自两个方面：(1)在开始迭代前，首先要初始化形成个矩阵元素(必需)；(2)在形成正交基的每次迭代过程中，计算矩阵与向量积需要次运算。可见，提高GMRES 算法计

17、算效率的关键是矩阵与向量积运算的加速。将FMM引入GMRES的基本思想，即隐式构造一个与原离散积分方程相联系的稠密矩阵的稀疏表示，再将这一稀疏表示用于GMRES方法所需的矩阵向量积的快速计算中，以减少运算量及计算所需内存，提高计算速度。在这里，我们将矩阵与向量的乘积理解为计算给定强度的奇点分布在空间任意一点所产生的势。根据前面对FMM使用条件介绍，给定奇点对空间任意一点的影响应首先分成近场影响和远场影响两部分，之后远场部分影响可采用有效的近似计算势的的方法近似计算。根据多极展开定理，在足够远场的近似计算中，该影响只与给定奇点的强度有关，而与其具体分布无关。采用用于处理含有奇点的势流边值问题的F

18、MM方法，将问题域逐层次剖分，再对矩阵元素进行形式处理，使其转化为适用于多极展开公式的形式后，递归运用势函数在球坐标系下的多极展开和局域展开表示，可将矩阵中一部分元素与运算集中起来，转化为向量与相对固定的多极系数的计算，而无需再分别生成并存储这些元素。这一稀疏表示被用于快速近似势的远场影响部分，即核函数与远场边界变量乘积的边界积分；近场的影响采用直接计算的方法得到，两者相加就得到了考虑所有源点影响的场点的势。由此，便隐式构造一个与离散边界积分方程相联系的稠密系数矩阵的稀疏表示。这样，由于计算过程中便不再需显式生成，并存储矩阵中的全部的影响系数，多极展开GMRES算法可使占用内存的乘积最小化，进而可使大型稠密矩阵线性方程组的求解具有较低的存储需求。在现代计算机上，内存容量的增长速度低于CPU和硬盘的增长速度。内存是计算机上最“紧俏”的