代数方程的求解

1、1(五五) 代数方程的求解代数方程的求解 5.1 代数方程系统 5.2 直接法 5.3 主要迭代法 5.4 其他迭代方法25.1 代数方程系统 有限差分(体积)离散格式提供一个网格点(单元）的代数方程, 以线性代数方程为例： P点和周围邻居点构成计算模板(比差分基架还大） 计算模板(计算分子;解元SE)(1) lPllPPQAA35.1 代数方程系统: 计算模板2D 2阶模板2D 3阶模板3D 2阶模板45.1 代数方程系统: 整体方程系统 流场中每一点都有一个方程(小组), 整个计算域就有一个大型稀疏方程系统向东，向上。从西南角开始，向北，在结构网格上的排序：的排序。其结构依赖于稀疏方阵,

2、:(2) AQA55.1 代数方程系统: 系数矩阵的存储只存储非零的对角元素2维5点格式： 5 Ni *Nj3维7点格式： 7 Ni *Nj*NkAl,l-Nj=WAl,l-1 =SAl,l =PAl,l+1 =NAl,l+Nj=E PQEENNPPSSWWAAAAA65.2 直接法5.2.1 Gauss elimination5.2.2 LU decomposition5.2.3 Tridiagonal system5.2.4 Cyclic reduction75.2.1 Gauss EliminationBy backward substitution, we havefromRequir

3、e O(n3/3) arithmetic operationBackward substitution O(n2/2)Pivoting Rarely used in CFDforward elimination85.2.2 LU decompositionQA 的所有元素）（可求出UL, ALU QLU YU QLY whereletthenRequire O(2n2) arithmetic operationBasis of other iterative methods95.2.3 Tridiagonal system (TDMA)*Gives upper bi-diagonal matr

4、ix. By backward substitution, we getelimination:*105.2.3 Tridiagonal system:块三对角方程组 1*1*11*11*11:onsubstitutiback and:neliminatioiiEiiPiiiPiWiiiEiPiWiPiPiiiEiiPiiWAQAQAAQQAAAAAQAAA115.2.3 Tridiagonal system (cont) 计算量 O (n) 周期三对角方程组 三对角方程组的并行化解法 cyclic reduction, recursive doubling, SPP 五对角方程组（类似三对角

5、） 125.3 迭代法 5.3.1 基本概念 5.3.2 收敛速度 5.3.3 一些基本方法 5.3.4 不完全LU 分解方法 5.3.5 ADI 和其他分裂方法 5.3.6 Conjugate gradient methods 5.3.7 Bi-conjugate gradients,CGSTAB, GMRES 5.3.8 Multigrid methods13迭代误差迭代解的收敛：Matrix A is sparse 设n次迭代的近似解为 , 不满足上述方程，带入上述方程后有残量 ：nn5.3.1 基本概念基本概念0or , 0实际计算中:速度预处理矩阵，加速收敛PPQPA145.3.2

6、收敛性 Consider an iterative scheme for a linear system上两式相减或n 这里M称为迭代矩阵15设特征向量完备，则1is the largest eigenvalue迭代次数：5.3.2 收敛性（续）趋于零的充要条件：1k165.3.2 收敛性：收敛速度0,NAMNMAQA要想收敛快，要求：17Jacobi method:iikikikiAR1nixAxAQRnijkjijijkjijiki,.,2 , 1,111Gauss-Seidel Method:iikikikiAR1nijkjijijkjijikixAxAQR1111Successive

7、Over-relaxation (SOR if w1):Useful for solving linear systems occurring in certain PDEsiikikikiAR1nijkjijijkjijikixAxAQR111For positive definite matrix, the SOR converges for . 20Converge slow2 times as fast as Jacobi5.3.3 一些基本迭代方法18GS 和SOR的一般形式ijijiiijijiipppppppaaiiaaXUXbLXSORUXbLXUXbLXGSULAbAX有且至

8、少对一个收敛条件：,)1 (:,1111119GS迭代法的应用：LU-SGSpppppppppppppppXUXDXDeiXUXLAXbXDXLAXbXDAXbXUDLXXXUDLAbAX*1 .:SweepU:SweepL)(,则奇次迭代步从左下角开始，偶次迭代步就从右上角开始20GS迭代法的应用：线-SGSjinjinjnjinjinjnjifuuuBuuuAybBxaAfyubxua,11,111, 111, 1222222)2()2(GS, ：线定义对于方程21GS迭代法的应用：并行的Red-black225.3.4 不完全LU 分解方法 (ILU)在PDE中的应用：SIP方法LU m

9、ethod是通用方法，但没有利用原矩阵的稀疏性质；ILU: 非精确分解，i.e. M=LU =A+N;在ILU中,如果迭代矩阵M尽量接近原矩阵A，则收敛快.ILU method for CFD is Strongly Implicit Procedure (SIP)，by Stone. N 含有 两个对角线的非零元素，而在 A却为零.M中的元素由矩阵相乘得出：M=LU专用的2D五点格式：LM=A+NU23Standard ILU:收敛慢！24Stone (1968):SIP N在7条对角线都可以有元素 N和向量相的结果尽量接近零N* :要求：25SIP: (cont)带入 (5.39)，并等于

10、（5.38)，可以得到N的所有元素，并令M=A+N,可得到SIP的LU.(5.40)仅对PDE的点离散格式有效。SIP求解用更新变量： SIP求解由L-sweep和U-sweep组成 收敛所用迭代次数少,但计算L和U的工作量大，总体效率较高3D 七对角线和2D 九对角线(九点格式）的程序见Peric书附件。ppAQLU1265.3.5 ADI 和其他分裂方法 主要解对多维抛物型方程，也可以解拟时间的抛物型方程- 椭圆形方程Crank-Nicholson Discretizationwhere2D抛物型方程27改写成The last term is proportion to and can b

11、e neglected.3)( t只需求解两个坐标坐标方向的三对角线方程。2D无条件稳定。3D有条件稳定。特殊形式可以无条件稳定。增量形式ADI称为 approximate factorization (AF)。优点: 收敛性快, 计算量不大，缺点：中间变量的边界条件不知道。285.3.6 Conjugate gradient methods 线性代数方程和极小化： 对于对称正定矩阵A，求解共轭 :0212211021PPPPAAFFPP是共轭的：关于矩阵如果这两个方向意一个方向做极小化，我们只需要针对其中任：平面内极小化假设在等价于找到 , 使F极小化： i295.3.6 Conjugate

12、 gradient methods (cont)最速下降法：收敛慢，搜索方向可能重复共轭梯度法：新的搜索方向要求和过去所有的搜索方向共轭 n*n矩阵，n次搜索就可以收敛 CG的收敛速度依赖于A的条件数 CFD问题的条件数 Ni*2 改进（其实对所有方法都有效): 预处理minmaxCACCQCCACC精心选择共轭梯度法应用于11111(1) 30M=C-1, C为pre-conditioning matrix. The choice of M is incomplete cholesky LU对称正定矩阵方程的 Conjugate gradient method(Golub and van L

13、oan, Matrix computation, 1990) 31非对称矩阵方程的 Bi-conjugate gradient method CG 方法只适用于对称系统(如Poisson方程) 把非对称转化为对称：0AQAT第一个方程：原始方程第二个方程：转置方程，与解无关。When preconditioned CG method is applied to above system , the following bi-conjugate gradient method results:32适用于非对称 矩阵的Bi-conjugate gradient 算法如下:2倍于CG的计算量，相同的

14、收敛速度，鲁棒性好33其他解法 CGSTAB (稳定化的CG) GMRES (Saad and Shultz, 1986)345.3.8 Multigrid methods 大多数迭代法在细网格上可以很快消除误差的高频分量，但低频分量相当顽固。可以在粗网格上消除这些低频分量。10031)(2)2()2(GS, ,1,111, 111, 1222222GGBABeAeBABeAeGfuuuBuuuAybBxaAfyubxuaiiiijinjinjnjinjinjnji，有，对于低频分量，有且对误差放大因子迭代法：点中心差分定义对于方程35典型V循环式多重网格法的粗网格、限制和插值算子36两级线性

15、多重网格法步骤多级多重网格法：继续向更粗的网格限制，直到无更粗的网格为止。在最粗网格上精确求解修正方程。37公式描述:线性方程(3) )(I,)2(.(2) )(R,(1) cfcfoldfnewcffhfcfchfhffffffffhfcffLELLQQL可以得到，将此更新插值到细网格的近似解是设的误差方程该方程逼近于细网格上网格上求解差比较光滑，可以在粗当收敛较慢时，表明残残差为更新近似解细网格上需解方程：38公式描述:非线性方程(6) )(I(5) )(R)(R:(4) foldcfcfcfoldfnewfcfhfcfchfhffhhffhRLLLQLLQLccff网格上的误差方程差比较

16、光滑，可以在粗当收敛较慢时，表明残细网格上残差误差方程细网格上需解方程：39限制和插值算子:对于eq (5.63)1/2*eq(i-1)+*eq(i+1)+eq(i) results in:40Comparison of count for convergence On 2D Poisson equation, k*k grid, n=k2, unknown Gaussian elimination O(n3)GS O(n2logn)CG O(n1.5)FFT/Cyclic reduction O(nlogn)Multigrid O(n)optimal 41选择solver MG+SIP or

17、 MG +GS ICCG SIP ADI GS GMRES+MG 没有MG时， ICCGSIP ADI GS425.4 其他迭代法coupled equations (system of nonlinear equations) Simultaneous approach: All equations are considered part of a single system. Sequential approach: Each equation is solved for its dominant variable, treating the other as known, and one

18、 iterates through the equations until the solution of the coupled system is obtained. Iterations performed on each equation are called inner iteration. In order to obtain a solution which satisfy all equations, the coefficient matrices and source vector must be updated after each cycle ad the process repeated. The cycle are called outer iteration.43Sequential solution: Under-RelaxationOn the nth iteration the equation for generic variable is Patankar 1980对SIM