差分格式的截断误差【紧致差分格式的构造和验证】

题 目：
紧致差分格式的构造和验证 摘 要 目前，紧致差分格式已逐渐成为差分方程的数值方法的主要方向。具有良好特性的高精度的紧差分格式相继构造出来并能够应用到一些特殊的问题的数值求解，显现出了良好的效果。本课题针对紧致差分格式这一研究方向，希望能够通过MATLAB等软件的辅助以及前人对紧致差分格式的研究帮助对紧致差分格式进行构造一种差分格式，并且通过解微分方程的数值解实验对紧致差分格式进行验证其稳定性、收敛性以及误差等特性，最终能够比较直观了解这类紧致格式差分方法的精度等。

关键词:有限差分；
差分格式；
构造 ABSTRACT At present, compact difference schemes have gradually become a main research direction of the numerical method of differential equations, and the compact difference schemes with high precision and good characteristics have been constructed one after another and applied to the numerical solution of some specific problems, and good results have been achieved. This topic for compact difference scheme, the research direction of hope can through MATLAB software such as aided and previous study of compact difference scheme to help to construct a compact difference scheme difference scheme, and by solving the differential equation numerical solution of experiments to verify its compact difference scheme features such as stability, convergence and error, finally can more intuitive understanding of the compact format the precision of the finite difference method, etc. Key words：Finite difference; Difference scheme; Structure 目录 摘 要 2 ABSTRACT 3 1 引言 5 1.1 有限差分方法简介 5 1.2 紧致差分法研究概况 5 1.2.1 抛物线方程 5 1.2.2椭圆型方程 6 1.2.3双曲线方程 6 1.3 本文研究内容 6 2 常见差分格式 7 2.1 显式差分格式 7 2.1.1 古典显式格式的推导 7 2.2 隐式差分格式 8 2.2.1 古典隐式格式的推导 8 2.3 Crank-Nicolson隐式格式 10 2.4 交替方向隐式格式 11 2.4.1 Peaceman-Rachford格式 12 2.4.2 Douglas-Rachford格式 12 2.4.3 Mitchell-Fairweather格式 12 2.4.4 交替方向隐式格式算法步骤 12 3 紧致差分格式分析 13 3.1 抛物线方程 13 3.1.1 抛物线方程的一种高精度紧致差分方法 13 3.2 椭圆型方程 13 3.2.1一维椭圆型方程的解法 13 3.2.2 二维椭圆型方程的解法 14 3.3双曲型方程 15 3.3.1双曲线方程一种解法 15 3.3.2双曲线方程的常见数值解法 16 4实例分析与结果分析 17 4.1 数值算例 17 4.1.1 已知有精确解的热传导问题 17 4.1.2 未知精确解的热传导问题 18 4.2 结果分析 18 4.3 r变化对稳定性的探究 19 4.3.1 P-R格式格式的稳定性 19 4.4本文研究的热传导方程 20 5 总结 25 参考文献 26 1 引言 1.1 有限差分方法简介 重要的数值离散方法其中有有限差分方法（FDM），在研究、计算中有着广泛运用。有限差分法具有计算少、格式小、程序少等很多长处，尤其适合于偏微分方程的近似求解。有限差分法基本问题有：构造逼近微分方程定解问题的差分格式；
研究差分解的存在唯一性、收敛性及稳定性；
差分方程的解法等。创造出差分的方法有多种形式，目前主要采用的是泰勒方法。总的来说，一阶向前差分、一阶向后差分、二阶中心差分等都是比较常见的差分公式。这些年，有限差分方法有了很大的发展。在研究传统差分格式的基础上，学者们开始关注高精度的有限差分格式及其应用。高精度格式的采用放松了对网格步长的要求，采用计算区域内较少的格点得到较高精度的数值解，改进了求解效率，具有十分重要的理论和现实意义微分方程的数值作为求解方程的主要的近似解法，可以给出解在一些离散点上的近似值。现在求解微分方程的数值解法主要有两类，一类是有限差分方法，一类是有限元法。两者的基本思想基本上都是把连续问题（即微分方程的初边值问题）离散化，化为有限形式的线性代数方程组，求出原问题的离散解，再应用数值逼近的方法求出原问题的逼近解（连续解）。我们求得一般步骤是：开始，对求解区域进行网格剖分，用有限个网格节点代替连续区域；
然后，将微分方程离散化，从而把微分方程的定解问题化为线性代数方程组的求解问题。紧随之后微分方程数值解法研究的益深入，其在数值分析中也发挥着越来越重要的作用。因此，进行数值解法方面的研究，有着十分重要的理论和现实意义。

1.2 紧致差分法研究概况 伴随计算物理的发展，为了能更加方便研究与剖析差分问题，人们对于数值结果的精度要求越来越高，而普通的有限差分格式在很大步长下的精度偏低，为了提高精度不得不缩小步长，增加更多的网格点数量，这肯定增加了计算机的计算量，使得计算速度变慢，计算时间拉长。为了克服这一计算难题，人们开始寻找计算少精度高的差分方法。紧致差分方法正是在这种情形下应运而生。在1991年，lelei总结好了对称型紧致差分格式。相比于传统的差分格式，在相同的计算网格中，紧致差分格式有着更高的精度和分辨率。1993年傅德薰在紧致差分格式中引入迎风机制，1997年提出了五格点五阶精度的迎风紧致格式。紧致差分格式中迎风机制的引入能够有效抑制非物理振荡，更适合于多尺度复杂流场的计算。以上的对称型紧致差分格式和迎风型紧致差分格式是建立在均匀网格的基础上的。

1.2.1 抛物线方程 偏微分方程可用于描述出现在工程问题中的许多数学模型。抛物线问题作为偏微分方程中的一种频繁处理，在诸如扩散，渗流和热传导等问题中具有大量应用。因此，在很多情况下，偏微分方程找不到精确解，只能通过求解其数值解来研究问题。因此，找到高精度，小计算，稳定性好的数值计算方法具有重要意义。对于抛物问题，有限差分法经常用于解决数值问题。有限差分方法的基本思想是将连续解区域分离为有限点网格，然后使用网格上定义的离散变量差。该方程近似于连续解区域上的连续变量微分方程。最后，原始微分方程可以转换为代数方程组。求解方程可以得到离散点处原始微分方程的数值解。在众多差分方法中，紧凑差分格式使用较少的网格基点，但计算格式精度很高。因此，与传统的髙精度格式相比，紧凑的差分格式具有许多优点，例如计算量。小，对单元敏感，易于处理边界条件等。.因此，紧凑差分格式的研究是高精度格式研究的主要方向之一。目前，关于抛物问题的高精度紧致差分方法的研究就显而易见，由于稳定条件的过度限制，显式紧致差分方法难以被广泛使用。虽然隐式紧致差分方法具有更好的稳定性，但在求解方程时计算量大，计算效率低。

1.2.2椭圆型方程 针对一般椭圆型方程的数值求解，目前已有大量的相关文献，例如六阶精度差分格式用到了五个点，而十阶精度差分格式用到了七个点等等，并且推导出了高阶导数的紧致差分格式，同时分析了这些格式能够正确数值模拟的波数范围；
马延文和傅德薰两位学者构造了迎风紧致差分格式这样做的目的是为了克服数值解在激波附近的非物理高频振荡的难题，为了进一步提高精度，但不增加网格节点，他们推导出了超紧致差分格式沈孟育和蒋莉推导出了仅含有三个点的三阶精度紧致差分格式，但不足在于格式的系数推导求解过于复杂；
沈孟育等学者还将以上紧致差分格式做了进一步的改善，继而得到了广义紧致差分格式，上述所提到的不论是对称型紧致差分格式、或者迎风紧致差分格式还是超紧致差分格式都是它的特例；
田振夫和刘明会分别提出了一种求解两点边值问题的四阶紧致差分格式；
金涛等人提出了一种求解两点边值问题的的高阶隐式格式；
对于二维椭圆型方程，田振夫利用待定系数法和截断误差余项修正方法构造了数值求解二维泊松方程的四阶和六阶紧致差分格式；
陈国谦等人提出的求解对流扩散方程的四阶指数型格式，指数型摄动格式及迎风变换方法。

1.2.3双曲线方程 在偏微分方程中线性双曲型方程里，对于物理和生物领域的一些非线性现象，均可用此类方程来描述。例如，对流、扩散和反应、扩散之间的相互作用等。由于很难找到问题本身的精确解决方案，它是具有深远的意义，并通过数值方法来解决这样的方程实际应用价值。近年来，许多学者在国内外已经提出了求解一般线性双曲型方程的数值解的方法。莫汉蒂提出了无条件稳定的一维，二维和三维有限差分方法用于求解线性方程。刘等人提出了在狄利克雷和诺依曼边界条件下，求解一维电报方程的无条件稳定的两层紧致差分格式。空间上采用四次样条方法，时间上应用广义梯形公式。丁等人为了求解一维和二维电报方程，提出了一种条件稳定的高精度紧致差分格式。

1.3 本文研究内容 根据前人构造的格式，现在本文针对热传导方程建立了传统的差分格式p-r格式，做出数值实验；
基于理查德松外推法构造一种数值求解二维热传导方程的高阶紧致差分方法．在此基础上建立了一种紧差分格式。该方法首先利用时间二阶、空间四阶精度的紧致交替方向隐式差分格式在不同网格上对原方程进行求解，然后利用外推一次，最终得到了二维热传导方程时间四阶、空间六阶精度的数值解，数值实验验证了该方法的高阶精度及有效性。

2 常见差分格式 2.1 显式差分格式 现在，对于二维抛物线方程构建有限差分的显式格式 （2-1） 二维热传导方程 （2-2） 在热传导，磁扩散等许多领域有重要的应用。实际导热问题必然涉及边值条件，在有限差分法中它们也必须差分化。因此，我们需要研究的不仅是差分方程本身，而且是包括全部内部区域和所有边界上的差分方程所组成的代数方程组。又称为差分格式。

常系数热传导方程的古典显式格式首先考虑热传导方程的边值问题 离散化 [15] 2.1.1 古典显式格式的推导 最简单的显式差分为经典的显式格式，现在由热传导方程导出了 故 式中右边如果仅保留二阶导数项，且以替代，替代，则得差分格式 或者 （2-3） 这是一个显式格式（四点格式），如图（2-1）所示。

图2-1：古典显式格式 格式（2-3）应用在一维热传导相关方面中，则有古典显式格式：
（2-3） 通过一个方程组求解得到每一层各个节点上的值，显一隐格式区域分解方法就是以显式计算出比邻子域交界内的近似值的一种方法。显一隐格式区域分解方法综合了二者的优点，借助前一层数值解的信息，用显格式给出在这—层的子问题的我们所不知道德内边界条件，把一个整体区域上的问题化为很多个子区域上的子问题，在每个子域上采用古典隐式求解，从而实现了并行。这还能从下面的计算过程看出。

2.2 隐式差分格式 不同于显式方法，隐式差分格式中包括了（n+1）时间层上二个或二个以上结点处的未知值（例如 ），使用隐式差分格式和使用显式差分格式求解完全不同。相对而言，使用隐式差分格式求解，每时间层包含有较多的计算工作量。从后面对差分格式的稳定性分析可知，隐式格式的优点在于，其稳定性规则对步长比的限制大为放宽，我们正这么期望。对于二维抛物线方程，从公式开始，构建有限差分的显式格式。差分格式是由于热传导方程（2-1）构成。我们需要求解的不仅是差分方程本身，而且还包括所有的内部子区域，各边界的差分方程的代数方程。也被称为差分格式。

2.2.1 古典隐式格式的推导 现在对热传导方程推导其最简单的隐式差分逼近——古典隐式格式。由 故 式中右边如果仅保留二阶导数项，且以替代，替代，则得差分格式 或者 （2-4） 将格式（2-3）应用于一维热传导方程，古典显式格式为：
（3-2） 显式与隐式的比较如下：
（1）同一阶数下，隐式的局部截断误差的系数的绝对值比显式的要小；

（2）显式的计算工作量比隐式的小；

（3）隐式的稳定范围比显式的大。

局部截断误差的阶最高是3，式（2-5） 是允许函数任意变化情况下截断误差最小的二阶方法。要再提高阶就必须增加计算函数值的次数。上述式（2-5）又称为古典隐式格式。

故格式用下图3-1表示，我们知道截断误差阶为,与上述古典差分格式一样。这是一个四点的差分格式，正如下图3-1所示。

图3-1：古典隐式格式 为了求得第（n+1）时间层上的的值，必须通过解线性代数方程组。

这是一个隐式差分格式，必须联合其初边值条件求解。格式（2-5）通常称为古典隐式格式。我们也可以通过直接用差分算子代替 的方即：
代入微分方程，得到格式（2-4）。

古典隐式格式的方程组和矩阵形式如下：
当知道第n层上的时，要确定第n+1层上各点值必须通过求解一个线性代数方程组。以m=1为例：
其系数矩阵如下：
2.3 Crank-Nicolson隐式格式 Crank-Nicolson隐式差分格式是解热传导方程（2-1）的基本的差分格式，近年来，关于抛物方程的区域分解算法作为并行计算的一种有效工具，引起了很多学者的注意。这类算法的主要困难在于怎样定义内边界点的值和在子域上选取合理的计算解去近似。为了推导它，有 由 得 （2-6） 两边仅保留前二项，用代替，替代，则得差分格式 这是一个隐差分格式，称为Crank-nicolson差分格式，截断误差阶为，也可写为 （2-7） 由于格式（2-7）中包括六个结点，故也可称为六点格式（如图3-2所示）。

图3-2：
Crank-Nicolson隐式格式 也可将 代入微分方程，得到Crank-Nicolson格式。解微分方程，根据Crank-Nicolson格式得到的方程组：
 其矩阵表达式为：
2.4 交替方向隐式格式 交替方向隐式格式是x方向的显式格式和y方向的隐式格式。在n + 1/2层上以经典显式格式计算超出值，并且在n + 1层上以经典隐式格式校正预测值。获得的稳定性和收敛性是相似的。然而，仅讨论了z在常系数热传导方程方向上具有相同时间步长的情况，因此我们试图在此进行一些扩展。我们首先考虑内边界点使用一个方向作为显式格式而另一个方向作为隐式格式的情况，然后检查在z和y方向中使用显式格式的情况，我们讨论z方向的时间步骤的情况。

2.4.1 Peaceman-Rachford格式 考虑二维抛物方程（2-1）的差分，该显式格式的上述稳定性分析认为，要求比一维情况下更加苛刻。分析表明，尺寸越大，越小所需的时间步骤和更大的计算工作量。ADI格式从第n层到第n + 1层中，P-R格式在两个步骤中进行，每个步骤只需要求解方程的线性系统，而其具有三个对角系数。

P-R格式为 （2-8） （2-9） ，是中心差算子，， P-R格式对任意无条件稳定。但当我们预测P-R格式对三个变量空间问题，无条件稳定性不再成立。

2.4.2 Douglas-Rachford格式 道格拉斯于1956年提出另一个隐式差分格式，即Douglas-Rachford格式。第一次提出的D-R格式能被推广到三维情形的交替方向隐式格式，二维D-R格式是 （2-10） （2-11） 2.4.3 Mitchell-Fairweather格式 Mitchell和Fairweather在1964年推导了一个高精度ADI差分格式，称为M-F格式。M-F格式截断误差达较之P-R格式和D-R格式更好，且无条件收敛。

2.4.4 交替方向隐式格式算法步骤 ADI的格式的算法步骤类似于那些经典明确的格式，所不同的是ADI格式是一个隐式格式，需要通追赶方法来解决。追赶法是求解一个三对角矩阵的线性方程组的方法，而不是适用于其它类型的矩阵。

第一步，根据给出的初边值条件，得出t=0时; 第二步，根据各个交替方向隐式格式中的第一个公式，由第求得出如上式， 采用追赶法，确定追赶因子a=1+r/2，b= r/2，c=1+r/2，d=, =b-a*,=（d-*a）/w然后根据因子确定。=,= d -*; 第三步，再由追赶法由第求得出。

该算法，计算的近似解的过程。下文中将通过一个具体的数值例子说明了计算的方法，体现了这种格式的实用性和优越性。

3 紧致差分格式分析 3.1 抛物线方程 3.1.1 抛物线方程的一种高精度紧致差分方法 考虑一维抛物型方程的初边值问题 （3-1） 给定初始条件为 （3-2） 给定边界条件为 （3-3） 其中 u( x，t) 是待求未知函数，v 为扩散项系数，a( t) 、b( t) 、φ( x) 均为已知函数． 构造了差分格式：
（3-4） 上式即为所构造的差分格式，通过对空间方向的离散过程不难发现，其空间具有四阶精度． 另外，由文献中的研究结果我们能，时间方向也具有四阶精度，因此，本文中所构造的差分格( 3-4) 式的截断误差为O（t^4+h^4） 3.2 椭圆型方程 3.2.1一维椭圆型方程的解法 针对一维椭圆型方程，首先建立一个混合型髙精度紧致BCD差分格式，即格式整体精度为四阶；
然后分别分析此格式的截断误差。

考虑如下两点边值问题：
其中：均为己知函数，为常数．计算求解区间进行N等分，计算节点为，则网格步长，定义如下差分算子 （） 第一步：对u和u_{x}赋初值0，u的边界已知 第二步：通过u算出，即 观察得出，格式所形成的代数方程组系数矩阵为三对角线型的，本文采用追赶法进行求解． 第三步：由第二步算出的的值和（已知）算出新的u的内部值，即 由于边界点己知，因此该格式形成的代数方程组系数也是三对角线型的，故也可采用追赶法进行计算. 第四步：转入第二步的计算，循环进行下去，直到本次迭代算出的u值与上一次迭代算出的u值之间的误差按某种范数小于事先给定的收敛准则，则计算终止，输出u值，本章的数值实验部分取收敛准则为，为范数，定义为，k为迭代次数。

3.2.2 二维椭圆型方程的解法 考虑如下两点边值问题：
边界条件：
其中，为R^2上的矩形区域，为的边界，，均为已知函数，且假设其具有充分的光滑性． （） 第一步：对和赋初值0，的边界已知 第二步：通过算出和，即 观察得出，格式所形成的代数方程组系数矩阵为三对角线型的，本文利用追赶法进行求解. 第三步：由第二步计算出来的的值和f（已知）算出新的u的内部值，即 观察得出，该格式为９点模板，对该格式我们采用髙斯－赛德尔（Gauss－seider）迭代法进行求解． 第四步：转入第二步的计算，循环进行下去，直到本次迭代算出的值与上一次迭代算出的 值之间的误差按某种范数小于事先给定的收敛准则，则计算终止，输出值．本章的数值实验部分取收敛准则为，为范数，定义为，为迭代次数。

3.3双曲型方程 3.3.1双曲线方程一种解法 考虑如下一维线性双曲型方程的初边值问题：
其中是常数，均已知，u（x，t）未知。

首先，通过对时间二阶导数采用中心差分并保留截断误差主项，对空间导数项利用四阶紧致差分公式，之后采用四阶Ｐaｄé紧致差分格式进行逼近。整理差分格式为 其中 上式即为求解一维线性双曲型方程的四阶紧致差分格式，格式整体具有四阶精度。

3.3.2双曲线方程的常见数值解法 考虑如下的四阶双曲方程初边值问题：
（3-1） 其中为已知函数，且它们都满足一定的光滑性． 我们引入一个中间函数，使问题(3-1) 变成如下形式: （3-2） 可得问题(3-2) 的差分格式为 （3-3） 由前面的推导过程可知，上述格式（3-3）的局部截断误差为 4实例分析与结果分析 4.1 数值算例 4.1.1 已知有精确解的热传导问题 例1是二维Dirichlet边值问题的热传导方程 方程的精确解为 以P-R（ADI）格式，以相同的步长(，)求解同一处()的与解析解间的差值，列在误差表中如表所示。

表4-1 P-R（ADI）格式解例1的误差（，，） 误差 y=0 y=0.1 y=0.2 y=0.3 y=0.4 y=0.5 x=0 0 0 0 0 0 0 x=0.1 0 8.86E-07 1.19E-06 1.47E-06 6.44E-07 1.18E-06 x=0.2 0 1.19E-06 2.05E-06 2.19E-06 9.66E-07 9.93E-07 x=0.3 0 1.47E-06 2.19E-06 1.48E-06 9.21E-07 3.37E-09 x=0.4 0 6.44E-07 9.66E-07 9.21E-07 -8.6E-07 -1.4E-06 x=0.5 0 1.18E-06 9.93E-07 3.37E-09 -1.4E-06 -1.4E-06 x=0.6 0 6.44E-07 9.66E-07 9.21E-07 -8.6E-07 -1.4E-06 x=0.7 0 1.47E-06 2.19E-06 1.48E-06 9.21E-07 3.37E-09 x=0.8 0 1.19E-06 2.05E-06 2.19E-06 9.66E-07 9.93E-07 x=0.9 0 8.86E-07 1.19E-06 1.47E-06 6.44E-07 1.18E-06 x=1 0 0 0 0 0 0 续表4-1 P-R（ADI）格式解例1的误差（，，） 误差 y=0.6 y=0.7 y=0.8 y=0.9 y=1 x=0 0 0 0 0 0 x=0.1 6.44E-07 1.47E-06 1.19E-06 8.86E-07 0 x=0.2 9.66E-07 2.19E-06 2.05E-06 1.19E-06 0 x=0.3 9.21E-07 1.48E-06 2.19E-06 1.47E-06 0 x=0.4 -8.6E-07 9.21E-07 9.66E-07 6.44E-07 0 x=0.5 -1.4E-06 3.37E-09 9.93E-07 1.18E-06 0 x=0.6 -8.6E-07 9.21E-07 9.66E-07 6.44E-07 0 x=0.7 9.21E-07 1.48E-06 2.19E-06 1.47E-06 0 x=0.8 9.66E-07 2.19E-06 2.05E-06 1.19E-06 0 x=0.9 6.44E-07 1.47E-06 1.19E-06 8.86E-07 0 x=1 0 0 0 0 0 4.1.2 未知精确解的热传导问题 例2是二维Dirichlet边值问题的热传导方程。

采用P-R格式，由0变化到1时，不同的t, 关于，的图像（，，）对比图如图4-1。

图4-1：P-R格式解例2，变化 关于x，y的图像（，） 4.2 结果分析 由例1，可看出P-R格式精度较高，且P-R格式接近解析解。由例2，由图4-1可看出随时间的推移，温度渐渐靠近边界的温度。可见，数值方法也能一定程度上反应解的情况。

4.3 r变化对稳定性的探究 4.3.1 P-R格式格式的稳定性 采用P-R格式解例2，由0变化到1时，不同的时， 关于，的图像（，，）对比图如图4-2。采用P-R格式解例1，由0变化到1时，不同的时， 数值解与解析解间误差（，，变化）对比如表4-2。

图4-2：P-R格式解例2，变化 关于，的图像（，） 表4-2 P-R格式解例1的误差（，） 误差 r=1/8 r=2/8 r=3/8 r=4/8 x=0，y=0 0 0 0 0 x=0.1，y=0.1 8.86E-07 1.19E-06 -0.000464 -0.000489 x=0.2，y=0.2 2.05E-06 1.19E-06 -0.000883 -0.000928 x=0.3，y=0.3 1.48E-06 9.21E-07 -0.001214 -0.001275 x=0.4，y=0.4 -8.6E-07 9.66E-07 -0.001424 -0.001497 x=0.5，y=0.5 -1.4E-06 1.47E-06 -0.001497 -0.001574 x=0.6，y=0.6 -8.6E-07 1.18E-06 -0.001424 -0.001497 x=0.7，y=0.7 1.48E-06 1.47E-06 -0.001214 -0.001275 x=0.8，y=0.8 2.05E-06 2.05E-06 -0.000883 -0.000928 x=0.9，y=0.9 8.86E-07 1.47E-06 -0.000464 -0.000489 x=1.0，y=1.0 0 0 0 0 续表4-2 P-R格式解例1的误差（，） 误差 r=5/8 r=6/8 r=7/8 r=1 x=0，y=0 0 0 0 0 x=0.1，y=0.1 -0.000464 -0.000396 -0.000288 -0.000152 x=0.2，y=0.2 -0.000883 -0.000753 -0.000548 -0.000288 x=0.3，y=0.3 -0.001214 -0.001034 -0.000753 -0.000396 x=0.4，y=0.4 -0.001424 -0.001214 -0.000883 -0.000464 x=0.5，y=0.5 -0.001497 -0.001275 -0.000928 -0.000489 x=0.6，y=0.6 -0.001424 -0.001214 -0.000694 -0.000464 x=0.7，y=0.7 -0.001214 -0.001034 -0.000753 -0.000396 x=0.8，y=0.8 -0.000883 -0.000753 -0.000548 -0.000288 x=0.9，y=0.9 -0.000464 -0.000396 -0.000288 -0.000152 x=1.0，y=1.0 0 0 0 0 由图4-2，表4-2可验证，对P-R格式而言，r取0～1中的任何数，该格式都是稳 定的。

4.4本文研究的热传导方程 我们考察如下二维热传导方程的初边值问题：
式子中为边界；
为待求函数；
f（x，y，t）为源项；
g和为已知函数，且具有充分的光滑型；
a为热扩散系数或者导温系数，用表示时间步长，空间取等距网格，步长用ｈ表示．由文献［６］可得方程的高精度紧致ADI格式为 过渡变量为，边界条件可以精确表示为 其中和表示二阶中心差分算子，该格式的截断误差为。

下面采用Von Neumann分析法对格式进行稳定性分析。设源项没有误差，特征项表示为：
其中，I=为第n个时间层上的波幅，则放大因子可以表示成：
其中 很容易证明 因此，即格式是无条件稳定的。

采用Richardson外推法解下列问题：
选取正整数M和计算所需要到达的时刻Ｔ，并且令h=1/M,N=T/，而四阶紧致格式 可以写为 截断误差为：
假设为上述定解问题的解，为上述差分格式的解，记，则得到误差方程：
假设r（x，y，t）和s（x，y，t）分别满足二维热传导方程的初边值问题 和 其中 则 分别为紧致差分格式 和 的解，且 记 由上述定义，可得 其中 即为 同理，有 可得 ， 即分别在区域和上得到解和，其中。

综上，用外推法，可以得到方程（1）在上截断误差的近似值，即为 。

利用本文提出的方法求解如下有精确解的初边值问题，初始条件和边界条件均由精确解给出. 其精确解为：
下表4-3列出了当，时，在不同网格步长下的最大误差Error，收敛阶Rate和CPU时间以及与高精度紧致ADI方法计算结果的比较．其中可以看出，高精度紧致ADI格式在空间上的精度达到了四阶，而本文 Richardson外推法对于提高计算精度，效果是非常显著的，在空间上达到了六阶．在相同的网格等分数下，本文格式较文献格式[16]需要更多的CPU 时间，是因为本文格式需要在2个不同的网格层上进行计算，并且要用到Richardson外推，而文献[16]格式只需要在单层网格上进行计算．但如果我们限定计算需要达到的精度，比如1.44e-8，采用文献[16]的格式，需要取h=1/32，CPU 时间为0.34s，而采用本文方法，只需要取h=1/8，误差即可达到3.96e-9，而CPU 时间仅为0.04s． 表4-4给出了该问题当时刻，本文方法的最大误差Error，CPU 时间和收敛阶Rate 及与文献[17]中结果的比较．文献[17]采用紧致的ADI方法进行求解．通过对比可以看出，本文方法的精度明显高于文献[17]格式的精度．如果限定计算需要达到的精度，采用本文方法可以大大节约CPU时间． 表4-3 τ=h2,t=0.125时刻，最大误差、收敛阶和CPU时间 h 文献 本文方法 CPU Error rate CPU Error rate 1/8 0.01 3.52e-6 0.04 3.96e-9 1/16 0.06 2.21e-7 4.00 0.15 6.16e-11 6.01 1/32 0.34 1.44e-8 3.93 0.72 9.30e-13 6.05 表4-4 τ=h2,t=1时刻，最大误差、收敛阶和CPU时间 h 文献 本文方法 CPU Error rate CPU Error rate 1/10 0.14 1.71e-5 0.34 1.09e-9 1/20 0.58 1.22e-6 4.00 1.28 1.60e-11 6.01 1/40 2.46 7.90e-8 3.93 5.72 2.06e-13 6.05 总结 本文主要讨论有限差分法与紧致差分法，系统阐述了对于它们的数值解法，介绍了有限差分法和紧致差分方法的基本思路以及解题过程。着重讨论了各种方程的紧致差分解法，利用了理论计算进行分析了这些方法的收敛性和稳定性。以热传导方程作为一个例子，编程实现了使用P-R格式和紧致格式解二维抛物方程。

通过数值算例以及程序，验证了这两种格式的稳定性条件，当步长比，对任意步长比r，P-R格式都是稳定的。而本文紧致差分方法的精度明显比文献精度格式要更高．如果需要限制计算的精确度，采用本文方法可以大大节省CPU时间． 进一步确定了自己的研究方向：
实现差分格式步长自适应调整，可以提高计算结果的准确性，在不知道准确解的情况下，得到较准确的结果。

探究新的稳定性好，计算简单的紧致差分格式，争取构造更好的紧致方法。并进一步学习有限元法，有限体积法，有限单元法，无网格方法等算法。

参考文献 [1]陆静颖,葛永斌.数值求解一维波动方程的四阶紧致差分方法[J/OL].宁夏大学学报(自然科学版):1-9[2019-05-19] [2]韩俊茹,葛永斌.求解一维线性双曲型方程的高精度紧致差分格式[J].重庆师范大学学报(自然科学版),2018,35(06):64-69. [3]罗传胜,李春光,董建强,景何仿.求解对流扩散方程的一种高精度紧致差分格式[J].西南大学学报(自然科学版),2018,40(09):91-95. [4]李青,杨青.一维四阶双曲方程的紧致差分格式[J].山东师范大学学报(自然科学版),2018,33(02):154-158. [5]张含笑. 求解抛物型问题的混合型高精度紧致差分格式[D].宁夏大学,2018. [6]祁应楠.两点边值问题基于非均匀网格的高阶紧致差分格式[J].宁夏师范学院学报,2018,39(04):9-15+41. [7]梁昌弘. 椭圆型方程的混合型高精度紧致差分格式[D].宁夏大学,2017. [8]祁应楠.求解一维抛物型方程的高精度有限差分方法[J].西北师范大学学报(自然科学版),2016,52(06):29-34. [9]梁昌弘,马廷福,葛永斌.两点边值问题的混合型高精度紧致差分格式[J].宁夏大学学报(自然科学版),2017,38(01):1-4. [10]姜蕴芝. 求解波动方程的高精度紧致显式差分格式[D].宁夏大学,2016. [11]杨晓佳,魏剑英.求解抛物型方程的一种高精度紧致差分格式[J].湖北大学学报(自然科学版),2016,38(02):160-167+171. [12]王慧蓉.求解对流扩散方程的紧致二级四阶Runge-Kutta差分格式[J].云南民族大学学报(自然科学版),2015,24(05):382-385. [13]贾伟. 几种高精度紧致差分格式的构造与应用[D].西北师范大学,2015. [14]王小光. 紧致差分格式的优化和初步应用[D].河南师范大学,2012. [15]王芳芳. 紧致差分格式的理论及其分析[D].东北大学,2010. [16]Liao W,Zhu J,KHALIQAQM.An Efficient High-Order Algorithm for Solving Systems of reaction-Diffusion Equation[J].Numer Methods Partial Differential E q,2002,18(3):340-354 [17]DAI W,NASSAR R.Compact ADI Method for Solving Parabolic Differential Equations[J].Numer Methods Partial Differential E q,2002,18(2):129-146 [18]田振夫．非齐次热传导方程的高精度隐式格式[J].宁夏大学学报:自然科学版，199617(3)34-38．

相关热词搜索：构造 验证 差分 紧致差分格式的构造和验证 差分格式的截断误差 紧差分格式