下载: -
开放科学(资源服务)标志码(OSID):
-
抛物型方程是偏微分方程中的一种重要方程，它广泛应用于力学、天文学、物理学、生态学及工程技术等各个领域中，对于上述各领域中的许多问题，为了定量或定性地对它们展开系统的研究，需要建立各种数学模型，而建立起来的这些数学模型最后都可归结为求解相应的抛物型方程，这些抛物型方程是没有精确解的，我们只能通过计算机求它们的数值解，常见的数值解法有有限元法、有限体积法、边界元法和有限差分法等. 这些方法中，有限差分法仍然是求解抛物型方程的重要方法. 经典的差分格式，有古典显式格式、古典隐式格式、Crank-Nicolson格式等，这些差分格式形式简单且稳定性条件好，但它们的截断误差(精度)低，古典显式格式与古典隐式格式的精度仅为O(Δt+Δx2)，Crank-Nicolson格式的精度为O(Δt2+Δx2)[1-2]，与精确解的误差都较大. 因此，高精度且稳定性好的抛物型方程的差分格式的构造成了许多学者研究的问题.
生态学中提出的各种数学模型、神经轴突中电脉冲的传导及燃烧理论等许多物理现象，以及热的传导、流体在多孔介质中的运动规律等可归结为如下经典的抛物型方程-热传导方程，本文讨论二维的情形：
对方程(1)的求解，许多学者提出了一些改进的差分格式，大部分格式精度达到了O(Δt2+Δx4)[3-12]，如文献[12]构造了一个求解二维抛物型方程的三层显式格式，格式的精度为O(Δt2+Δx4)，稳定性范围是
$r<\frac{1}{2}$ . 在文献[11]中也构造了一个三层显式差分格式，格式的精度也是O(Δt2+Δx4)，稳定性范围是$\frac{1}{12} \leqslant r \leqslant \frac{1}{6}$ . 而在本文中则构造了一个精度为O(Δt3+Δx4)的高精度差分格式，格式的稳定性条件为$r<\frac{1}{6}$ ，精度比上述文献所构造的格式精度更高.
全文HTML
-
设Δt为时间步长，Δx与Δy为空间x，y方向的步长，为简便起见，设Δx=Δy，xj=jΔx，yk=kΔx，tn=nΔt，记ujkn=u(xj，yk，tn)，将ujkn+1与ujkn-1在节点(jΔx，kΔy，nΔt)处作Taylor展开得
化简得
为方便起见，记
它是一阶偏导数
$\left(\frac{\partial u}{\partial t}\right)_{j k}^{n}$ 的一个差分近似式.
-
利用上面的差分近似式和已知的差分近似式构造如下的差分格式来逼近方程(1)
其中：□ujkn=uj+1，k+1n+uj-1，k+1n+uj+1，k-1n+uj-1，k-1n-4ujkn，◇ujkn=uj+1，kn+uj-1，kn+uj，k+1n+uj，k-1n-4ujkn，
$\varDelta_{t} u_{j k}^{n}=\frac{u_{j k}^{n+1}-u_{j k}^{n}}{\varDelta t}, \eta_{t} u_{j k}^{n}=\frac{\frac{3}{2} u_{j k}^{n+1}-2 u_{j k}^{n}+\frac{1}{2} u_{j k}^{n-1}}{\varDelta t}$ ，θi(i=1，…，5)为待定参数.将(2)式中各节点上的u在节点(jΔx，kΔy，nΔt)处作Taylor展开，整理可得
令
$r=\frac{\varDelta t}{\varDelta x^{2}}$ ，为了使格式的截断误差达到O(Δt3+Δx4)，下面方程组成立解方程组(3)可得
将各参数的值代回(2)式中，即得到一个截断误差为O(Δt3+Δx4)的三层显式差分格式
-
下面利用Fourier分析法分析格式(4)的稳定性，首先写出与格式(4)等价的两层格式组
取
$u_{j k}^{n}=U^{n} \mathrm{e}^{\mathrm{i}(j \theta+k \varphi)}, v_{j k}^{n}=V^{n} \mathrm{e}^{\mathrm{i}(j \theta+k \varphi)}$ ，其中$\mathrm{i}=\sqrt{-1}$ ，通过计算可知□ujkn=-4s2ujkn，◇ujkn=-4s1ujkn，代入(5)式并整理得其中：
传播矩阵G(s1，s2)的特征方程为
引理1[13] 特征方程(6)的根满足|λ1，2|≤1的充要条件是
引理2[13] 差分格式(4)稳定，即矩阵族Gn(s1，s2)(0≤s1，s2≤2，n=1，2，…)一致有界的充要条件是
1) |λ1，2|≤1(λ1，2是方程(6)的两个根).
2) 使
$1-\frac{g_{11}^{2}}{4}=g_{11}^{2}+4 g_{12}=0$ 成立的(s1，s2)或不存在，或不属于区域[0, 2]×[0, 2].定理1 当
$r<\frac{1}{6}$ 时差分格式(4)稳定且收敛.证 由引理2的条件2)可知，当g12≠-1时，
$1-\frac{g_{11}^{2}}{4}=g_{11}^{2}+4 g_{12}=0$ 对所有的s1和s2都不成立，因此由(7)式知，格式(4)稳定的条件为-1+g12≤g11≤1-g12＜2.先讨论g11≤1-g12，可得
为确定起见，不妨假定
该式成立的条件是
当(9)式成立时，(8)式可化简为
该式成立的条件为
综合考虑(9)，(10)式，可得当
$r \leqslant \frac{1}{6}$ 时有g11≤1-g12.而当
$r \leqslant \frac{1}{6}$ 再由1-g12＜2可得由s1，s2的取值范围可知，(11)式成立的一个充分条件为
$r<\frac{1}{6}$ .最后由-1+g12≤g11可得
由s1，s2∈[0, 2]和
$r<\frac{1}{6}$ ，(12)式成立的一个充分条件为也即
而当
$r<\frac{1}{6}$ 时，(13)式恒成立.综上所述，由Lax的稳定性与收敛性定理，定理1得证.
-
在本节中，通过两个数值算例，利用MATLAB软件进行数值模拟，将本文提出的格式(4)与文献[11]和文献[12]格式进行比较，进一步证明本文提出的格式是一种高精度的差分格式.
例1 考虑如下带有初边值问题的二维热传导方程
在本例中取Δx=Δy=0.05，
$r=\frac{1}{12}, \frac{1}{10}, \frac{1}{7}, \frac{1}{5}$ ，由Δt=rΔx2可知Δt的取值分别为$\frac{1}{4800}, \frac{1}{4000}$ ，$\frac{1}{2800}, \frac{1}{2000}$ ，为方便计算，第一层的值ujk1利用方程(14)的精确解u(x，y，t)=e-2t sin (x+y)计算得到. 将本文格式、文献[11]格式及古典显格式进行比较，分别计算层数n=800时，部分节点处3种格式的数值解与精确解之间的绝对误差及3种格式所花费的时间(表 1).从表 1中我们可以看到，r在稳定性范围内时，本文格式比文献[11]格式和古典显式格式具有更高的精度，而当
$r=\frac{1}{5}$ 时，本文格式与文献[11]格式均不稳定，故数值解溢出. 在机器耗时方面，古典显式格式是一个两层格式，计算效率最高，而本文格式与文献[11]格式均为三层格式，计算效率差别不大，但都比古典显式格式低.例2
为验证差分格式(4)的空间精度与时间精度，取
$r=\frac{1}{8}$ ，Δx分别等于$\frac{\pi}{10}, \frac{\pi}{20}, \frac{\pi}{40}, \frac{\pi}{80}$ ，由Δt=rΔx2=$\frac{\varDelta x^{2}}{8}$ ，可知Δt分别等于$\frac{\pi^{2}}{800}, \frac{\pi^{2}}{3200}, \frac{\pi^{2}}{12800}, \frac{\pi^{2}}{51200}$ ，为方便计算，第一层的值ujk1利用方程(15)的精确解u(x，y，t)=e-t sin x sin y来计算.利用本文格式、文献[12]格式计算到层数n=500时，在部分节点处将这两种格式的数值解与精确解进行比较，结果如表 2. 从表 2可以看出，当$r=\frac{1}{8}$ ，本文格式相比文献[12]格式精度更高，随着Δx的逐渐减小，本文格式的数值解与精确解越来越接近，也验证了本文格式的收敛性.从以上两个数值算例可以看出，本文格式的收敛性及稳定性与理论分析一致，说明本文格式是一种有效的差分格式.