470 likes | 660 Vues
( 五 ) 代数方程的求解. 5.1 代数方程系统 5.2 直接法 5.3 主要迭代法 5.4 其他迭代方法. 5.1 代数方程系统. 有限差分 ( 体积 ) 离散格式提供一个网格点 ( 单元)的代数方程 , 以线性代数方程为例: P 点和周围邻居点构成计算模板 ( 比差分基架还大) 计算模板 ( 计算分子 ; 解元 SE). 5.1 代数方程系统 : 计算模板. 2D 2 阶模板. 2D 3 阶模板. 3D 2 阶模板. 5.1 代数方程系统 : 整体方程系统.
E N D
(五) 代数方程的求解 • 5.1 代数方程系统 • 5.2 直接法 • 5.3 主要迭代法 • 5.4 其他迭代方法
5.1 代数方程系统 • 有限差分(体积)离散格式提供一个网格点(单元)的代数方程, 以线性代数方程为例: • P点和周围邻居点构成计算模板(比差分基架还大) • 计算模板(计算分子;解元SE)
5.1 代数方程系统: 计算模板 2D 2阶模板 2D 3阶模板 3D 2阶模板
5.1 代数方程系统: 整体方程系统 • 流场中每一点都有一个方程(小组), 整个计算域就有一个大型稀疏方程系统
5.1 代数方程系统: 系数矩阵的存储 • 只存储非零的对角元素 • 2维5点格式: 5 Ni *Nj • 3维7点格式: 7 Ni *Nj*Nk • Al,l-Nj=W • Al,l-1 =S • Al,l =P • Al,l+1 =N • Al,l+Nj=E
5.2 直接法 5.2.1 Gauss elimination 5.2.2 LU decomposition 5.2.3 Tridiagonal system 5.2.4 Cyclic reduction
5.2.1 Gauss Elimination from forward elimination By backward substitution, we have Require O(n3/3) arithmetic operation Backward substitution O(n2/2) Pivoting Rarely used in CFD
5.2.2 LU decomposition where let then Require O(2n2) arithmetic operation Basis of other iterative methods
5.2.3 Tridiagonal system (TDMA) elimination: * * * Gives upper bi-diagonal matrix. By backward substitution, we get *
5.2.3 Tridiagonal system (cont) • 计算量 O (n) • 周期三对角方程组 • 三对角方程组的并行化解法 • cyclic reduction, recursive doubling, SPP… • 五对角方程组(类似三对角)
5.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 methods
5.3.1 基本概念 Matrix A is sparse 设n次迭代的近似解为 , 不满足上述方程,带入上述方程后有残量 : 迭代误差 迭代解的收敛: 实际计算中:
5.3.2 收敛性 M称为迭代矩阵 • Consider an iterative scheme for a linear system 上两式相减 或
5.3.2 收敛性(续) 设特征向量完备,则 趋于零的充要条件: is the largest eigenvalue 迭代次数:
5.3.3 一些基本迭代方法 • Jacobi method: Converge slow • Gauss-Seidel Method: 2 times as fast as Jacobi • Successive Over-relaxation (SOR if w>1): Useful for solving linear systems occurring in certain PDE’s For positive definite matrix, the SOR converges for
GS迭代法的应用:LU-SGS 奇次迭代步从左下角开始,偶次迭代步就从右上角开始
5.3.4 不完全LU 分解方法 (ILU)在PDE中的应用:SIP方法 • LU method是通用方法,但没有利用原矩阵的稀疏性质; • ILU: 非精确分解,i.e. M=LU =A+N; • 在ILU中,如果迭代矩阵M尽量接近原矩阵A,则收敛快. • ILU method for CFD is Strongly Implicit Procedure (SIP),by Stone . 专用的2D五点格式: M=A+N U L N含有 两个对角线的非零元素,而在 A却为零. M中的元素由矩阵相乘得出:M=LU
Standard ILU: 收敛慢!
Stone (1968):SIP • N在7条对角线都可以有元素 • N和向量φ相的结果尽量接近零 N* φ : 要求:
SIP: (cont) • 带入 (5.39),并等于(5.38),可以得到N的所有元素,并令M=A+N,可得到SIP的LU. • (5.40)仅对PDE的5点离散格式有效。 • SIP求解用更新变量: • SIP求解由L-sweep和U-sweep组成 • 收敛所用迭代次数少,但计算L和U的工作量大,总体效率较高 • 3D 七对角线和2D 九对角线(九点格式)的程序见Peric书附件。
5.3.5 ADI 和其他分裂方法 主要解对多维抛物型方程,也可以解拟时间的抛物型方程-> 椭圆形方程 2D抛物型方程 Crank-Nicholson Discretization where
改写成 The last term is proportion to and can be neglected. 只需求解两个坐标坐标方向的三对角线方程。 2D无条件稳定。3D有条件稳定。特殊形式可以无条件稳定。 增量形式ADI称为 approximate factorization (AF)。 优点: 收敛性快, 计算量不大,缺点:中间变量的边界条件不知道。
5.3.6 Conjugate gradient methods • 线性代数方程和极小化: • 对于对称正定矩阵A,求解 • 共轭 : 等价于找到 , 使F极小化:
5.3.6 Conjugate gradient methods (cont) • 最速下降法:收敛慢,搜索方向可能重复 • 共轭梯度法:新的搜索方向要求和过去所有的搜索方向共轭 n*n矩阵,n次搜索就可以收敛 CG的收敛速度依赖于A的条件数 CFD问题的条件数~ Ni**2 • 改进(其实对所有方法都有效): 预处理
对称正定矩阵方程的 Conjugate gradient method (Golub and van Loan, Matrix computation, 1990) M=C-1, C为pre-conditioning matrix. The choice of M is incomplete cholesky LU
非对称矩阵方程的 Bi-conjugate gradient method • CG 方法只适用于对称系统(如Poisson方程) • 把非对称转化为对称: 第一个方程:原始方程 第二个方程:转置方程,与解无关。 When preconditioned CG method is applied to above system , the following bi-conjugate gradient method results:
适用于非对称 矩阵的Bi-conjugate gradient 算法如下: 2倍于CG的计算量,相同的收敛速度,鲁棒性好
其他解法 • CGSTAB (稳定化的CG) • √GMRES (Saad and Shultz, 1986)
5.3.8 Multigrid methods • 大多数迭代法在细网格上可以很快消除误差的高频分量,但低频分量相当顽固。可以在粗网格上消除这些低频分量。
典型V循环式多重网格法的粗网格、限制和插值算子典型V循环式多重网格法的粗网格、限制和插值算子
两级线性多重网格法步骤 多级多重网格法:继续向更粗的网格限制,直到无更粗的网格为止。在最粗网格上精确求解修正方程。
限制和插值算子: 对于eq (5.63) 1/2*eq(i-1)+ ½*eq(i+1) +eq(i) results in:
Comparison of count for convergence • On 2D Poisson equation, k*k grid, n=k2, unknown Method Cost Gaussian elimination O(n3) GS O(n2logn) CG O(n1.5) FFT/Cyclic reduction O(nlogn) Multigrid O(n) optimal
选择solver • MG+SIP or MG +GS > ICCG > SIP >ADI> GS • GMRES+MG • 没有MG时, ICCG>SIP >ADI >GS
5.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 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.
Sequential solution: Under-Relaxation On the nth iteration the equation for generic variable is Patankar 1980对SIMPLE采用, 稳定求解,但可能降低收敛率 时间相关法就是一种松驰法。
5.4.2延迟修正办法deferred-correction approaches 对于高阶差分离散,如果左端项用高阶差分,则计算复杂 如果左端项只保留相邻点的项,远邻点移到右端,则计算可能发散 为克服上述困难,可用延迟修正: 高阶差分移到右端,同时在左右两端加仅涉及相邻点的低阶差分: 用途:可以处理隐式高阶差分、交叉导数项、非正交相等。但不能处理高阶导数项。
第5次课阅读提示 • 傅《计算流体力学》第5章全部。 • Peric《书》Chapter 5 全部。
第五次课后作业 • 实践Peric《书》附的代数方程求解程序 • (待具体 化)