有限差分方法

第三章有限差分方法§3.1线性微分方程3－1第三章有限差分方法求解常微分方程和偏微分方程是物理学中最常见的问题，而用数值计算方法求方程的解已经是一个发展得相当成熟的领域，这样的方法主要有有限差分和有限元方法。除了物理学的问题以外，数值方法在工程和气象等众多研究领域中也得到广泛应用，由于数值求解微分方程的课题众多，本章中只能涉及针对物理学中一些重要内容所采用的主要普适的方法。对一些标准的微分方程形式，人们已经开发了相应的程序库和软件包。§3.1线性微分方程3.1.1微分方程的分类3.1.1.1偏微分方程的类型微分方程是表达物理量、它对其变量的导数以及变量之间的一个关系。只有一个变量的情形是常微分方程，多个变量下是偏微分方程。如有两个独立变量时，二阶偏微分方程的一般形式是，222220abcdefgxxyyxyφφφφφφ∂∂∂∂∂++++++=∂∂∂∂∂∂，(3.1.1.1-1)其中的系数,,,,,,abcdefg可以是独立变量x和y的函数。如果24bac，称为椭圆形偏微分方程，24bac=是抛物线型，24bac双曲线型。例如，波动方程2222210xvtφφ∂∂−=∂∂，(3.1.1.1-2)是双曲线型，而二维的Poisson方程，()2222,xyxyφφρ∂∂+=−∂∂，(3.1.1.1-3)是椭圆形。抛物线型方程的例子是扩散方程，(),DSxttxxφφ∂∂∂⎛⎞=+⎜⎟∂∂∂⎝⎠。(3.1.1.1-4)在数值计算中，方程类型之间的差别不是那么重要。在（3.1.1.1-1）式中，我们假定了系数,,,,,,abcdefg不包含φ及其高阶导数（否则方程就不是二阶的），这是线性的偏微分方程，否则就是非线性的偏微分方程。在很大的程度上，方程的数值解法要视方程的阶数和是否线性而定。3.1.1.2初始值和边界值微分方程还可以用初始值问题和边界值问题进行分类，对于数值计算来说，这两类之间的差别最为重要。这是因为，初始值问题中，我们是在起始端给定函数值，然后计算在后续各点的函数值。而在边界值问题中，我们是给函数加了限第三章有限差分方法§3.1线性微分方程3－2制条件，要求它在始端和终端必须是一定的值。这些问题中，物理学中的时间变量和空间变量是等同看待的，它们的不同是物理含义不同，但用数学的数值方法处理它们时，它们的地位是等价的，因此，初始值问题不一定就必定是含时的，而边界值就是指空间问题。我们之所以考虑初始值问题和边界值问题的差异，主要是它们对解所加的限制条件的区别。某些偏微分方程是混合型的，即对某个变量是初始值型的，而对另外的变量则是边界值型的。例如，对于谐振子的运动，2220ddtφωφ+=。(3.1.1.2-1)这里只有一个变量，是时间t，我们要找在时间段[]0,Ntt的解。因为它是一个二阶微分方程，需要给定两个初始条件以确定解中的常数。如果指定0tt=时的φ和一阶导数'φ的话，这就是初始值问题。如果给出的两个条件是0tt=和Ntt=时的φ值的话，它就是边界值问题。3.1.2初始值问题3.1.2.1有限差分在有限差分法解微分方程的过程中，独立变量是取离散值的，我们将连续变量的微分方程化成离散变量网格中的差分方程，以初始值或边界值确定每个网格点上的函数值。其具体做法是，首先将微分方程离散化，然后求解差分方程组。对于一个变量的常微分方程，如以步长()hxbaN=Δ=−将x轴上的区间[],ab等间距分割成N等分，则在第n个格点可以构造差值，1nnφφ+−，1nnφφ−−和11nnφφ+−−，分别称为一阶向前、向后和中心差分。由于()()()()2331'''2!3!nnnnnhhhxxxφφφφφ+=++++；(3.1.2.1-1)()()()()2331'''2!3!nnnnnhhhxxxφφφφφ−=−+−+，(3.1.2.1-2)两式相减和相加后分别得，()1111'2nnnnnnnxhhhφφφφφφφ−++−−−−；(3.1.2.1-3)()1122''nnnnxhφφφφ+−+−。(3.1.2.1-4)将微分的差分表示代换到微分方程中即可得到差分方程。对于偏微分方程，我们需要在构造网格面（二维、2个变量）或格点阵（三维），其构造方式与变量的坐标系或边界形状有关。如对于圆形边界，自然用极坐标系下的等角度和等径距比较好，对于方形边界，则采用等边长正方形分割，若两个变量分别是时间和坐标的话，则等时间间隔和等空间间距分割相当于二维平面上的矩形分割。因此，对于函数(),xyφ，（3.1.2.1-1）－（3.1.2.1-4）式中的第三章有限差分方法§3.1线性微分方程3－3导数应该换成对x和对y的偏微分，如()()221,1,,,1,1,2222221,1,,1,1,22,4mnmnmnmnmnmnmnmnmnmnmnxyxyhhhφφφφφφφφφφφφφφ+−+−+−+−+−+−∂∂∇=+≈+∂∂=+++−。(3.1.2.1-5)3.1.2.2Euler法求解初始值问题的最简单的方法是Euler折线法，即简单地用（3.1.2.1-3）－（3.1.2.1-4）式代换进微分方程即可，虽然它的精度有限，因而很少使用，但这个方法反映了初始值问题的基本思想。例如，对于常微分方程()',fxφφ=，Euler折线公式是，()1,nnnnnnhfxhfφφφφ+≈+=+，(3.1.2.2-1)它也可以被看成是简单地将微分方程两端进行积分后得到的，()()11,,nnxnnnnnxdxfxhfxφφφφφ++=+≈+∫。(3.1.2.2-2)该式中的f也可以用1nx+点处的值，得到后退的Euler公式，()()1111,,nnxnnnnnxdxfxhfxφφφφφ++++=+≈+∫。(3.1.2.2-3)（3.1.2.2-2）与（3.1.2.2-3）式有重要的差别：因为（3.1.2.2-2）式中只有左边有1nφ+，称它为显式的，即向前差分Euler法。而（3.1.2.2-3）式左右两边均有1nφ+，称它是隐式的，即向后差分Euler法，通常必须用迭代的方法求解关于1nφ+的非线性方程。隐式法尽管效率很低，但它比显式法稳定。例如，对于谐振子的方程（3.1.1.2-1）式，可得()221120nnnhφφωφ+−+−−=。(3.1.2.2-4)因此由1nφ−和nφ可以算出1nφ+。对于它的初始值，我们给定()001tφφ≡==，以及()'00tφ==。由于第二个条件，有101φφ≈=，将此式代入到（3.1.2.2-4）式后就得到以后各个时刻的23,,φφ。对于扩散方程（3.1.1.1-4）式，设时间分割间隔为τ，空间间隔为h，记(),,mnmnxtφφ=，离散化后得()2,1,1,1,,,2mnmnmnmnmnmnDhSφφτφφφτ++−⎡⎤=++−+⎣⎦。(3.1.2.2-5)求解该方程时需要知道初始值(),0xφ，即对每个m值都知道,0mφ，因此可以由方程陆续地求出,1,2,,mmφφ。值得注意的是，当22Dhτ时，该算法是不稳定的。在数字计算时，误差主要来源于两个方面。一是舍入误差，初始值问题中的误差会逐渐积累，特别是对于小h或大N值时，每步中的舍入误差累计后不再是可以忽略不计的。更重要的误差是截断误差，这是因为我们在（3.1.2.1-3）式中第三章有限差分方法§3.1线性微分方程3－4用两点间的差值近似了导数，忽略了2h以上的高阶项。通常当0h→时，这个截断误差是很小的，但是'φ的变化值或''φ值较大时可能导致计算结果偏离真值，这时我们说结果是不稳定的。控制计算误差使其不能恶性增长是计算方法的稳定性理论所考虑的中心问题。3.1.2.3多步法为了提高计算的准确性，应该在近似的公式中包含2h项。在（3.1.2.2-1）式中，1nφ+只依赖于nφ，这是单步法。将f作插值()()11nnnnxxxxfffhh−−−−≈−，(3.1.2.3-1)代入到积分表达式后得Adams两步法公式，113122nnnnhffφφ+−⎛⎞=+−⎜⎟⎝⎠。(3.1.2.3-2)更一般地在多步法中，第nk+点的值是之前的k个点的函数，()0111111,;,;;,nknnknknnnnnknkaaahFxxxφφφφφφφ++−+−++++=++++，(3.1.2.3-3)其中01,,kaa−是系数，F是某一已知的函数。当它用前面k个格点的值表示出来时是线性多步法()110,;,;;,knnnnnknkiniiFxxxfφφφβ+++++==∑。(3.1.2.3-4)如果0kβ≠，差分方程式（3.1.2.2-2）式右边也包含有nkφ+项，这时方程式是隐式的。3.1.2.4蛙跳法Euler法中，截断误差是()Oh。采用更为对称的中心差分（即蛙跳）法可使截断误差降为()2Oh。由（3.1.2.1-3）式，()()211'2nnnxOhhφφφ+−−=+，(3.1.2.4-1)可得()112,nnnnhfxφφφ+−=+。(3.1.2.4-2)3.1.2.5预报－校正法对于常微分方程还有其它若干高精度的方法。预报－校正法也是有限差分法，其最简单形式包含有为“预报”的向前Euler法和随后对预报结果进行修正补值的“校正”方法。它与Verlet算法结合起来就是分子动力学中积分求解运动方程的最常用算法。第一步是预报步，由显式Euler公式给出φ在1nx+处的第一个预报值*1nφ+，第三章有限差分方法§3.1线性微分方程3－5()*1,nnnnhfxφφφ+=+。(3.1.2.5-1)第二步是用隐式Euler公式给出*1nφ+进行修正，()*111,nnnnhfxφφφ+++=+。(3.1.2.5-2)如果反复多次使用第二步，则该方法就变成了迭代计算，它并不改变有限差分法的阶数，称为迭代预报－校正法：()()*1*111*11,,nnnnnnnnnnhfxhfxrepeatφφφφφφφφ++++++=+⎧=+⎪⎨=⎪⎩。(3.1.2.5-3)3.1.2.6Crank-Nicholson法Crank-Nicholson方法是一种二阶有限差分法，它是对nx和1nx+的f函数求平均，即将差分中一点的显式表示用显式和隐式的线性组合表示，()112nnnnhffφφ++=++，(3.1.2.6-1)这一方法相当于一半x间隔在nx处确定的导数上，另一半在1nx+处确定的导数上，相当于梯形法则。如对于扩散方程（3.1.1.1-4）式，取扩散系数D为常数，源0S=，通常的差分方程是，(),1,1,1,,22mnmnmnmnmnDhτφφφφφ++−−=+−，(3.1.2.6-2)它是显式的，而按照Crank-Nicholson方法，差分方程则变为，()()1,1,,1,11,1,1,1,22222mnmnmnmnmnmnmnmnDhφφφφφφτφφ+−++−+++⎧⎫+−+−⎪⎪−=+⎨⎬⎪⎪⎩⎭。(3.1.2.6-3)Crank-Nicholson方法也可以变为预报－校正法，其中，*1nφ+的第一个预报值由显式Euler公式给出，第二个预报值可利用Crank-Nicholson步的平均化得到，()()()*1*111*11,,,2nnnnnnnnnnnnhfxhfxfxrepeatφφφφφφφφφ++++++=+⎧⎡⎤=++⎪⎣⎦⎨⎪=⎩。(3.1.2.6-4)3.1.2.7Runge-Kutta法由显式Euler方法和Crank-Nicholson方法组合而成的迭代预报－校正法即为Runge-Kutta方法。或者，也可将上述的Crank-Nicholson法看成是广义Runge-Kutta方法的一种特殊二阶情形，事实上，所有以Euler法为基础的方法都归类为推广第三章有限差分方法§3.1线性微分方程3－6的Runge-Kutta方法。Crank-Nicholson法的二阶Runge-Kutta公式（3.1.2.6-4）式可以改写为，()()()1121212,,nnnnnnhFFFfxFfhFxhφφφφ+⎧=++⎪⎪=⎨⎪=++⎪⎩，(3.1.2.6-4)函数F的