偏微分方程数值解第二次大作业.docxVIP

偏微分方程数值解——第二次大作业作业要求：零狄氏边界条件；非齐次狄氏边界条件；非齐次其他边界条件。零D氏边界条件采用C-N差分格式，则有所用记号，且令，最终，C-N差分格式可以表示如下：解上述方程使用FST方法，对系数矩阵进行快速对角化处理，加快计算的速度，具体如下：记，则C-N差分格式与下式有同样的时间空间精度，因此可以用下式代替：求解上式的算法步骤如下：对于零D氏边界条件，有其中和的可以利用FST技巧进行计算，如下所示： 根据FST的原理，有由此可以快速求解得到（g为矩阵或者向量）的值。计算时，取源项函数，即这个方程的精确解为具体程序见附录（程序1.1），计算结果如下所示：表一、零D氏边界条件下的FST法计算结果（）N*2012345误差1.51E-057.47E-067.10E-073.73E-064.26E-07时间(s)0.34631621284832.7205176.93984图1a. N与计算误差的双对数图（FST法）图1b. N/20与相对误差关系图程序中，取空间步长；本来 C-N格式对网格比没有要求，但是由于采用了以上算法对其的改造，特别是迭代第一步采用了显式格式，由于其稳定性对时间步长的要求，这里取网格比，也即时间步长。此外计算误差取为近似解与精确解之差的二范数模长，相对误差定义为计算误差与精确解二范数的比值，精确解取为时的函数值。如图1a所示，在N100范围内，计算误差随着N的增大逐渐减小，但毫无规律可循。然而值得庆幸的是，计算误差全部都在10-5量级以下，计算精度还是比较高的。此外，图2b还反应了结果的相对误差随N 的变化规律，表明随着空间网格规模的增大，相对误差逐渐减小，这从另外一个角度说明：随着空间步长的减小，计算的精度在不断的提高。此外，在处理二维抛物型偏微分方程时，交替方向隐格（ADI算法）的应用非常广泛。为此，本人将ADI算法与上面所用的FST方法进行了对比：仍然使用原来的方程组，ADI算法的程序见附录（程序1.2），计算时两个程序都取空间步长，时间步长（即网格比），精确解取为时的函数值，对比的结果如下：表二、ADI算法与FST方法的对比N*2012345ADI误差5.33E-062.66E-062.48E-061.33E-061.49E-06ADI时间(s)0.651943829498125.6153281.737FST误差1.51E-057.47E-067.10E-073.73E-064.26E-07FST时间(s)0.34631621284832.7205176.93984 由表二可以看出在相同的计算规模下，FST算法的计算效率明显的高于ADI算法，同时当计算规模较大时，FST算法的计算误差也呈小于ADI的趋势。非齐次D氏边界条件如上所述，解决此类问题，可以采用ADI算法（交替方向隐格式），但边界条件的处理会比较麻烦且不能达到很高的计算精度，为此尝试采用课堂介绍的一种新算法，使得边界条件的处理变得相对简单一些，同时继续采用FST方法进行对称三对角矩阵求逆的快速处理。具体方法如下：对于带源项的控制方程，其C-N格式如下：由于其在时间上的精度为，下式与上式截断误差一致，可以认为是等效的：算法如下： 其中，第二步边界条件为，在本节中我们将之称为算法1。如此迭代，即可以求得各个时刻未知函数u的近似值，算法1对任意形式的源项函数f都适用。但这里我们考虑一种比较特殊的非齐次方程——本人在实际应用中经常遇到的情况，源项函数。也即其中，c为常数。虽然这类方程也可以使用上述方法，但我们也可以将源项处理到系数矩阵中，且通过实测，发现这种处理是可以提高计算精度的，具体推导如下：C-N差分格式类似第一题的结果：写成矩阵形式，即： （*）等价于又 ，因此，（*）式等价于也即上式和C-N格式有相同的时间精度，即。利用这个式子，就可以得到处理此类方程的完整算法： 其中，第二步仍然有边界条件为，这种算法称为算法2。如此迭代，既可以求得各个时刻位置函数u的近似值。计算时，采用的边界条件和初始条件如下：此方程的精确解为：空间步长取为，由于迭代的第一步使用的是显式欧拉方法，由格式的稳定性要求，可知网格比，程序中取，因此时间步长就是，计算时间点为，分别采用算法1和算法2对此方程进行了求解。具体程序见附录（程序2.1为算法1对应的程序，程序2.2为算法2对应程序），计算结果如下：表三、非齐次D氏边界条件下的计算结果（）N*201234568算法1误差0.0013920.0006670.0014060.0003230.0008270.0006850.000585时间s0.30818826755140.00886790185315