第三章椭圆方程迭代法介绍

第三章椭圆型问题的差分法§3－1流体力学中的椭圆型问题·无旋流场中速度势02（LaplaceEq.）·二维不可压定常流动，利用涡－流函数表示：方程PoissonVt22·不可压分离流问题中，扰动压力场：02p·定常的N－S方程求解问题·在网格自动生成中，求解椭圆型方程的网格生成方法由于椭圆型方程的数学性质：求解域内部任何一点的解函数依赖于所有边界上的边界条件，因此从数值计算方法来看，就不能从一部分边界起步进行推进计算到另外的边界，这与发展方程的求解方法有很大的差别，椭圆型方程的数值求解方法，只能是在整个流场中进行迭代计算来求解。§3－2椭圆型问题的迭代法求解（一）迭代法的基本概念例：方程2（Poisson方程）二维2222yx差分离散jijijijijijijiyx,21,1,1,2,1,,1)(2)(2………………..(*)写成矩阵形式代数方程组为：BA……………………………………..(1)其中AJIji,,2,11,1一般地，对于线性方程组有BA，欲求未知函数的解矢量若A为非奇异矩阵，即：1A，则BA1由于A是个阶数甚大的矩阵（非三对角），直接求解，或利用Gauss消去法求逆矩阵，计算量及所需计算机的内存都将十分巨大，所以在实际计算中不希望采用直接法求解。迭代法的基本思想是：定义一个序列)(k，当k时，BAk1)(，从而得到方程(1)的解。迭代法设法给出)()2()1()(,,,,rkkkkkBAF的迭代关系。（通常为计算方便，迭代法采取1r，使之简单）)1()(,,kkkBAF若kF（即迭代关系式）与迭代步k无关，则称为平稳迭代；若kF是)1(k的线性函数关系，则称为线性迭代。例如最简单的线性迭代关系可设为：)()1()()(kkkkVH………………………(2)若迭代是有效的，则)()(kkVH即)(1)(1kkVBAHBA…………………………(3)BMBAHIVkMkkk)(1)()()()(即BMHkkk)()()(且1)()(AHIMk或IMAH研究迭代的收敛性：引入误差：BAEkk1)()(而由(2)-(3)得：)(1)1()(1)(BAHBAkkk即)1()()(kkkEHE或有递推关系式：)0()2()1()()(EHHHHEHEHEkkkkk个由于)0(E是初始解与精确解的误差，应是一个有界的任意函数，故迭代矩阵H应具有：当k时，0limZHHHkk，Z为任意的有界向量函数。可以证明：（参阅“偏微分方程的有限差分方法”P239）对于任意的向量Z，0)1()1()(ZHHHkk的充分必要条件是H的所有的特征值i的绝对值（即谱半径）都小于1。推论当k很大时，)(~)()1(HEEkkiiHmax)(所以若1~，则迭代法的收敛速率很慢。二、几种迭代法介绍1．Jacobi迭代(简单点迭代)由方程BA将矩阵分解为：A=L+D+UL：主对角线以下的元素ija（ij时等于A，其余为零）D：主对角线元素U：主对角线以上的元素ija（ij时等于A，其余为零）BUDL)(BUDLkkk)1()()1(BDULDkk1)1(1)()()(1DLDH,1DM,BDV1H，M可以验证满足迭代有效性条件，即IMAH2、Gauss-Seidel点迭代类似1但是BUDLkkk)1()()(BDLUWDLkk1)1(1)()()(在实际计算中L中（ij）只要遵循已有新值时，用新值，没有新值时用旧值，即为G-S。*往返扫描的Gauss-Seidel迭代，即step1:BUDLkkk)1()()(step2:BUDLkkk)1()1()(3、SOR（逐点松弛迭代）step1.用G－S迭代法求中间值，即BUDLkkk)1()()(…………………………………(a)step2.)1((*))()1(kk….…………………………….(b)消去中间结果(*)即)1((*))()1()(kkDDDDb将(a)代入)1()1()()()1(][kkkkDULBDBUDDUBLDkkkk)1()1()1()(])1[()1()(BLDUDLDkk1)1(1)()(])1[()(其中为松弛因此，10为亚松弛，21时为超松弛。4、线迭代和线松弛迭代将A分解为ULDA但D保留主对角元素在D中，L，U则仍为余下元素的上三角与下三角矩阵则BUDLkkk)1()()1(导出线迭代而BUDLkkk)1()()(导出线D－S迭代*往返扫描的G－S线迭代（线松弛迭代）而)1((*))()1((*))()1(kkkkDDBUDL导出松弛迭代三、迭代法的收敛性及松弛因子的选择1．迭代法收敛的几个充分条件对于方程BA①若矩阵A满足强对角优势条件，则Jacobi迭代和G－S迭代均收敛②若矩阵A满足对角占优条件，且矩阵A为不可约矩阵，则Jacobi迭代和G－S迭代均收敛③若矩阵A是对称正定矩阵，则G－S迭代收敛。④若),2,1(0Niaii且有12)(12iiNijjijaa则Jacobi，G–S迭代收敛⑤若对于BAJacobi迭代收敛的，则21的松弛迭代也总是收敛的。只证明①（余略）A为强对角占优，及11Nijjijiiaa将A的每一行元素均用该行的主对角线元素去除，可得到主对角元素为1，且不改变将对角占优的性质，BA，然后将A分解为UDLA，且有：1ijijUL对于Jacobi迭代，BULBDULDkkk)1(1)1(1)()()(即)(ULG利用矩阵的特征值分布定理(Gerschgorin圆盘定理)，可知A的所有特征值均在单位圆内，证毕！2、对于Poisson方程Jacobi迭代矩阵的特征分析结论：1、Jacobi迭代矩阵的特征值为：（参考苏煜诚，吴启光，偏微分方程数值解）1,,2,11,,2,1cos21cos21,lrmslrmsrsx方向总网格数为s+1(0,1,2,…,s)，0,s为边界y方向总网格数为l+1(0,1,2,…,l)，0,l为边界2、逐点松弛迭代法中迭代矩阵的特征值（G－S或SOR）由})1{()()(1UDLDG设)(G的特征值为则结论为：0)1(,,,srsrsr(1)3、SOR方法中松弛因子的最优化迭代)(,,,)(maxMin)(srsroptsrMax＝谱半径由0)1(，使1的充要条件是：11（20）)1(1结论：2*2*2*112)11(2opt（2）rsrs,,*maxJacobi迭代矩阵的谱半径并有2*2*11111)(optopt但是由于实际过程中)]([*JG，未知，所以opt不能预先知晓。4、优选松弛因子的两个近似方法方法1：利用的关系仍为上面之（1）（2）两式步骤①取200用SOR迭代计算若干步，然后用下面的计算近似的)(0ijkjikjiijkjikjiuuuu)1(,)(,)(,)1(,0)(②以0，)(0代入（1）式求20202)0()1(][B③根据（2）2)1(112Bopt即为)1(opt的第一次近似值可以类似求出)1(，)2(，…，直至)(k，)1(k之差小于为止。方法2、令)1()(,)(maxkijkijjikuu由1开始，近似认为)2()1()1()(kkkk直至akkkk)2()1()1()(~取2112aopt5、几种主要的迭代算法的收敛速度比较设而为问题求解域为yx00nh内点共有2)1(n个a.Jacobi迭代：))cos()(cos(21shrhrs21,1,,211coshmaxhsrsrJ所以收敛速度2121lnhRJb.G－S迭代由0)1(当1时，即为G－S，2所以22maxmaxJrsSG2ln2lnhRJSGSGc．SORhGopt21sinh1sinh11111)]([2*2*hR2ln§3－3定常问题的迭代法求解与（伪不定常）时间推进法计算的一致性讨论一、概述例1，定常方程xuxua22(*)采用Jacobi迭代，差分格式用中心差分2)(1)1()(1)(1)(122xuuuxuuakjkjkjkjkj)(21)(4)(1)(1)(1)(1)1(kjkjkjkjkjuuuuxau)2(21)(4)(1)()(1)(1)(1)()1(kjkjkjkjkjkjkjuuuuuxauu2)(22)(1)(12)()1(22)(4xuxxuuxauukjxkjkjkjkj其中可以视为虚拟的时间步长)(2)(2)(222222)(xOxuxxOxuxaOukj222222xuxxuxau即从方程2222xuxuaux出发的FTCS格式，与从方程（*）出发的Jacobia迭代得到相同的差分。或稳定条件：1222xaxxaxtac212222xxxts12'2sc12xa稳定条件即取：2xaRxe例22tkijkijykijxkijkijyxt22221为简单起见令yx]4[2,1,!,,!,!21kijkjikjikjikjikjikijkijt若取41s即412t][41]4[4121,!,,!,!2,1,!,,!,!1kijkjikjikjikjikijkjikjikjikjikjikijkij将k，k+1视为相邻两个迭代步的解，则上式是原方程2的Jacobi迭代，即t可视为一个虚拟时间项（时间相关！）。例30采用线松弛迭代求解椭圆方程(y方向是隐式求解，x方向是G－S迭代)step1:0222(*)1,(*),(*)1,2)1(,1(*),)1(,1yxjijijinjijinji(1)step2:)1(,(*),)(,)1(njijinji(2)由(2)式)1(,)(,(*),)11(1njinjiji