本帖最后由 nagami 于 2016-1-18 19:07 编辑
数值方法(一)：经典数值算法整理 文/Nagami 2016-1-17
前言 小编整理了个人觉得非常重要的几种经典数值方法，希望这些东西能起到抛砖引玉的作用。无论对于计算化学、计算物理还是计算数学，这些都是基本的砖头。觉得计算比纯理论来得有趣刺激很多，难道不是么！ 参考的书籍如下： [1]Numerical.Recipes (C++)(3rd)(2007) [2]Numerical Recipes in Fortran 90 (2nd 1996) [3][Jaan_Kiusalaas] Numerical Methods in Engineering 这篇帖子可认为是对Numerical Recipes系列书籍的部分算法做简单概括，但排序、随机数等不涉及。最后小编给了一个实例，用这些算法代码求解一个二元线性偏微分方程组，这是某网友在实际科研中出现的问题。
1. LU分解、Cholesky分解 LU分解的计算量近似2n^3/3,一次分解可适用于任意右端项求解，可认为是Gauss消元的存储策略。故LU分解算法可应用于矩阵求逆、矩阵与矩阵逆相乘计算、行列式计算等。部分选主元的LU算法，又能保证实际使用时的数值稳定性，可将LU分解算法作为一般线性方程组求解的常用方法。对于奇异性不是很强的问题，LU分解也能很好的适用，配合不动点线性迭代改进，效果会很好。由于选主元打乱了矩阵行列的顺序，如此在LU分解上做rank-1（uv‘）更新分解，算法难度上会大大增加。 Cholesky分解由于考虑了矩阵的对称性，计算量原则上是LU分解的一半。同时正定性是保证传统Cholesky算法能顺利进行所必须的，故如果Cholesky分解失败，有必要怀疑该矩阵的正定性缺失。Cholesky分解有rank-1更新算法，见wiki
2. 三对角与带状矩阵分解 三对角矩阵问题，无论是求解，还是特征值计算，在实际中都出现得极为平凡，比如三次样条插值，一维薛定谔方程特征值求解，微分方程的有限差分离散等。它的一般推广形式称为带状矩阵，考虑到零值的分布，采取了紧凑存储，对应的LU分解算法是一种很好的选择，由于存储了分解，同样可对多右端项求解。 一般的，三对角矩阵可用所谓的追赶法来求解，在分解的同时顺带完成回代。寻求三对角矩阵求解的并行算法也是颇为有意义的。
3. QR分解、SVD算法 QR分解的计算量近似4n^3/3，是LU分解的两倍，这是由Hourseholder变换产生的，所以一般不用于线性方程求解，但有些情况下QR分解也是一种可选的方案。理由是QR分解不选主元，也能适当保证分解的稳定性，QR分解的rank-1更新算法，只需O(n^2)就可获得对应矩阵QR分解，这可应用于多维求根Broyden方法之中。 SVD算法，也称为奇异值分解算法，比QR分解计算量更大。它的优点是能处理其他算法失效的极为病态的矩阵问题。比如在偏微分方程的反演问题中，出现很多内在病态矩阵问题，SVD分解是一种策略，自然，结合正则化技术效果会更好。SVD算法在数据压缩领域表现得也很风骚。
4. 稀疏矩阵存储 一般在网格上离散偏微分方程的数值方法会出现各种形式的稀疏矩阵，比如有限元法、有限体积法等，如何紧凑的存储阶数很大且非零值所占比例又小的矩阵，是一个难题，也是快速求解问题的关键点。目前常用的有两种，一种是行（列）压缩存储；另一种是标记横纵指标存储。前者所占内存少，后者并行效果好。稀疏矩阵求解，有分解算法，但绝大多数考虑迭代算法，比如CG迭代。
5. CG迭代、不动点迭代 不动点迭代比较有名的是超松弛迭代算法，但需要选择合适的超松弛因子才能体现它的优点，使用起来也是颇为麻烦。现在的共轭梯度（CG）系列算法基本取代了它的地位。一般而言，不动点迭代收敛着眼于所构造矩阵的谱半径。CG迭代的收敛率依赖矩阵的条件数，条件数越接近1越好。也正是预条件处理技术的出现，使得CG算法如此普及。
6. Vandermonde矩阵与Toeplitz矩阵的快速求逆技术 在数值方法中，会出现很多与多项式有关的问题，比如多项式插值中求多项式系数、数值积分代数精度优化等问题。这类问题往往与Vandermonde矩阵相关。 Toeplitz矩阵在电磁、信号处理等领域经常出现，比如螺线管自电容计算、DFT计算等。 求解这些特殊矩阵，都有对应的O（n^2）算法。用LU分解也是可行的，尤其是这些矩阵比较病态的时候。 Vandermonde矩阵可利用Lanrange插值多项式的性质，再结合多项式的符号除法，可建立对应的快速求解算法。 Toeplitz矩阵快速求解有很多种，见wiki，一般的，Levinson迭代可作为学习的范本。
7. 多项式插值、有理式插值 插值问题从属于反演问题，本身的病态性质不可忽视。对于多项式插值，一种稳健的算法称为Neville方法，可同址进行计算。 它的缺点是算法本身是针对单点计算设计的，不适合做多数据计算。如果插值问题的病态性不是很高，求解对应多项式系数也是可选的方法，可多值计算。 有理插值是多项式插值的自然推广，Bulirsh和Stoer发现了一种与Nevile方法类似的有理插值算法，值得一看。 进一步的推广，称为BaryCentric Rational Interpolation，这是一种精度可调的有理插值算法。可见如下文献 Berrut, J.-P., and Trefethen, L.N. 2004, “Barycentric Lagrange Interpolation,” SIAM Review vol. 46, pp. 501–517.
8. 三次样条插值 众所周知，多项式插值存在Runger现象，有理插值表面上削弱了这种现象。但根本上，是需要全局平滑的插值算法。三次样条插值是一种平民级的选择。 插值系数建立后，也可应用于数值积分。比如对于非等间距的数值积分，这就是一种不错的选择。
9. 扩展梯形求积、Romberg积分、积分的变量替换技术(Variable Transformation) 扩展梯形积分算法，是一种等间距细分积分法，本身精度不高O(h^2),但在应用中却颇具欢迎，原因在于Euler-Maclaurin公式给出了该公式的渐进阶，发现其是偶次提升的。对此，有两种应用上的策略可以被考虑到： 第一种是，做合适的外推，使得大量偶次阶被消去，收敛速度很快，比如Simpson积分、Romberg积分等，此法对于光滑非震荡函数效果不错。 第二种是，做变量替换，由于Euler-Maclaurin公式的截尾误差与函数的端点性质有极大关系，选择一些特殊的变换，收敛效果很好。此法也可处理奇异积分。比如TANH、DErule、IMT等变换。 同时注意一种与三角函数相关的积分方法：Censhaw-Curtis Rule。
10. 经典权Gauss积分、任意权Gauss积分问题 Gautschi等人大力发展了基于正交多项式的数值积分方法，利用对称三对角矩阵的特征值求解，可在O（n^2）内计算出任意非经典权的nodes和weights。比如量化的Rys积分，Multexp积分，以及统计物理中的Fermi-Dirac积分等。这些算法本身并不复杂，但思想却非常深刻。 如果获得了该正交多项式的递推公式，同时根的近似公式也知道。基于牛顿迭代改进的求权算法非常不错，其比一般的算法速度上可快6-10倍。现在经典权的Gauss积分求权程序非常普及了。
11. 振荡积分、级数加速收敛算法 对于一些振动比较剧烈的函数考虑数值积分，按照经典的做法，不仅收敛速度很慢，同时计算量也非常大。 这种情况，自适应数值积分是一种可选的方法。它的思路是自迭代思想，比如：在区间[a,b]上考虑两种不同精度的求积公式，做差继而判断是否可接受，不接受则对半细分，在[a,(a+b)/2]，[(a+b)/2,b]上继续上述迭代。 振荡积分的另一种考虑是把区域做无穷细分（对于无限域积分），使得在每个子区域上使用一般算法就可获得足够精度，继而原问题转化成了级数求和问题。这样如何加速级数收敛（或者说序列加速），就自然产生了。这方面有Euler变换，Levin变换等。不得不说Levin变换是一种值得花时间去了解的序列加速算法。
12. 多维积分与蒙特卡罗积分 对于多维积分问题，可考虑的因素会很多，应该针对不同问题采取不同的策略： 1）多维积分是否可降维处理，比如转一维积分 2）被积函数是否奇异，光滑性如何，注意到Gauss积分计算量少，如果在精度上可接受那再好不过，但对于奇异问题其效果或许会很差，除非做对应的奇异权Gauss积分。 3）被积区域是否规则，可否转换到合适的多重积分处理。对于不规则区域，蒙特卡罗积分是一种不错的选择，虽然其收敛阶不高。 4）边界奇异的数值积分问题，不妨采用变量替换技术试一试。 等等
13. 切比雪夫多项式逼近、Pade逼近 切比雪夫多项式逼近是一种最佳上确界范数逼近，整体误差界在同阶多项式最佳，对于一个很难计算的函数，建立逼近的切比雪夫多项式，可方便处理导数计算和积分计算。 Pade逼近是把多项式改造成有理式，通过其各阶导数相等建立等式约束，最终转换到求解对应的线性方程组，该问题可能会很病态。如果逼近函数相对平坦，Pade逼近在“远区”的效果比多项式会来得更优，显然Taylor展开的多项式的增长速度较快，它只具备局部逼近。
14. Gamma函数、Beta函数及其导出函数 以下15-19特殊函数计算，以后有时间详细描述 15. 不完全Gamma函数、不完全Beta函数、Gauss误差函数、指数积分 16. Bessel柱函数系列（整数阶、分数阶、球Bessel函数，Airy函数） 17. 球谐函数求值 18. Fresnel积分、Dawson积分与采样定理 19. 椭圆积分与椭圆函数 20. 二分法求根、Ridder方法、一维Newton法 二分法是一维求根中利用界定来细分逼近的方法，线性收敛速度，不过稳健性很高。区间的连通性和函数的连续性在原理上确保最终能逼近某个零点。Ridder方法是一种不依赖导数的求根算法，超线性收敛，据说稳健性也不错，但不具备多维推广性。在这方面，所谓的Open Method，比如Newton方法非常著名，具备局部的二次收敛率，依赖导数。遗憾的是Newton方法对初值选择极为敏感。所以一些结合性方法也是一种不错的选择，比如Brent求根方法：结合了二分法、抛物线反插值、弦截法。
21. 多项式求根算法：Laguerre方法 Laguerre方法是多项式求根的平民级算法，可针对复系数多项式求根。一些软件包里面，可能有更加优秀的取代算法。 一般数值方法中，Laguerre算法的原理并不十分严格，但在实践中它很成功。 还有一种，是对二次项（ax**2+bx+c）进行迭代改进。
22. 全局收敛策略的高维Newton法与Broyden方法 J. E. Dennis, Robert B. Schnabel-Numerical methods for unconstrained optimization and nonlinear equations-Society for Industrial and Applied Mathematics (1987)
23. 黄金比寻优法、Brent方法 与二分法类似，黄金比寻优法是通过界定极小值所在区间的一种搜索算法。Brent方法是抛物线插值算法。显然，一维问题界定一个极小值需要三个点。而求根只需两个点。
24. 单纯下山法 单纯下山法，一种多维寻优算法，也称为Nelder-Mead算法，是matlab默认的寻优算法。该算法只需要计算函数值，通过对比单纯形各顶点上的函数值，按照算法策略对单纯形进行伸缩反射变换，最终收敛到目标值。就像一只变形虫，触手上带有高度传感器，通过判断，逐渐爬向低谷。
25. Powell方向集方法、共轭梯度法 以后有时间详细描述
26. 拟牛顿法 以后有时间详细描述
27. 解对称矩阵特征值的Jacobi方法与QL算法 Jacobi算法是比较古老的求解对称矩阵特征值算法，在矩阵阶数大于10以后，速度上由于一个影响因子大大减慢。之后QL算法出来，其被取代，但Jacobi rotation的特性仍然有用武之地，比如QR的rank-1 update。 QL算法，或者对应的QR算法（依据数值稳定性做对应的选择），目前无论在理论还是在程序设计上都极为漂亮，Lapack中tred2和tqli可作为平民级的对称矩阵特征值求解器。其中tred2是利用Househodler变换转到对称三对角矩阵，然后tqli再进行隐式移位QL迭代，每次迭代收敛甚至达到3阶以上。特征值的QR算法被誉为20世纪最伟大的十大算法之一。
28. 快速Fourier算法 以后有时间详细描述
29. 常微分方程：自适应步长控制的Runge-Kutta方法 RK方法是比较有名的常微分方程求解的单步长积分方法。对于非刚性问题，RK方法效果不错。而对于刚性问题，必须考虑隐式格式。在实践中，自适应步长控制是需要的，单步长RK方法如何进行误差估计？必须要加以考虑，一种改进的方法称为Embedded Runge-Kutta方法。嵌入式RK方法需要准备两套同nodes，不同weights的步进公式，由于它们的精度不同，做差可作为误差估计使用，进一步分析可获得步长控制数据，比如：多大的步长可在接受范围内、下次步进时步长又设置为多少合适等。 要想嵌入式RK获得等间距的数据，需要考虑用插值去做，有相关文献。
30. 常微分方程：Richardson外推的Bulirsch-Stoer方法 将Romberg积分的思想推广到ODE求解中，就是所谓的Bulirsch-Stoer方法。BS-ODE求解通过建立偶次阶步进公式（修正中点公式），利用外推获得对应步长的积分值。如果ODEs右端项极为光滑，BS-ODE方法会比RK-ODE方法效果更佳。BS-ODE的原理是很精细，算法上借鉴了通信原理思想。比如如何进行外推、多项式外推or有理式外推哪个好、如何判断外推失败等。
31. 常微分方程两点边值问题：打靶法与迭代松弛法 打靶法是把ODE求解器与多维求根结合来求解ODE两点边值问题。 迭代松弛法是通过有限差分离散，建立迭代格式来求解。
32. 积分方程：Nystrom方法与偏微分方程反演问题 积分核越光滑，该问题病态性就会增加。但在实际中，比如库伦势，它的奇异性反而削弱了这种病态性。这方面的理论需要泛函分析的理论。
33. 偏微分方程 一）偏微分方程定性方面 1.Hilbert Space Method
关于微分算子的一种自然的表述是广义函数下的解释。因而PDEs被理解为广义函数空间中的一个元，在检验函数空间的归纳极限拓扑下会加宽微分运算的使用条件。由于大量的函数空间的对偶空间可以被嵌入在广义函数空间里面，PDEs问题可被理解为一种对偶形式（有限元法的变分格式就是对偶式）。Hilbert space 一方面是自对偶的，另一方面有一个自然的内积结构与物理的能量观念对应。大量的力学及物理问题可由变分不等方程描述。热学，弹性，粘弹性，塑性，各种摩擦问题，Bingham流体，电磁，etc
当然一般不能在广义函数空间下面讨论这类问题，一个重要的困难在于乘法运算。局部可积函数对应的广义函数，加入范数结构，克服了这个问题，一个著名的函数空间称为Sobolev space。确定在什么样的空间下面讨论PDEs的解，这是首先面对的问题，继而讨论非齐次，边界的几何条件等其他问题。
2.Weak Convergence Method
脱离Hilbert space的良好结构，这是Banach space中求解PDEs的另一种方法。首先无限维Banach space中单位球紧性缺失，要考虑各种弱拓扑。通常是选择使得连续线性泛函族（范拓扑下）保持连续的最弱拓扑。一系列的讨论，比如紧性，可分性就可展开。再在其上建立PDEs理论
3.Fixed Point Method
为了说明某些存在性问题，除了紧性，fixed point 的观点另辟蹊径。一般有metrix fixed point theorem，但最有影响的还是topological fixed point theorem。拓扑化已经是一个趋势，与拓扑相关的结果一般具有可遗传性，就是在同胚下保持某种不变性，当然这里说的是不动点存在性。各种延拓的技术在此也表现得更为自然。PDEs问题从另一种角度所解释，我们可以在normed space下考虑这类问题，利用代数拓扑工具，直接分析算子方程。
4.Lie Group Method
李群分析与PDEs，有一些书讲这个，比如“瑞典人Nail H. Ibragimov编写的李群分析与偏微分方程关系式之手册”，以代数角度来分析和求解一些特殊的PDEs类型。
二）偏微分方程数值方面：FDM、FEM、FVM、Relaxation方法、Multigrid方法、谱方法etc 以后有时间详细描述， 见一个具体例子： 这是一个二元线性偏微分方程组的求解，方程及公式见附件图片，由于边界条件很特殊，可考虑利用三角函数对角度部分做有限维逼近，将PDEs系统转化到ODEs系统。程序中角度部分采用了5阶逼近计算。最后转化到XY坐标输出所需的数据，用matlab的三角化网格画图。
| 
发表于 Post on 2016-1-17 19:37:31
收藏 Add to favorites
