下载: -
化学、生物和物理领域的很多现象都可以用扩散反应方程来描述.由于物理问题本身精确解不易求出，所以寻求此类方程精确、高效、实用的数值方法有着重要的理论价值和实际意义.
求解扩散反应方程的数值方法包括有限差分法、有限体积法、有限元法等[1-6].目前已有的一维扩散反应方程的差分格式中时间与空间方向均达到四阶精度且绝对稳定的格式比较少见.因此，本文构造一种隐式高精度紧致差分格式，空间二阶导数项采用四阶紧致差分公式进行离散，时间导数项用四阶向后欧拉公式进行逼近，建立一种新的无条件稳定的隐式四阶紧致差分格式，然后通过数值算例验证本文方法的精确性和可靠性.
全文HTML
-
考虑如下一维扩散反应方程的初边值问题：
初始条件为：
边界条件为：
其中：u(x，t)为待求未知量，a＞0为扩散项系数，f(u)为非线性反应项，且φ(x)，g0(t)，gl(t)均为已知函数.不失一般性，设u(x，t)具有充分的光滑性.
设时间步长为τ，空间步长为
$ h = \frac{l}{N} $ ，网格节点为(xi，tn)，其中：xi=ih，tn=nτ，i=0，1，…，N，n＞0.为了便于推导，定义差分算子：
利用Kreiss四阶紧致差分公式[7]：
对方程(1)空间内部节点进行离散：
整理可得：
将(4)式代入(7)式整理可得：
考虑n+1时刻的值：
对式(9)时间导数uti+1n+1，utin+1，uti-1n+1用四阶向后欧拉公式[8]：
将式(10)代入式(9)可得：
令
$ r = \frac{\tau }{{{h^2}}} $ ，化简整理后，略去高阶项可得：其中n≥3，式(12)即为求解方程(1)的隐式四阶紧致格式.
从式(12)中可发现，所求未知量涉及5层.起动步除了初始时间层(第0个时间步)已知外，还须求得第1，2，3时间步的值，然后方可利用(12)式计算第4时间步及以后的值.因此还须构造前3个时间步的计算格式，为此考虑方程(10)在
$ n + \frac{1}{2} $ 时刻的值：对时间项利用中心差分格式离散，可得：
令
$ r = \frac{\tau }{{{h^2}}} $ ，(14)式可化简整理为：当n=0，1时，利用(15)式即可分别求出第1，2时间步的值.
接下来计算第3时间步的值，利用如下三阶向后欧拉公式[8]：
将式(16)代入式(9)中可得：
(17)式可化简整理为：
令n=2，即可通过(18)式得第3个时间步的计算格式.
由上述构造过程可知，本文所提方法在前2个时间步的计算时间仅具有二阶精度，第3个时间步时间具有三阶精度，第4个时间步以后，时间均具有四阶精度.尽管前3个时间步不能达到四阶精度，但是由于推进步数少，累积误差相对小，而在第4个时间步以后的计算中误差逐步减小，因此并不影响方法的整体精度，下面的数值实验将验证这一点.另外，文献[8]中提到，对于向后欧拉格式而言，当精度不高于六阶时，格式是稳定的；当精度高于六阶时，格式是不稳定的.而本文是四阶向后欧拉格式，故由文献[8]中的结论可知，本文格式是无条件稳定的.
-
为了验证本文格式的精确性和可靠性，考虑以下3个具有精确解的数值算例.所有计算采用Fortran77语言进行编程，且在PC机上采用的是双精度.分别给出不同空间步长、不同时间步长、不同网格比下的最大绝对误差Error及收敛阶Rate，其定义如下：
其中Error1和Error2为空间网格步长分别为h1和h2时的最大绝对误差.
例1精确解为
$ u\left( {x, t} \right) = {\left( {x - \frac{4}{3}} \right)^6}{{\rm{e}}^t} $ ，当t=1时的最大绝对误差及收敛阶见表 1.例2
例2精确解为u(x，t)=e-tcos(x)，在t=2，τ=h时的最大绝对误差及收敛阶见表 2.
例3
例3的精确解为u(x，t)=ex+t，在t=1，τ=h时的最大绝对误差及收敛阶见表 3.
例1是一个非齐次扩散问题，表 1给出了在t=1时刻的最大绝对误差和收敛阶，并与文献[1]、文献[9]和文献[10]的方法进行了比较.取相同参数下，文献[9]精度只有二阶.文献[10]条件稳定且精度只有二阶.文献[1]尽管达到了四阶精度，但由于格式时间仅具有二阶精度，因此必须取较小的时间步长τ=h2.本文格式时间具有四阶精度，因此可以取较大时间步长τ=h，同样可以达到四阶精度，且计算误差略优于文献[1].例2和例3分别为线性和非线性的扩散反应问题，对C-N格式收敛准则取tol=10-8，本文格式取tol=10-12.利用本文格式进行了计算并与Crank-Nicolson(C-N)格式所得结果进行了比较，取相同网格时，C-N格式只有二阶精度，而本文可以达到四阶精度，且计算误差比C-N格式小好几个数量级，通过对本文格式CPU时间的计算，可以发现文献[1]与C-N格式都是两层格式，而本文是5层格式，尽管比其它格式计算时间长，但是在计算允许误差给定的情况下，本文格式可以取较粗的网格，相比较之下，本文格式时间效率仍然很高.另外，由于本文格式计算均取τ=h，当
$ h = \frac{1}{{128}} $ 时，网格比$ r = \frac{\tau }{{{h^2}}} $ 最高可达128，故本文格式是无条件稳定的.
-
本文提出了数值求解一维扩散反应方程初边值问题的隐式高精度紧致差分格式，对空间二阶导数采用四阶紧致差分格式进行离散，时间导数项采用四阶向后欧拉公式进行计算，格式的截断误差为Ο(τ4+h4)，即格式的时间和空间均具有四阶精度.目前对扩散反应方程所提高精度问题的研究多为两层或三层隐格式，而本文采用的隐格式时间层是5层，随着时间步的增加，计算误差在逐渐减小，且时间和空间可以同时达到四阶精度.通过与文献[1, 9-10]中方法和经典的C-N格式进行比较可知，本文计算结果更精确.本文方法可以推广到二维甚至三维空间去，比如文献[11]和文献[12]所考虑的模型方程，对此我们将另作讨论.