求解一维对流扩散方程的高精度紧致差分格式

求解一维扩散反应方程的隐式高精度紧致差分格式1概述一维扩散反应方程是描述许多物理过程的数学方程之一，如化学反应、热传导等。

在求解这样的方程时，我们需要寻找适合的数值解法。

本文将介绍一种隐式高精度紧致差分格式，用于求解一维扩散反应方程。

2一维扩散反应方程一维扩散反应方程可表示为：$$\frac{\partial u}{\partial t}=D\frac{\partial^2u}{\partial x^2}+\rho u(1-u)$$其中，$u(x,t)$表示物理量的变量，$D$为扩散系数，$\rho$为反应速率常数。

初始条件为$u(x,0)=u_0(x)$，边界条件为$u(0,t)=u(L,t)=0$，其中$L$为区间长度。

3差分方法为了求解上述方程的数值解，我们需要使用差分方法。

差分方法可以将连续的偏微分方程转化为离散的方程，从而得到数值解。

这里我们采用一阶差分法和二阶差分法分别对时间和空间进行离散化。

时间离散化：$$\frac{\partial u(x,t)}{\partialt}\approx\frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}$$空间离散化：$$\frac{\partial^2u(x,t)}{\partialx^2}\approx\frac{u(x+\Delta x,t)-2u(x,t)+u(x-\Deltax,t)}{\Delta x^2}$$将上述两个式子带入到原方程中，得到离散化形式：$$\frac{u_i^{n+1}-u_i^n}{\Delta t}=D\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2}+\rho u_i^n(1-u_i^n)$$其中，$n$表示时间步长，$i$表示空间位置。

4隐式高精度紧致差分格式在上述差分方法中，我们采用了一阶差分法和二阶差分法，这种方法的精度有限。

为了提高求解的精度，可以采用更高阶的差分方法。

第48卷第8期西南师范大学学报(自然科学版)2023年8月V o l.48N o.8 J o u r n a l o f S o u t h w e s t C h i n aN o r m a lU n i v e r s i t y(N a t u r a l S c i e n c eE d i t i o n)A u g.2023D O I:10.13718/j.c n k i.x s x b.2023.08.006求解对流扩散反应方程的高阶指数型组合紧致差分格式①王明镜,田芳,郭亚妮宁夏大学数学统计学院,银川750021摘要:针对一维对流扩散反应方程,基于对流扩散方程的四阶指数型紧致差分格式,采用四阶混合差分逼近算子和P a dé公式,间接地构造了两种求解一维对流扩散反应方程的四阶指数型组合紧致差分格式;针对二维对流扩散反应方程,采用降维法,结合高阶混合差分逼近算子和P a dé公式构造了求解二维对流扩散反应方程的四阶指数型组合紧致差分格式.本文所提差分格式较经典四阶格式和文献中的组合型格式具有更低的耗散性,因此对于对流占优等边界层问题的求解计算精度更高.最后给出数值算例验证了本文格式的精度.关键词:对流扩散反应方程;高阶紧致差分格式;对流占优;边界层中图分类号:O241.82文献标志码:A文章编号:10005471(2023)08004114AH i g h-O r d e rE x p o n e n t i a l C o m b i n a t i o nC o m p a c t F i n i t eD i f f e r e n c e S c h e m e f o rC o n v e c t i o n-D i f f u s i o m-R e a c t i o nE q u a t i o nWA N G M i n g j i n g, T I A NF a n g,G U O Y a n iS c h o o l o fM a t h e m a t i c sa n dS t a t i s t i c s,N i n g x i aU n i v e r s i t y,Y i n c h u a n750021,C h i n aA b s t r a c t:F o r t h e o n e-d i m e n s i o n a l c o n v e c t i o nd i f f u s i o n r e a c t i o ne q u a t i o n,b a s e do n t h e f o u r t h-o r d e r e x p o-n e n t i a l c o m p a c t d i f f e r e n c e s c h e m e o f c o n v e c t i o n d i f f u s i o n e q u a t i o n,t w o f o u r t h-o r d e r e x p o n e n t i a l c o m b i n e d c o m p a c t d i f f e r e n c e l a t t i c e s f o r s o l v i n g o n e-d i m e n s i o n a l c o n v e c t i o nd i f f u s i o n r e a c t i o n e q u a t i o n a r e c o n s t r u c-t e d i n d i r e c t l y b y u s i n g f o u r t h-o r d e rm i x e dd i f f e r e n c ea p p r o x i m a t i o no p e r a t o ra n dP a déf o r m u l a;F o r t h e t w o-d i m e n s i o n a l c o n v e c t i o nd i f f u s i o nr e a c t i o ne q u a t i o n,a f o u r t h-o r d e r e x p o n e n t i a l c o m b i n e dc o m p a c t d i f f e r-e n c e s c h e m e f o r s o l v i n g t h e t w o-d i m e n s i o n a l c o n v e c t i o nd i f f u s i o n r e a c t i o n e q u a t i o n i s c o n s t r u c t e db y u s i n g t h e d i-m e n s i o n r e d u c t i o nm e t h o d,c o m b i n i n g t h e h i g h-o r d e rm i x e dd i f f e r e n c e a p p r o x i m a t i o no p e r a t o r a n dP a déf o r m u l a. T h e d i f f e r e n c e s c h e m e p r o p o s e d i n t h i s p a p e r i s l e s s d i s s i p a t i v e t h a n t h e c l a s s i c a l f o u r t h-o r d e r s c h e m e a n d t h e c o m-b i n e d s c h e m e i n t h e l i t e r a t u r e,s o i t h a s h i g h e r a c c u r a c y f o r s o l v i n g t h e c o n v e c t i o n d o m i n a t e d b o u n d a r y l a y e r p r o b-l e m s.F i n a l l y,a n u m e r i c a l e x a m p l e i s g i v e n t o v e r i f y t h e a c c u r a c y o f t h e p r o p o s e d s c h e m e.K e y w o r d s:c o n v e c t i o n d i f f u s i o n r e a c t i o n e q u a t i o n;h i g h o r d e r c o m p a c t d i f f e r e n c e s c h e m e;c o n v e c t i o n d o m-i n a t e s;b o u n d a r y l a y e r①收稿日期:20221020基金项目:宁夏自然科学基金项目(2020A A C03059);国家自然科学基金项目(11902170;11772165;12161067);宁夏自治区青年拔尖人才培养工程项目.作者简介:王明镜,硕士研究生,主要从事偏微分方程数值解法的研究.通信作者:田芳,副教授.对流扩散反应方程已被广泛用于研究模拟热流㊁污染物的扩散㊁化学反应等物理现象[1].求解对流扩散反应方程的数值方法有很多[2-19].针对此类方程,研究者们发展了许多高精度有限差分格式[12-19].常见的差分格式依据差分格式系数的函数类型不同可分为多项式型有限差分格式[14-15]和指数型有限差分格式[2,16],依据有限差分格式结构的区别则可以分为方程型有限差分格式[11-17]和导数型有限差分格式[18].本文的主要工作是基于文献[2]提出的对流扩散方程的四阶指数型有限差分格式,间接得到求解一维㊁二维对流扩散反应方程的四阶指数型组合紧致差分格式,所提出的差分格式不同于导数型差分格式,其优点在于采用高阶差分算子分步迭代计算待求量的导数值,避免求解大型代数方程组,同时还兼具了指数型格式高精度的优点.1一维对流扩散反应方程的四阶指数型组合紧致差分格式1.1对流扩散方程的四阶指数型紧致差分格式对于如下对流扩散模型方程-εu x x+μu x=S(x),xɪ[X1,X2](1)文献[2]中给出了形如-αδ2x u i+μδx u i=S i+c1S x i+c2S x x i(2)的四阶指数型紧致差分格式,系数为α=μh2c o t hμh2εæèçöø÷,μʂ0ε,μ=0ìîíïïïc1=ε-αμ,μʂ00,μ=0ìîíïïïc2=ε(ε-α)μ2+h26,μʂ0h212,μ=0ìîíïïïï(3)此处δx u i和δ2x u i分别定义为δx u i=u i+1-u i-12h(4)δ2x u i=u i+1-2u i+u i-1h2(5)该格式对于求解对流占优等大梯度变化的对流扩散问题具有高精度㊁高分辨率的优点.本文将基于该格式构造对流扩散反应方程的四阶指数型组合紧致差分格式.1.2对流扩散反应方程的四阶指数型组合紧致差分格式本文考虑如下对流扩散反应方程-εu x x+μu x+b(x)u=f(x),xɪ[X1,X2](6)其中:ε为扩散项系数且ε>0,μ为对流项系数;b(x)和f(x)是关于x的足够光滑的函数;u(x)是待求未知量;当b(x)=0时,模型方程(6)即为文献[2]中的对流扩散方程.首先将求解区间[X1,X2]等分为N个子区间:x i=X1+i h,i=0,1,2, ,N,h=X2-X1N.在点x i处由泰勒级数展开得到u x i=δx u i-h26∂3u∂x3æèçöø÷i+O(h4)(7)u x x i=δ2x u i-h212∂4u∂x4æèçöø÷i+O(h4)(8)由(7)式得u x x i=δx u x i-h26∂4u∂x4æèçöø÷i+O(h4)(9)再由(8)式和(9)式可得待求量的二阶导数u x x i的四阶组合差分逼近公式u x x i=2δ2x u i-δx u x i+O(h4)(10)下面对对流扩散反应模型方程(6)的四阶指数型组合紧致差分格式进行推导.将模型方程(6)改写为如24西南师范大学学报(自然科学版)h t t p://x b b j b.s w u.e d u.c n第48卷下形式-εu x x +μu x =S (x )S (x )=f (x )-b (x )u {(11)对(11)式中的第一个方程考虑其在点x i 处形如(2)式的四阶指数型紧致差分格式-αδ2x u i +μδx u i =S i +c 1S x i +c 2S x x i(12)其中α,c 1,c 2由(3)式给出.对(11)式中的第二个方程两端的x 求导得S x =f x -b x u -b u x ,S x x =f xx -2b x u x -b x x u -b u x x (13)将(13)式代入(12)式中,整理得-αδ2x u i +μδx u i +(b i +c 1b x i +c 2b x x i )u i =f i +c 1f x i +c 2f x x i +(-c 1b i -2c 2b x i )u x i -c 2b i u x x i (14)由(10)式得f x x i =2δ2x f i -δx f xi +O (h 4)并将其代入(14)式右端,略去高阶项整理得-αδ2x u i +μδx u i +(b i +c 1b x i +c 2b x x i )u i =(1+2c 2δ2x )f i +(c 1-c 2δx )f xi +(-c 1b i -2c 2b x i )u x i -c 2b i u x x i (15) 下面通过对(15)式右端项中待求量的二阶导数采用不同的处理方法,得到以下求解模型方程(6)的两种四阶指数型组合紧致差分格式.1.2.1 本文格式1将(10)式代入(15)式中,消去u x x 项整理得-(α-2c 2b i )δ2x u i +μδx u i +(b i +c 1b x i +c 2b x x i )u i =(1+2c 2δ2x )f i +(c 1-c 2δx )f xi +(-c 1b i -2c 2b x i +c 2b i δx )u x i 即得本文格式1:-P δ2x u i +Q δx u i +R u i =F i(16)其中P =α-2c 2b i ,Q =μ,R =b i +c 1b x i +c 2b x x iF i =(1+2c 2δ2x )f i +(c 1-c 2δx )f xi +(-c 1b i -2c 2b x i +c 2b i δx )u x i (17)1.2.2 本文格式2假设εʂ0时,将原模型方程(6)改写为u x x i =-1ε(f i -μu x i -b i u i )将其代入(15)式右端项,消去u x x 项整理得-αδ2x u i+μδx u i +b i +c 1b x i +c 2b x x i +c b 2i εæèçöø÷u i =1+c 2b i ε+2c 2δ2x æèçöø÷f i +(c 1-c 2δx )f x i +-c 1b i -2c 2b x i -c 2μb i εæèçöø÷u x i 即得本文格式2:-P δ2x u i +Q δx u i +R u i =F i(18)其中P =α,Q =μ,R =b i +c 1b x i +c 2b x x i +c 2b 2iεF i =1+c 2b i ε+2c 2δ2x æèçöø÷f i +(c 1-c 2δx )f x i +-c 1b i -2c 2b x i -c 2μb i εæèçöø÷u x i (19)1.3 本文格式右端项F i 的计算若使(16)式和(18)式逼近模型方程(6)的精度达到四阶,则需要对F i 中的u x i ,f x i 采用四阶差分逼近.此处采用如下四阶P a d é型紧致差分公式进行离散16(U x )i -1+23(U x )i +16(U x )i +1=δx U i ,i =1,2, ,N (20)其中U 代表u 或f .离散u x i 和f x i ,边界条件采用如下四阶逼近公式[20]34第8期 王明镜,等:求解对流扩散反应方程的高阶指数型组合紧致差分格式(U x )0+89(U x )1=1h -22190U 0+433108U 1-196U 2+4318U 3-2527U 4+27180U 5æèçöø÷(U x )N -89(U x )N -1=1h 2110U N -641108U N -1+12118U N -2-256U N -3+4127U N -4-43180U N -5æèçöø÷ìîíïïïï(21)1.4 F o u r i e r 误差分析令u ~=e i k x,i 为虚数单位,则δx u ~=i s i n k h h u ~,δ2xu ~=2h2(c o s k h -1)u ~且由(20)式得u ~x=i 3s i n k h h (2+c o s k h )u ~,f ~x =i 3s i n k h h (2+c o s k h )f~又令ω=k h ,P e =μεh ,C o u n t =b μh ,η=P e ㊃C o u n t =b εh 2,则微分方程(6)精确的特征函数为λE x a c t =εh2(ω2+η)+i (ω㊃P e )[](22)经过计算可得文献中的F O C 格式[14]㊁文献[18]格式㊁MHO C 格式[19]㊁本文格式1和本文格式2的特征函数均可统一地表示为如下形式λ=εh 2M 1M 2(23)其中相应的M 1和M 2的表达式详见表1.表1 差分格式的特征函数差分格式M 1M 2F O C 格式[14]-2-P e 26()(c o s ω-1)+η6㊃(c o s ω+5)[]+i ㊃s i n ωP e -η㊃P e12()5+c o s ω6-i P e ㊃s i n ω12MHO C 格式[19]4(c o s ω-1)+η+3(s i n ω)22+c o s ω[]+i ㊃P e ㊃3s i n ω2+c o s ω1文献[18]格式24(1-c o s ω)+η2-ηP e 2+12η[]+i ㊃(12P e -P e 3)s i n ωi ㊃P e 2+2(c o s ω+5)[]-i ㊃P e ㊃s i n ω本文格式1-2γ(c o s ω-1)+4σ2η(c o s ω-1)+η+3σ2η㊃(s i n ω)22+cos ω[]+i ㊃P e ㊃s i n ω+3σ1η㊃s i n ω2+cos ω()1+4σ2(c o s ω-1)+3σ2(s i n ω)22+c o s ω[]+i ㊃3σ1s i n ω2+c o s ω本文格式2-2γ(c o s ω-1)+η+σ2η2[]+i ㊃P e ㊃s i n ω+3σ1η㊃s i n ω2+c os ω+3σ2P e ㊃η㊃s i n ω2+c o s ω()1+ησ24σ2(c o s ω-1)+3σ2(s i n ω)22+c o s ω[]+i ㊃3σ1s i n ω2+c o s ω图1表示了在区间0ɤk h ɤπ上,当C o u n t 数确定为常数1.0时,对于不同的网格雷诺数P e ,数值波数和精确波数之间的关系.图中κ2h 2表示耗散误差.结果表明:在P e =1.0时,MHO C 格式耗散误差很大,而F O C 格式和文献[18]格式的耗散误差曲线完全重合,并且相比其他格式,本文格式1和本文格式2的耗散误差较小;当网格雷诺数P e 从10增大到1000时,本文格式1的耗散误差逐渐增大,表现出强耗散性,而本文格式2的耗散误差要明显小于其它格式的耗散误差.44西南师范大学学报(自然科学版) h t t p ://x b b jb .s w u .e d u .c n 第48卷图1 C o u n t =1.0的耗散误差比较图2表示了在区间0ɤk h ɤπ上,当网格雷诺数P e 确定为常数1.0时,对于不同的C o u n t 数,数值波数和精确波数之间的关系.结果表明:对于不同的C o u n t 数取值,本文格式1和本文格式2的耗散误差均明显小于F O C 格式的耗散误差,并且在高波段上本文格式2的耗散误差小于本文格式1的耗散误差.图2 P e =1.0时的耗散误差比较54第8期 王明镜,等:求解对流扩散反应方程的高阶指数型组合紧致差分格式图3表示了在区间0ɤk h ɤπ上,当P e 数和C o u n t 数相等时,对于不同的取值,数值波数和精确波数之间的关系.结果表明:对于不同的取值,在整个波段上,本文格式2的耗散误差很小,而其它4种格式的耗散误差随着P e 数和C o u n t 数的增大而增大,表现出强耗散性.特别地,文献F O C 格式和文献[18]格式的耗散误差曲线完全重合.图3 P e =C o u n t =c 时的耗散误差比较1.5 算法描述将(16)式和(18)式中的算子展开,差分方程统一写为A -1u i -1+A 0u i +A 1u i +1=B -1f i -1+B 0f i +B -1f i +1+C -1f x i -1+C 0f x i +C 1f x i +1+D -1u x i -1+D 0u x i +D 1u x i +1(24)其中对应(24)式的系数见(17)式和(19)式.求解算法流程图如图7所示.图4 算法流程图64西南师范大学学报(自然科学版) h t t p ://x b b jb .s w u .e d u .c n 第48卷2 二维对流扩散反应方程的四阶指数型组合紧致差分格式2.1 对流扩散反应方程的四阶指数型组合紧致差分格式考虑如下二维对流扩散反应问题-ε1u x x -ε2u y y +μ1u x +μ2u y +b (x ,y )u =f (x ,y ),(x ,y )ɪ[X 1,Y 1]ˑ[X 2,Y 2](25)其中:ε1,ε2为扩散项系数且ε1,ε2>0;μ1和μ2为对流项系数;b (x ,y )为反应项系数;f (x ,y )为源项;u (x ,y )是未知函数,函数u (x ,y ),b (x ,y )和f (x ,y )在计算区域上足够光滑.将求解区间[X 1,X 2]等分为N 个子区间:x i =X 1+i h ,y j =Y 1+j h ,i ,j =0,1,2, ,N ,h =X 2-X 1N.将(25)式改写为如下等价的方程组形式-ε1u x x +μ1u x +b (x ,y )u =f 1(x ,y )-ε2u y y +μ2u y +b (x ,y )u =f 2(x ,y ){(26)其中f 1(x ,y )=f (x ,y )+ε2u y y -μ2u y f 2(x ,y )=f (x ,y )+ε1u x x -μ1u x {(27)对(26)式和(27)式应用差分格式(14)得-P 1δ2x u i j +Q 1δx u i j +R 1u i j =F 1i j-P 2δ2y u i j+Q 2δy u i j +R 2u i j =F 2i j{(28)其中F 1i j =f 1i j +c 1f 1x i j +c 2f 1x x i j +(-c 1b i j -2c 2b x i j )u x i j -c 2b i j u x x i j(29)F 2i j =f 2i j +d 1f 2y i j +d 2f 2y y i j +(-d 1b i j -2d 2b y i j )u y i j -d 2b i j u y y i j(30)P 1=α1,Q 1=μ1,R 1=b i j +c 1b x i j +c 2b x x i j α1=μ1h 2c o t h μ1h 2ε1æèçöø÷,μ1ʂ0ε1,μ1=0ìîíïïïï,c 1=ε1-α1μ1,μ1ʂ00,μ1=0ìîíïïïï,c 2=ε1(ε1-α1)μ21+h 26,μ1ʂ0h 212,μ1=0ìîíïïïï(31)P 2=α2, Q 2=μ2, R 2=b i j +d 1b y i j +d 2b y y i j α2=μ2h 2c o t h μ2h 2ε2æèçöø÷,μ2ʂ0ε2,μ2=0ìîíïïïï,d 1=ε2-α2μ2,μ2ʂ00,μ2=0ìîíïïïï,d 2=ε2(ε2-α2)μ22+h 26,μ2ʂ0h 212,μ2=0ìîíïïïï(32)将(28)式中两项相加得-P 1δ2x u i j -P 2δ2y u i j +Q 1δx u i j +Q 2δy u i j +(R 1+R 2)u i j =F i j(33)其中F i j =F 1i j +F 2i j(34)由f 1(x ,y )=f (x ,y )+ε2u y y -μ2u y (35)直接计算得f 1x =f x +ε2u y y x -μ2u y x f 1x x =f x x +ε2u y y x x -μ2u y x x (36)将(35)式和(36)式代入(29)式中得F 1i j =(f +ε2u y y -μ2u y )i j +c 1(f x +ε2u y y x -μ2u y x )i j +c 2(f x x +ε2u y y x x -μ2u y x x )i j +(-c 1b i j -2c 2b x i j )u x i j -c 2b i j u x x i j =74第8期 王明镜,等:求解对流扩散反应方程的高阶指数型组合紧致差分格式(f+c1f x+c2f x x)i j+(ε2u y y i j-μ2u y i j)+c1(ε2u y y x i j-μ2u y x i j-b i j u x i j)+c2(ε2u y y x x i j-μ2u y x x i j-b i j u x x i j-2b x i j u x i j)(37)由f2(x,y)=f(x,y)+ε1u x x-μ1u x(38)直接计算得f2y=f y+ε1u x x y-μ1u x yf2y y=f y y+ε1u x x y y-μ1u x y y(39)将(38)式和(39)式代入(30)式中得F2i j=(f+ε1u x x-μ1u x)i j+d1(f y+ε1u x x y-μ1u x y)i j+d2(f y y+ε1u x x y y-μ1u x y y)i j+ (-d1b i j-2d2b y i j)u y i j-d2b i j u y y i j=(f+d1f y+d2f y y)i j+(ε1u x x i j-μ1u x i j)+d1(ε1u x x y i j-μ1u x y i j-b i j u y i j)+d2(ε1u x x y y i j-μ1u x y y i j-b i j u y y i j-2b y i j u y i j)(40)将(37)式和(40)式相加并整理得F i j=F1i j+F2i j=(f+c1f x+c2f x x+d1f y+d2f y y)i j+b u i j+c1(ε2u y y x i j-μ2u y x i j-b i j u x i j)+d1(ε1u x x y i j-μ1u x y i j-b i j u y i j)+c2(ε2u y y x x i j-μ2u y x x i j-b i j u x x i j-2b x i j u x i j)+d2(ε1u x x y y i j-μ1u x y y i j-b i j u y y i j-2b y i j u y i j)(41)对其中的c2(-μ2u y x x i j-b i j u x x i j-2b x i j u x i j)+d2(-μ1u x y y i j-b i j u y y i j-2b y i j u y i j)中各项采用二阶中心差分算子离散,对其中的c1(ε2u y y x i j-μ2u y x i j-b i j u x i j)+d1(ε1u x x y i j-μ1u x y i j-b i j u y i j)中各项采用以下四阶差分算子离散u x y i j=δx u y i j+δy u x i j-δxδy u i j+O(h4)(42)u x x y i j=δ2x u y i j+δ2xδy u i j-δxδy u x i j+o(h4)(43)u y y x i j=δ2y u x i j+δ2yδx u i j-δxδy u y i j+o(h4)(44)得F i j=(f+c1f x+c2f x x+d1f y+d2f y y)i j+b u i j+(c2ε2+d2ε1)δ2xδ2y u i j+(-c2μ2+d1ε1)δ2xδy u i j+(c1ε2-d2μ1)δ2yδx u i j-(-c1μ2-d1μ1)δxδy u i j-c2b i jδ2x u i j-d2b i jδ2y u i j-2c2b x i jδx u i j-2d2b y i jδy u i j+ (c1ε2δ2y-d1ε1δxδy+(-c1μ2-d1μ1)δy-c1b i j)u x i j+(d1ε1δ2x-c1ε2δxδy+(-c1μ2-d1μ1)δx-d1b i j)u y i j (45)将(45)式代入(33)式中并整理得(-P1+c2b i j)δ2x u i j+(-P2+d2b i j)δ2y u i j+(Q1+2c2b x i j)δx u i j+(Q2+2d2b y i j)δy u i j+(-c2ε2-d2ε1)δ2xδ2y u i j+(c2μ2-d1ε1)δ2xδy u i j+(d2μ1-c1ε2)δ2yδx u i j+(-c1μ2-d1μ1)δxδy u i j+(R1+R2-b i j)u i j= (f+c1f x+c2f x x+d1f y+d2f y y)i j+(-c1b i j+c1ε2δ2y-d1ε1δxδy+(-c1μ2-d1μ1)δy)u x i j+(-d1b i j+d1ε1δ2x-c1ε2δxδy+(-c1μ2-d1μ1)δx)u y i j (46)从而得到模型方程(25)的离散差分格式-Aδ2x u i j-Bδ2y u i j+Cδx u i j+Dδy u i j+Eδ2xδ2y u i j+Gδ2xδy u i j+Hδ2yδx u i j+Kδyδx u i j+R u i j=F i j(47)其中A=P1-c2b i j,B=P2-d2b i j,C=Q1+2c2b x i j,D=Q2+2d2b y i j,E=-c2ε2-d2ε1G=c2μ2-d1ε1,H=d2μ1-c1ε2,K=-d1μ1-c1μ2,R=R1+R2-b i j F i j=f i j+c1f x i j+d1f y i j+c2f x x i j+d2f y y i j+(-c1b i j-(d1μ1+c1μ2)δy+c1ε2δ2y-d1ε1δxδy)u x i j+(-d1b i j-(d1μ1+c1μ2)δx+d1ε1δ2x-c1ε2δxδy)u y i j (48)84西南师范大学学报(自然科学版)h t t p://x b b j b.s w u.e d u.c n第48卷(48)式中右端项中的f x i j ,f y i j ,f x x i j ,f y y i j 分别采用一阶导数和二阶导数的中心差分算子计算得到.而u x i j ,u y i j 采用如下的四阶Pa d é型紧致差分公式进行离散16(u x )i -1,j +23(u x )i j +16(u x )i +1,j =δx u i j16(u y )i ,j -1+23(u y )i j +16(u y )i ,j +1=δy u i j ìîíïïïïi =1,2, ,N (49)离散u x i j和u y i j ,边界条件采用文献[20]给出的如下四阶逼近公式,其中α代表x 或y .(u α)0+89(u α)1=1h -22190u 0+433108u 1-196u 2+4318u 3-2527u 4+27180u 5æèçöø÷(u α)N -89(u α)N -1=1h 2110u N -641108u N -1+12118u N -2-256u N -3+4127u N -4-43180u N -5æèçöø÷ìîíïïïï(50)2.2 算法描述1)给u 赋初值u 0i j ,i ,j =1,2, ,N ;2)采用中心差分算子计算f x i j ,f y i j ,f x x i j ,f y y i j ;3)采用四阶P a d é型紧致差分公式(49)和四阶一致边界条件(50)计算u x i j ,u y i j ;4)将以上各量值代入(48)式计算F i j ;5)将F i j 代入(47)式中计算u i j 的新值,记为u k +1i j ;6)如果e r r =u k +1i j -u k i j u k i j ȡ10-14,则将u k +1i j 赋值于u 0i j ,重复2)-5),直到e r r =u k +1i j -u k i j u ki j<10-14,输出u i j .3 数值算例本节将选取典型算例,采用F o r t r a n 95程序语言,在具有4G B 内存,I n t e l (R )C o r e (T M )i 5-7200U C P U 处理器的计算机上编制统一程序,采用本文格式和文献格式进行计算,并与精确解进行比较,验证格式的精度.收敛阶采用如下公式R a t e =l og E 1E 2æèçöø÷l o g N 2N 1æèçöø÷(51)进行计算,其中:E 1,E 2分别为网格数取N 1,N 2时采用差分格式计算得到的最大绝对误差.例1[15]-u x x +u x +u =2s i n x +c o s x ,0<x <π.该算例的边界条件为u (0)=u (π)=0.精确解为u (x )=si n x .它在计算区域上是一个光滑函数.表2给出了在不同网格数下采用以上5种格式计算的最大绝对误差和收敛阶.从表2中我们可以发现,5种差分格式在粗网格和细网格上均能达到理论上的计算精度,且本文格式1和格式2的误差略小于F O C ㊁文献[18]和MHO C3种四阶差分格式的计算误差.表2 算例1最大绝对误差和收敛阶比较网格数F O C 格式[14]最大绝对误差收敛阶文献[18]格式最大绝对误差收敛阶MHO C 格式[19]最大绝对误差收敛阶本文格式1最大绝对误差收敛阶本文格式2最大绝对误差收敛阶86.0149e -56.0149e -54.8660e -44.7342e -54.5602e -5163.8124e -63.983.8124e -63.981.8332e -54.733.0021e -63.982.9826e -63.93322.3734e -74.012.3734e -74.016.0042e -74.931.9034e -73.981.9005e -73.97641.4877e -84.001.4877e -84.001.9470e -84.951.1896e -84.001.1892e -84.001289.2959e -104.009.2960e -104.009.0068e -104.437.4354e -104.007.4345e -104.002565.8107e -114.005.8106e -114.004.7573e -114.244.6478e -114.004.6387e -114.0094第8期 王明镜,等:求解对流扩散反应方程的高阶指数型组合紧致差分格式例2 -εu x x -u x +u =0,0<x <1.边界条件为u (0)=0,u (1)=e m 1-e m 2(1+m 1)e m 1-(1+m 2)e m 2,精确解为u (x )=e m 1x -e m 2x(1+m 1)e m 1-(1+m 2)e m 2.其中m 1=-1+1+4ε2ε,m 2=-1-1+4ε2ε.如图5所示,当网格数N =16时,随着ε的减小,该算例的精确解在计算区域的左端点x =0处有一边界层.表3比较了当ε=10-1,10-2,10-3时本文格式和文献格式的计算结果.图5 例2的精确解表3 算例2最大绝对误差和收敛阶比较ε网格数F O C 格式[14]最大绝对误差收敛阶文献[18]格式最大绝对误差收敛阶MHO C 格式[19]最大绝对误差收敛阶本文格式1最大绝对误差收敛阶本文格式2最大绝对误差收敛阶0.183.7874e -41.7386e -32.3689e -37.6095e -54.4308e -5162.1887e -54.119.1882e -54.243.9102e -42.603.5697e -64.412.8621e -63.95321.4128e -63.955.7785e -63.992.9563e -53.731.9798e -74.171.9004e -73.91648.7848e -84.013.5728e -74.021.4512e -64.351.2019e -84.041.1913e -84.000.01162.8707e -29.4464e -1-3.6560e -45.8710e -4325.8496e -32.299.1507e -23.37--1.2391e -41.568.0193e -52.87646.0300e -43.284.2661e -34.42--1.5992e -52.959.6289e -63.061283.7333e -54.012.2505e -44.24--8.7058e -74.206.7542e -73.832562.2955e -64.021.3290e -54.08--4.5294e -84.264.3165e -83.970.001161.5129e -12.3159e -1-5.3801e -46.6504e -3321.2552e -10.272.4033e -1-0.05--2.5441e -41.082.8622e -31.22648.5562e -20.554.5989e -1-0.94--1.1368e -41.167.0721e -42.021284.0174e -21.099.4108e -1-1.03--4.8890e -51.221.1984e -42.562561.0205e -21.984.8920e -10.94--1.9168e -51.351.5517e -52.955121.2815e -32.991.0767e -25.51--3.4552e -62.472.0184e -62.9410249.0851e -53.825.8991e -44.19--2.3248e -73.891.6744e -73.5920485.4469e -64.063.3217e -54.15--1.1435e -84.351.0624e -83.9840963.3684e -74.022.0217e -64.04--6.7813e -104.086.6657e -103.9905西南师范大学学报(自然科学版) h t t p ://x b b jb .s w u .e d u .c n 第48卷从表3可知:1)在细网格下,所列格式都能到理论上的收敛阶,但本文两种差分格式的计算精度明显优于F O C 格式㊁MHO C 格式和文献[18]格式.2)随着ε的减小,本文两种差分格式优势更显著.比如当ε=10-3,N =512时,采用本文格式1计算得到的最大绝对误差为3.4552ˑ10-6,采用本文格式2计算得到的最大绝对误差为2.0184ˑ10-6;F O C 格式若要得到相同量级的计算精度(此时最大绝对误差为5.4469ˑ10-6),至少需要2048个网格点;文献[18]格式若要获得相同量级的计算精度(此时最大绝对误差为3.7697ˑ10-6),至少需要3500个网格点.3)当ε很小时,本文格式1在粗网格上的计算误差小于本文格式2,在细网格上的计算误差大于本文格式2.当ε=1和ε=10-1时,MHO C 格式可以达到理论上的四阶精度;但是当ε<0.065时,不论在粗网格还是在细网格上,MHO C 格式都是发散的.例3[16]-εu x x -u x +εu =s i n x ,0<x <π边界条件为u (0)=u (π)=0,精确解为u (x )=11+4ε2A e r 1x π+B e r 2xπ+2εs i n x +c o s x ().其中r 1=-π(1+1+4ε2)2ε,r 2=2πε(1+1+4ε2)2,A =(1+e r 2)(e r 1-e r 2),B =(1+e r1)(e r 2-e r 1).如图6所示,随着ε的减小,该算例的精确解在计算区域的左端点x =0处有一边界层.表4中比较了当ε=10-1,10-2,10-3时本文格式和文献格式的计算结果.从表4可知:图6 例3的精确解(ε=10-1,10-2,10-3,10-5,10-8,网格数N =16)1)在细网格下,所列格式都能达到理论上的收敛阶,但本文两种差分格式的计算精度明显优于F O C格式㊁MHO C 格式和文献[18]格式.2)随着ε的减小,本文两种差分格式的计算误差比所列其它文献格式高出2~6个数量级.当ε=10-3,N =32时,本文格式1的最大绝对误差为3.6641ˑ10-5,本文格式2的最大绝对误差为4.5278ˑ10-4;而F O C 格式(此时最大绝对误差为3.5389ˑ10-4)若要得到相同量级的计算精度,则至少需要4096个网格点;文献[18]格式若要得到相同量级的计算精度(此时最大绝对误差为5.2628ˑ10-4),则至少需要6000个网格点.3)当ε很小时,本文格式1在粗网格上的计算误差小于本文格式2,而在细网格上的计算误差大于本文格式2.当ε=1时,MHO C 格式可以达到理论上的四阶精度;但是当ε<0.23时,不论在粗网格还是在细网格上,MHO C 格式都是发散的.15第8期 王明镜,等:求解对流扩散反应方程的高阶指数型组合紧致差分格式Copyright ©博看网. All Rights Reserved.表4算例3最大绝对误差和收敛阶比较ε网格数F O C格式[14]最大绝对误差收敛阶文献[18]格式最大绝对误差收敛阶MHO C格式[19]最大绝对误差收敛阶本文格式1最大绝对误差收敛阶本文格式2最大绝对误差收敛阶0.189.5799e-21.7778e+0-1.8617e-31.5209e-3 161.2180e-22.981.8956e-13.23--3.3120e-42.491.9728e-42.95 328.7155e-43.802.2321e-23.09--2.2389e-53.891.6179e-53.61 645.2192e-54.064.6621e-32.26--1.1035e-64.341.0259e-63.980.01161.0685e+03.9618e+1-2.3825e-41.3615e-3 325.8212e-10.884.5983e+03.11--6.9312e-51.782.4469e-42.48 641.8500e-11.652.8066e+00.71--2.9283e-51.243.2805e-52.90 1282.9092e-22.673.2921e-13.09--7.3488e-61.994.3467e-62.92 2562.4340e-33.581.9755e-24.06--6.5607e-73.494.3379e-73.32 5121.4258e-44.091.5887e-33.64--3.0836e-84.412.7628e-83.970.001161.8075e+03.8920e+4-2.8944e-41.8337e-3 321.7625e+00.042.4552e+33.99--3.6641e-52.984.5278e-42.02 641.5638e+00.171.5510e+23.98--4.5424e-63.011.0677e-42.08 10245.7787e-29.1693e-1-1.2960e-78.4136e-8 20485.8109e-33.314.3509e-24.40--1.6025e-83.029.8384e-93.10 40963.5389e-44.042.3231e-34.23--8.5555e-104.236.7449e-103.87例4[21]-εu x x-εu y y+2u x+u y=f(x,y),0<x,y<1.边界条件为u(0,y)=e y+2-1ε(1+y)1+1ε,u(1,y)=e y-1+2-1ε(1+y)1+1ε,u(x,0)=e-x+2-1ε, u(x,1)=e1-x+2.精确解为u(x,y)=e y-x+1+y2æèçöø÷1ε(1+y).表5算例4最大绝对误差和收敛阶比较ε网格数文献[18]格式最大绝对误差收敛阶本文格式(47)最大绝对误差收敛阶18ˑ84.0956e-61.6858e-416ˑ162.5522e-74.002.5339e-52.53 32ˑ321.5953e-84.003.4814e-62.86 64ˑ649.9851e-104.004.5672e-72.930.116ˑ162.0305e-32.9483e-432ˑ321.1016e-44.203.7478e-52.98 64ˑ646.6592e-64.054.7185e-62.99 128ˑ1284.1305e-74.015.9340e-72.990.0164ˑ641.4384e+02.4627e-3128ˑ1286.5368e-37.781.5297e-44.01 256ˑ2563.2225e-44.349.7320e-63.97 512ˑ5121.9228e-54.076.3561e-73.940.001512ˑ512o v e r f l o w5.6155e-31024ˑ10241.8534e-2-3.4187e-44.04 2048ˑ20488.1406e-44.512.2305e-53.9425西南师范大学学报(自然科学版)h t t p://x b b j b.s w u.e d u.c n第48卷Copyright©博看网. All Rights Reserved.表5列出了当ε=1,10-1,10-2,10-3时本文格式和文献格式的最大绝对误差和收敛阶比较.从表5中数据不难发现,ε越小,本文格式的优势越明显.特别地,当ε=10-2,10-3时,本文格式较文献格式在计算精度上要高出1~2个数量级.例5 -εu x x -εu y y -u x -u y +εu =s i n x +s i n y ,0<x ,y <π.边界条件为u (0,y )=11+4ε2A 1+e r 1y π()+B 1+e r 2y π()+2ε(s i n y )+(1+c o s y )()u (π,y )=11+4ε2A e r 1+e r 1y π()+B e r2+e r 2yπ()+2ε(s i n y )+(-1+c o s y )()u (x ,0)=11+4ε2A e r 1x π+1()+B e r 2xπ+1()+2ε(s i n x )+(c o s x +1)()u (x ,π)=11+4ε2A e r 1x π+e r 1()+B e r 2xπ+e r2()+2ε(s i n x )+(-1+c o s x )()精确解为u (x ,y )=11+4ε2(A (e r 1x π+e r 1y π)+B (e r 2x π+e r 2yπ)+2ε(s i n x +s i n y )+(c o s x +c o s y )).其中r 1=-π(1+1+4ε2)2ε,r 2=2πε(1+1+4ε2)2,A =(1+e r 2)(e r 1-e r 2),B =(1+e r1)(e r 2-er 1).表6列出了当ε=1,10-1,10-2时本文格式和文献格式的最大绝对误差和收敛阶比较.从表6中数据不难看出,当ε=1时,本文格式和文献格式的计算精度相当.但是,当ε=10-1,10-2时,本文格式的计算精度远远高于文献格式.表6 例5最大绝对误差和收敛阶比较ε网格数文献[18]格式最大绝对误差收敛阶本文格式(47)最大绝对误差收敛阶18ˑ84.2758e -52.4427e -416ˑ162.7642e -63.952.2248e -53.4632ˑ321.7354e -73.991.5975e -63.8064ˑ641.0874e -84.001.0539e -73.920.116ˑ164.7910e -13.0075e -332ˑ321.7188e -24.801.7757e -44.0864ˑ649.4600e -44.189.1586e -64.28128ˑ1285.7461e -54.049.8648e -73.220.01128ˑ1282.0467e +06.4493e -5256ˑ2565.8375e -26.175.4368e -63.57512ˑ5122.8236e -33.333.1255e -74.121024ˑ10241.6520e -44.101.8342e -84.094 结论本文基于对流扩散方程的指数型差分格式[2],构造了求解一维和二维对流扩散反应方程的四阶指数型组合紧致差分格式,所选取的典型算例的计算结果充分表明对于对流(反应)占优问题,本文提出的差分格式的计算精度要明显高于文献中的多项式格式[14]和多项式型组合格式[18-19].参考文献:[1]MO R T O N K W.N u m e r i c a l S o l u t i o no fC o n v e c t i o n -D i f f u s i o nP r o b l e m s [M ].B o c aR a t o n :C R CP r e s s ,2019.[2] T I A NZF ,D A I SQ.H i g h -O r d e rC o m p a c tE x p o n e n t i a l F i n i t eD i f f e r e n c eM e t h o d s f o rC o n v e c t i o n -D i f f u s i o nT y peP r o b -35第8期 王明镜,等:求解对流扩散反应方程的高阶指数型组合紧致差分格式Copyright ©博看网. All Rights Reserved.45西南师范大学学报(自然科学版)h t t p://x b b j b.s w u.e d u.c n第48卷l e m s[J].J o u r n a l o fC o m p u t a t i o n a l P h y s i c s,2007,220(2):952-974.[3]陈文兴,戴书洋,田小娟,等.利用重心有理插值配点法求解一㊁二维对流扩散方程[J].西南师范大学学报(自然科学版),2020,45(8):35-43.[4] P E N GD Y.H i g h-O r d e rN u m e r i c a lM e t h o d f o rT w o-P o i n tB o u n d a r y V a l u eP r o b l e m s[J].J o u r n a l o fC o m p u t a t i o n a lP h y s i c s,1995,120(2):253-259.[5] A Y U S OB,MA R I N I LD.D i s c o n t i n u o u sG a l e r k i nM e t h o d s f o rA d v e c t i o n-D i f f u s i o n-R e a c t i o nP r o b l e m s[J].S I AMJ o u r-n a l o nN u m e r i c a lA n a l y s i s,2009,47(2):1391-1420.[6] P HO N G T HA N A P A N I C H S,D E C HA UM P HA IP.C o m b i n e dF i n i t eV o l u m ea n dF i n i t eE l e m e n t M e t h o df o rC o n v e c-t i o n-D i f f u s i o n-R e a c t i o nE q u a t i o n[J].J o u r n a l o fM e c h a n i c a l S c i e n c e a n dT e c h n o l o g y,2009,23(3):790-801. [7] S H I H Y,K E L L O G G RB,T S A IP.A T a i l o r e dF i n i t eP o i n tM e t h o df o rC o n v e c t i o n-D i f f u s i o n-R e a c t i o nP r o b l e m s[J].J o u r n a l o f S c i e n t i f i cC o m p u t i n g,2010,43(2):239-260.[8] A N G E L I N IO,B R E N N E R K,H I L HO R S TD.AF i n i t eV o l u m eM e t h o do nG e n e r a lM e s h e s f o r aD e g e n e r a t eP a r a b o l i cC o n v e c t i o n-R e a c t i o n-D i f f u s i o nE q u a t i o n[J].N u m e r i s c h eM a t h e m a t i k,2013,123(2):219-257.[9] B A U S E M,S C HW E G L E RK.H i g h e rO r d e r F i n i t eE l e m e n tA p p r o x i m a t i o n o f S y s t e m s o f C o n v e c t i o n-D i f f u s i o n-R e a c t i o nE q u a t i o n sw i t hS m a l lD i f f u s i o n[J].J o u r n a l o fC o m p u t a t i o n a l a n dA p p l i e d M a t h e m a t i c s,2013,246:52-64.[10]A D A K D,N A T A R A J A NE,K UMA RS.A N e w N o n c o n f o r m i n g F i n i t eE l e m e n tM e t h o d f o rC o n v e c t i o nD o m i n a t e dD i f-f u s i o n-R e a c t i o nE q u a t i o n s[J].I n t e r n a t i o n a lJ o u r n a lo f A d v a n c e si n E ng i n e e r i n g S c i e n c e sa n d A p p l i e d M a th e m a ti c s,2016,8(4):274-283.[11]朱海涛,欧阳洁.对流扩散反应方程的变分多尺度解法[J].工程数学学报,2009,26(6):997-1004.[12]兰斌,薛文强,葛永斌.对流扩散反应方程基于坐标变换的高阶紧致差分格式[J].青岛科技大学学报(自然科学版),2014,35(1):100-106.[13]J HA N,S I N G H B.E x p o n e n t i a lB a s i sa n dE x p o n e n t i a lE x p a n d i n g G r i d sT h i r d(F o u r t h)-O r d e rC o m p a c tS c h e m e s f o rN o n l i n e a rT h r e e-D i m e n s i o n a lC o n v e c t i o n-D i f f u s i o n-R e a c t i o n E q u a t i o n[J].A d v a n c e si n D i f f e r e n c e E q u a t i o n s,2019, 2019(1):1-27.[14]S P O T Z W F.H i g h-O r d e rC o m p a c t F i n i t eD i f f e r e n c e S c h e m e s f o rC o m p u t a t i o n a lM e c h a n i c s[D].A u s t i n:T h eU n i v e r s i t yo fT e x a s a tA u s t i n,1995.[15]S U N H W,Z HA N GJ.A H i g h-O r d e rF i n i t eD i f f e r e n c eD i s c r e t i z a t i o nS t r a t e g y B a s e do nE x t r a p o l a t i o nf o rC o n v e c t i o nD i f f u s i o nE q u a t i o n s[J].N u m e r i c a lM e t h o d s f o rP a r t i a lD i f f e r e n t i a l E q u a t i o n s,2004,20(1):18-32.[16]R A D HA K R I S HN AP I L L A IAC.F o u r t h-O r d e rE x p o n e n t i a l F i n i t eD i f f e r e n c eM e t h o d s f o rB o u n d a r y V a l u eP r o b l e m so fC o n v e c t i v eD i f f u s i o nT y p e[J].I n t e r n a t i o n a l J o u r n a l f o rN u m e r i c a lM e t h o d s i nF l u i d s,2001,37(1):87-106.[17]田芳,葛永斌.求解变系数对流扩散反应方程的指数型高精度紧致差分方法[J].工程数学学报,2017,34(3):283-296.[18]杨苗苗,葛永斌.求解对流扩散反应方程的一种高精度紧致差分方法[J].四川师范大学学报(自然科学版),2021,44(4):470-478.[19]梁昌弘,马廷福,葛永斌.两点边值问题的混合型高精度紧致差分格式[J].宁夏大学学报(自然科学版),2017,38(1):1-4.[20]王涛,刘铁钢.求解对流扩散方程的一致四阶紧致格式[J].计算数学,2016,38(4):391-404.[21]S A N Y A S I R A J U Y VSS,M I S H R A N.S p e c t r a l R e s o l u t i o n e dE x p o n e n t i a l C o m p a c tH i g h e rO r d e r S c h e m e(S R E C HO S)f o rC o n v e c t i o n-D i f f u s i o nE q u a t i o n s[J].C o m p u t e rM e t h o d s i nA p p l i e d M e c h a n i c s a n dE ng i n e e r i n g,2008,197(51/52):4737-4744.责任编辑张栒Copyright©博看网. All Rights Reserved.。

应用数学MATHEMATICA APPLICATA2019,32(3):635-642求解一维对流方程的高精度紧致差分格式侯波,葛永斌(宁夏大学数学统计学院,宁夏银川750021)摘要：本文提出数值求解一维对流方程的一种两层隐式紧致差分格式,采用泰勒级数展开法以及对截断误差余项中的三阶导数进行修正的方法对时间和空间导数进行离散.格式的截断误差为O(τ4+τ2h2+h4),即该格式在时间和空间上均可达到四阶精度.利用von Neumann方法分析得到该格式是无条件稳定的.通过数值实验验证了本文格式的精确性和稳定性.关键词：对流方程;高精度;紧致格式;无条件稳定;有限差分法中图分类号：O241.82AMS(2000)主题分类：65M06;65M12文献标识码：A文章编号：1001-9847(2019)03-0635-081.引言对流方程在生物数学、能源开发、空气动力学等许多领域都具有十分广泛的应用,因此求解该类方程具有非常重要的理论价值和实际意义.然而,由于实际问题通常十分复杂,往往难以求得精确解,因此研究其精确稳定的数值解法是十分必要的.针对对流方程国内外很多学者提出了很多的数值方法.如张天德和孙传灼[1]针对一维对流方程采用待定系数法,得到了两层四点格式和四阶六点格式,并且是无条件稳定的,该方法适用于在点数确定的前提下,得到精度高的差分格式;于志玲和朱少红[2]针对一维问题建立了中间层为两个节点的三层显格式,其截断误差为O(τ2+h2);曾文平[3]针对一维对流方程推导出了一种两层半显式格式，其截断误差为O(τ2+h2),该格式是无条件稳定的.姚朝辉等人[5]将二阶的迎风格式和中心差分格式进行加权得到了WSUC格式,该格式是无条件稳定的;但该格式时间方向和空间方向仅有二阶精度.汤寒松等人[6]通过立方插值拟质点方法(CIP方法),给出了一些保单调的CIP格式;Erdogan[9]针对一维的对流方程推导出了一种指数拟合的差分格式,其截断误差为O(τ2+h2);Bourchtein[10]构造了对流方程的三层五点中心型蛙跳格式,该格式的截断误差为O(τ4+h4);即该格式时间和空间均具有四阶精度,但是该格式是三层的,空间方向需要五个点,并且是条件稳定的;Kim[11]构造了多层无耗散的迎风蛙跳格式,即时间和空间分别具有二阶、四阶、六阶精度,但格式为三层甚至是四层的,并且六阶格式空间方向最多需要五个点,给靠近边界的内点的计算带来困难.综上所述,文献中已经有的数值计算方法大多为低阶精度的,而高精度方法涉及多个时间层,需要一个或多个时间启动步,或者空间方向的网格节点多于三个,这都给计算造成困难或不便.为此本文将构造一种紧致格式,这里紧致格式的定义为对时间导数项的离散采用不超过∗收稿日期：2018-08-10基金项目：国家自然科学基金(11772165,11361045),宁夏自然科学基金重点项目(2018AAC02003),宁夏自治区重点研发项目(2018BEE03007)作者简介：侯波,男,汉族,河南人,研究方向：偏微分方程数值解法.通讯作者：葛永斌.636应用数学2019三个时间层,而对空间导数项的离散采用不超过三个网格点,时间和空间即可以达到高阶精度（三阶及三阶以上）的格式.本文拟构造的格式时间方向仅用到两个时间层上的函数值,在每个时间层上仅涉及到三个空间网格点,格式时间和空间具有整体的四阶精度.该格式的优点是无须启动步的计算,并且在对靠近边界点的计算时,不会用到计算域以外的网格节点.此外该格式为无条件稳定的,可以采用比较大的时间步长进行计算.最后通过数值实验验证本文格式的精确性和稳定性.2.差分格式的建立考虑如下一维对流方程：∂u ∂t +a∂u∂x=f,b≤x≤c,t≥0,(2.1)给定初始条件为:u(x,0)=φ(x),b≤x≤c,(2.2)给定周期性边界条件为:u(b,t)=u(c,t),t≥0,(2.3)其中,u(x,t)为未知函数,f为非齐次项,a为对流项系数,φ(x)为已知函数.将求解区域[b,c]等距剖分为N个子区间：b=x0,x1,···,x N−1,x N=c,并且定义h=c−bN,时间也采用等距剖分,步长用τ表示.在本文中,我们利用u ni ,u n+1i,u n+12i分别表示u在(x i,t n),(x i,t n+1)和(x i,t n+12)点处的函数值.假设方程(2.1)在点(x i,t n+12)成立,简写表示为:(∂u ∂t )n+12i+a(∂u∂x)n+12i=f n+12i.(2.4)将u n+1i 和u ni在点(x i,t n+12)处做泰勒级数展开,可得:u n+1i=u n+12i+τ2(∂u∂t)n+12i+(τ2)22!(∂2u∂t2)n+12i+(τ2)33!(∂3u∂t3)n+12i+O(τ4),(2.5)u ni=u n+12i−τ2(∂u∂t)n+12i+(τ2)22!(∂2u∂t2)n+12i−(τ2)33!(∂3u∂t3)n+12i+O(τ4).(2.6)(2.5)-(2.6)可得:(∂u∂t)n+12i=δt u n+12i−τ224(∂3u∂t3)n+12i+O(τ4),(2.7)其中,δt u n+12i =u n+1i−u n iτ.同理可得:(∂u∂x)n+12i=δx u n+12i−h26(∂3u∂x3)n+12i+O(h4),(2.8)其中,δx u n+12i =un+12i+1−u n+12i−12h.将(2.7)和(2.8)代入(2.4)整理可得:δt u n+12i −τ224(∂3u∂t3)n+12i+aδx u n+12i−ah26(∂3u∂x3)n+12i=f n+12i+O(τ4+h4).(2.9)为了使该格式在时间方向和空间方向上均达到四阶精度,须对(2.9)式中的∂3u∂t3和∂3u∂x3项进行二阶的离散,同时为了保证本文格式的紧致性,即空间方向不超过三个网格点,我们对(2.1)式进行如下变形:∂u ∂t =−a∂u∂x+f,∂2u∂t2=a2∂2u∂x2−a∂f∂x+∂f∂t,第3期侯波等：求解一维对流方程的高精度紧致差分格式637∂3u ∂t 3=a 2∂3u ∂x 2∂t −a ∂2f ∂x∂t +∂2f ∂t 2,∂3u ∂x 3=−1a ∂3u ∂x 2∂t +1a ∂2f ∂x 2.(2.10)将上述∂3u ∂t 3和∂3u∂x 3的表达式(2.10)代入(2.9)并整理可得：δt u n +12i+aδx u n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t)n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+h 4).(2.11)如果对上式中的δx u n +12i 项采用时间方向算术平均,即δx u n +12i =δx u n +1i+u n i 2,则会导致格式时间退化为二阶精度,为此利用(2.5)+(2.6)可得：u n +12i =12(u n +1i +u n i )−τ28(∂2u ∂t2)n +12i +O (τ4).(2.12)从而可得:δx u n +12i =12δx (u n +1i +u n i )−τ28δx (∂2u ∂t2)n +12i +O (τ4).(2.13)将(2.13)代入(2.11)得:δt u n +12i +a 2δx (u n +1i +u n i )−aτ28δx (∂2u ∂t 2)n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t )n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+h 4).(2.14)由于δx (∂2u ∂t 2)n +12i =(∂3u ∂x∂t 2)n +12i+O (h 2),所以可得:δt u n +12i +a 2δx (u n +1i +u n i )−aτ28(∂3u ∂x∂t 2)n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t)n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+τ2h 2+h 4).又因为∂3u ∂x∂t 2=−a ∂3u∂x 2∂t +∂2f ∂x∂t ,所以有:δt u n +12i +a 2δx (u n +1i +u n i )−aτ28(−a ∂3u ∂x 2∂t +∂2f ∂x∂t )n +12i +124(4h 2−a 2τ2)(∂3u ∂x 2∂t )n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i +aτ224(∂2f ∂x∂t)n +12i =f n +12i +O (τ4+τ2h 2+h 4),即,δt u n +12i +a 2δx (u n +1i +u n i )+(a 2τ212+h 26)(∂3u ∂x 2∂t )n +12i −τ224(∂2f ∂t 2)n +12i −h 26(∂2f ∂x 2)n +12i −aτ212(∂2f ∂x∂t )n +12i =f n +12i +O (τ4+τ2h 2+h 4).由于(∂3u ∂x 2∂t )n +12i=δ2x (∂u ∂t )n +12i +O (h 2),所以有:u n +1i −u n i τ+a 4h(u n +1i +1−u n +1i −1+u ni +1−u n i −1)+(h 26+a 2τ212)δ2x u n +1i −u n i τ−τ224(f n +1i −2f n +12i +f n −1i (τ2)2)−h 212[(∂2f ∂x 2)n +1i +(∂2f ∂x 2)n −1i ]−aτ12[(∂f ∂x )n +1i −(∂f ∂x)n −1i ]=f n +12i +O (τ4+τ2h 2+h 4),其中,δ2xu i =u i +1−2u i +u i −1h 2,舍去O (τ4+τ2h 2+h 4),等式两边同时乘以τ,并令λ=τ/h ,整理可得:u n +1i +aλ4(u n +1i +1−u n +1i −1)+(16+a 2λ212)(u n +1i +1−2u n +1i +u n +1i −1)638应用数学2019=u n i−aλ4(u n i +1−u n i −1)+(16+a 2λ212)(u n i +1−2u n i +u ni −1)+τ6(f n +1i −2f n +12i +f n i )+τ12(f n +1i +1−2f n +1i +f n +1i −1+f n i +1−2f n i +f n i −1)+aτλ24(f n +1i +1−f n +1i −1−f n i +1+f ni −1)+τf n +12i,即,(23−a 2λ26)u n +1i +(16+aλ4+a 2λ212)u n +1i +1+(16−aλ4+a 2λ212)u n +1i −1=(23−a 2λ26)u n i +(16−aλ4+a 2λ212)u n i +1+(16+aλ4+a 2λ212)u n i −1+(τ12+aλτ24)f n +1i +1(τ12−aλτ24)f n +1i −1+(τ12−aλτ24)f n i +1+(τ12+aλτ24)f n i −1+2τ3f n +12i .(2.15)由推导过程可知,该格式的截断误差为O (τ4+τ2h 2+h 4),即格式(2.15)在时间和空间上均可达到四阶精度.我们注意到,格式为两层格式,并且格式每层仅用到三个网格点,形成的代数方程组系数矩阵为循环三对角矩阵,可采用追赶法进行求解[8],同时由于要求未知时间层上(第n +1层)中间点的系数不能等于0,即23−a 2λ26=0,因此aλ=2.3.稳定性分析下面采用von Neumann 方法分析本文所推导的差分格式(2.15)的稳定性.对于(2.15)式,舍掉非齐次项f ,即假设f 项精确成立,令u n i =ηn e Iσx i,其中,η为振幅,σ为波数,I =√−1为虚数单位,有(23−a 2λ26)ηn +1e Iσx i +(16+aλ4+a 2λ212)ηn +1e Iσx i +1+(16−aλ4+a 2λ212)ηn +1e Iσx i −1=(23−a 2λ26)ηn e Iσx i +(16−aλ4+a 2λ212)ηn e Iσx i +1+(16+aλ4+a 2λ212)ηn e Iσx i −1.(3.1)两边同时约掉e Iσx i ,并整理可得：(23−a 2λ26)ηn +1+(16+a 2λ212)ηn +1(e Iσh +e −Iσh )+aλ4ηn +1(e Iσh −e −Iσh )=(23−a 2λ26)ηn+(16+a 2λ212)ηn (e Iσh +e −Iσh )−aλ4ηn +1(e Iσh −e −Iσh ).(3.2)利用Euler 公式,即e Iσh =cos σh +I sin σh,e −Iσh =cos σh −I sin σh ,可得:(23−a 2λ26)ηn +1+[(13+a 2λ26)cos σh ]ηn +1+(aλI 2sin σh )ηn +1=(23−a 2λ26)ηn +[(13+a 2λ26)cos σh ]ηn −(aλI 2sin σh )ηn .(3.3)对上式进行化简整理有[(23−a 2λ26)+(13+a 2λ26)cos σh +aλI sin σh 2]ηn +1=[(23−a 2λ26)+(13+a 2λ26)cos σh −aλI sin σh 2]ηn .(3.4)从而可得格式(2.15)的误差放大因子为:G =ηn +1ηn =(23−a 2λ26)+(13+a 2λ26)cos σh −aλI sin σh2(23−a 2λ26)+(13+a 2λ26)cos σh +aλI sin σh2.(3.5)由von Numann 稳定性定理可知当|G |≤1时,格式是稳定的,由(3.5)可得|G |=1,因此,格式(2.15)是无条件稳定的.4.数值实验第3期侯波等：求解一维对流方程的高精度紧致差分格式639为了验证本文格式(2.15)的精确性和稳定性,现考虑以下三个具有精确解的初边值问题.分别采用Crank-Nicolson(C-N)格式,文[7]中格式和本文格式(2.15)进行计算;其中,最大绝对误差及收敛阶的定义为:L∞=maxi |u n i−u(x i,t n)|,Rate=log[L∞(h1)/L∞(h2)]log(h1/h2)L∞(h1)和L∞(h2)为空间网格步长分别为h1和h2时的最大绝对误差.问题1[7]:∂u ∂t +∂u∂x=0,0≤x≤2,t>0,u(x,0)=sin(πx),0≤x≤2,u(0,t)=u(2,t),t>0,该问题的精确解为:u(x,t)=sin[π(x−t)].表1问题1当λ=τ/h=0.5,t=1时刻的最大绝对误差及收敛阶h推进步数(n)C-N格式文[7]本文格式L∞误差Rate L∞误差Rate L∞误差Rate 1/510 2.217(-1) 4.865(-2) 1.993(-3)1/1020 5.752(-2) 1.95 1.263(-2) 1.95 1.208(-4) 4.041/2040 1.450(-2) 1.99 3.199(-3) 1.987.490(-6) 4.011/4080 3.631(-3) 2.008.038(-4) 1.99 4.672(-7) 4.001/801609.082(-4) 2.00 2.014(-4) 2.00 2.919(-8) 4.001/160320 2.271(-4) 2.00 5.041(-5) 2.00 1.824(-9) 4.00表2问题1当τ=λh,t=2时刻的最大绝对误差hτλC-N格式文献[7]本文格式1/160.050000000.8 5.290(-2) 1.292(-2) 1.574(-5) 0.10000000 1.69.013(-2) 5.095(-2) 3.198(-3) 0.20000000 3.2 2.307(-1) 1.941(-1) 6.055(-2) 0.40000000 6.4 6.874(-1) 6.597(-1) 1.746(-2)1/320.025000000.8 1.330(-2) 3.230(-3)9.814(-7) 0.20000000 6.4 2.041(-1) 1.950(-1) 1.575(-3) 0.4000000012.8 6.668(-1) 6.601(-1) 1.916(-2)图1问题1当N=32,τ=0.03125,t=0.2时刻的数值解与精确解640应用数学2019表1给出了针对问题1三种格式在不同空间步长h下,当λ=τ/h=0.5,t=1时的最大绝对误差和收敛阶.我们发现C-N格式在时间和空间上都为二阶精度,由于文[7]格式时间具有二阶精度,空间具有四阶精度,因此当取τ=O(h)时,格式空间仅有二阶精度,而本文格式时间和空间均为四阶精度.图1给出N=32,τ=0.03125,t=0.2数值解与精确解对比图,可以看出数值解与精确解吻合的很好.表2给出了当h=1/16和h=1/32时,τ=λh,t=2时刻对问题1采用三种格式计算的最大绝对误差.可以看出网格比λ最大取到12.8,计算仍然是稳定的,因此本文格式是无条件稳定的.并且本文格式在所有参数下,其计算结果比C-N格式和文[7]格式计算结果更加精确.问题2[7]:∂u ∂t +∂u∂x=0,0≤x≤2,t>0,u(x,0)=e cos(πx),0≤x≤2,u(0,t)=u(2,t),t>0,该问题的精确解为:u(x,t)=e cos[π(x−t)].表3问题2当λ=τ/h=0.5,t=1时刻的最大绝对误差及收敛阶h推进步数(n)C-N格式文[7]本文格式L∞误差Rate L∞误差Rate L∞误差Rate 1/510 6.754(-1) 1.428(-1) 5.567(-2)1/1020 2.310(-1) 1.55 3.099(-2) 2.20 3.041(-3) 4.191/2040 6.027(-2) 1.94 6.825(-3) 2.18 1.904(-4) 4.001/4080 1.492(-2) 2.01 1.658(-3) 2.04 1.165(-5) 4.031/80160 3.705(-3) 2.01 4.115(-4) 2.017.252(-7) 4.011/1603209.250(-4) 2.00 1.028(-4) 2.00 4.527(-8) 4.00表4问题2当τ=λh,t=2时刻的最大绝对误差hτλC-N格式文[7]本文格式1/160.050000000.8 2.171(-1) 5.372(-2) 3.897(-4) 0.10000000 1.6 3.450(-1) 2.056(-1)7.795(-3) 0.20000000 3.2 6.810(-1) 6.111(-1) 3.416(-1) 0.40000000 6.4 1.220 1.198 2.017(-1)1/320.025000000.8 5.575(-2) 1.325(-2) 2.449(-5) 0.20000000 6.4 6.302(-1) 6.109(-1) 2.350(-2) 0.4000000012.8 1.204 1.199 2.201(-1)表3和表4给出了针对问题2利用本文格式和C-N格式以及文[7]格式的计算结果.表3考察了格式的精度,表4验证了格式的稳定性.可以看出本文格式在时间和空间上均可达到四阶精度,并且是无条件稳定的.问题3∂u ∂t +a∂u∂x=f,0≤x≤2,t>0,u(x,0)=cos(πx),0≤x≤2,u(0,t)=u(2,t),t>0,f=π(1−a)sin[π(x−t)],该问题的精确解为:u(x,t)=cos[π(x−t)].第3期侯波等：求解一维对流方程的高精度紧致差分格式641表5问题3当λ=τ/h=0.5,a=0.5,t=1时刻的最大绝对误差及收敛阶h推进步数(n)C-N格式本文格式L∞误差Rate L∞误差Rate1/510 1.124(-1) 4.244(-4)1/1020 3.520(-2) 1.67 2.744(-5) 3.951/20409.957(-3) 1.82 1.739(-6) 3.981/4080 2.551(-3) 1.96 1.134(-7) 3.941/80160 6.413(-4) 1.99 1.351(-8) 3.07问题3为非齐次问题,由于文[7]的方程模型为齐次方程,不能计算非齐次问题,因此该问题我们采用本文格式和C-N进行计算和比较,表5给出了两种格式在不同空间步长h下,当t=1时的最大绝对误差和收敛阶.可以看出当λ=τ/h=0.5,a=0.5时,C-N格式在时间和空间上都为二阶精度,而本文格式时间和空间均为四阶精度.5.结论本文针对一维对流方程提出了一种两层隐式高精度紧致差分格式,时间和空间均采用泰勒级数展开法以及截断误差余项修正法进行处理,格式截断误差为O(τ4+τ2h2+h4),即该格式在时间和空间上均可达到四阶精度.并通过von Neumann方法分析得到该格式为无条件稳定的.最后通过三个数值算例验证了格式的精确性和稳定性.通过上述研究,我们可以得出如下结论:1.文献(如[10-11])中的高精度格式往往是时间多层格式,需要另外构造启动步的计算格式,如果采用低精度格式启动,必然会影响以后时间层的计算精度.而本文格式仅为两层格式,无须启动步的计算,时间即可达到四阶精度.2.文献(如[1,10-11])中的高精度格式空间方向上往往超过三个网格节点,导致靠近边界的内点计算困难,需要采用特殊处理,而本文格式仅用到三个网格节点,可以有效避免这一问题.3.尽管本文格式要求aλ=2,这是本文格式的一个缺陷,但是由于本文格式是无条件稳定的,从理论上讲可以采用任意网格比,因此可以很容易避开aλ=2的条件限制,使得这一缺陷并不太影响格式的使用.4.由于本文方法推导过程中涉及到∂2u∂t2,∂3u∂t3,∂3u∂x3的计算,需要用原方程进行多次求导并进行反复代入计算,在考虑对流项为变系数问题时,将涉及到a(x,t)关于x和t的二阶导数,由于我们考虑在时间半点处,即(x i,t n+12)处的函数值,即要用到(∂2a∂t2)n+12i,如果采用中心差分,则时间仅具有二阶精度,因此本文方法不适用于变系数问题.5.本文方法可直接推广到二维和三维问题中去,我们将另文报道.参考文献：[1]张天德,孙传灼.对流方程的差分格式[J].计算物理,1995,12(2):191-195.[2]于志玲,朱少红.关于对流方程一类三层显格式[J].南开大学学报(自然科学版),1998,31(3):27-30.[3]曾文平.解对流方程的加耗散项的差分格式[J].应用数学,2001,14(S1):154-158.[4]陆金甫,关治.偏微分方程数值解法[M].北京:北京大学出版社,1987.[5]姚朝晖,张锡文,任玉新等.一种低耗散、无伪振荡的实用差分格式[J].水动力学研究与进展(A辑),2001,16(02):195-199.[6]汤寒松,张德良,李椿萱.对流方程保单调CIP格式[J].水动力学研究与进展(A辑),1997(02):181-187.[7]赵飞,蔡志权,葛永斌.一维非定常对流扩散方程的有理型高阶紧致差分公式[J].江西师范大学学报(自然科学版),2014,38(4):413-418.642应用数学2019[8]李青,王能超.解循环三对角线性方程组的追赶法[J].小型微型计算机系统,2002(23):1393-1395.[9]ERDOGAN U.Improved upwind discretization of the advection equation[J].Numer.Meth.PartDiffer.Equ.,2014,30:773-787.[10]BOURCHTEIN A,BOURCHTEIN L.Explicitfinite schemes with extended stability for advectionequations[J]put.Appl.Math.,2012,236:3591-3604.[11]KIM C.Accurate multi-level schemes for advection[J].Int.J.Numer.Methods Fluids.,2003,41:471-494.A High-Order Compact Difference Scheme for Solving the1DConvection EquationHOU Bo,GE Yongbin(School of Mathematics and Statistics,Ningxia University,Yinchuan750021,China)Abstract:In this paper,a two-level implicit compact difference scheme for solving the one-dimensional convection equation is proposed.Taylor series expansion and correction for the third derivative in the truncation error remainder of the central difference scheme are used for the discretization of time and space.The local truncation error of the scheme is O(τ4+τ2h2+h4),i.e.,it has the fourth-order accuracy in both time and space.The unconditional stability is obtained by the von Neumann method. The accuracy and the stability of the present scheme are validated by some numerical experiments.Key words:Convection equation;High accuracy;Compact difference scheme;Unconditional sta-bility;Finite difference method。

一维对流扩散方程是指一维均匀的边界层上的传质过程的数学模型，常用于描述对流扩散过程中的温度、湿度、速度等场的分布情况。

一维对流扩散方程的数学形式为：∂φ/∂t+U∂φ/∂x=D∂^2φ/∂x^2其中φ表示传质物质的浓度，t表示时间，x表示空间坐标，U表示对流速度，D表示扩散系数。

二维对流扩散方程是指二维均匀的边界层上的传质过程的数学模型，常用于描述对流扩散过程中的温度、湿度、速度等场的分布情况。

二维对流扩散方程的数学形式为：∂φ/∂t+U∂φ/∂x+V∂φ/∂y=D∂^2φ/∂x^2+D∂^2φ/∂y^2其中φ表示传质物质的浓度，t表示时间，x和y分别表示两个空间坐标，U和V分别表示两个方向上的对流速度，D表示扩散系数。

单调差分格式是一种常用的数值求解方法，它通过进行差分运算来求解微分方程的数值解。

在求解一维和二维对流扩散方程时，可以使用单调差分格式来解决。

具体来说，可以将空间坐标和时间分别离散化，将对流扩散方程转化为一个线性方程组，然后使用单调差分格式来解决。

单调差分格式的具体形式取决于方程的类型和离散化的方式，但一般来说，它都是将微分方程的差分形式写成一个线性方程组的形式。

例如，在求解一维对流扩散方程时，可以使用下面的单调差分格式：φ_i^{n+1}=φ_i^n+Δt(D(φ_{i+1}^n-2φ_i^n+φ_{i-1}^n)/Δx^2+U(φ_ {i+1}^n-φ_{i-1}^n)/2Δx)其中φ_i^n表示第i个网格点在时间步n的浓度值，Δx和Δt分别表示网格的空间步长和时间步长。

同样的，在求解二维对流扩散方程时，可以使用下面的单调差分格式：φ_i^n=φ_i^n+Δt(D(φ_{i+1,j}^n+φ_{i-1,j}^n+φ_{i,j+1}^n+φ_{i,j-1}^ n-4φ_i^n)/Δx^2+U(φ_{i+1,j}^n-φ_{i-1,j}^n)/2Δx+V(φ_{i,j+1}^n-φ_ {i,j-1}^n)/2Δy)其中φ_i^n表示第(i,j)个网格点在时间步n的浓度值，Δx和Δy分别表示网格在x和y方向上的空间步长，Δt表示时间步长。

一、上机目的用数值方法计算一维对流方程在A 、B 、C 三种差分格式下的解。

x ∆取为0.05. x t ∆∆取值为0.5，1，2。

并作相关讨论。

二、实验原理1x 1 x 00) (x,1x 0 10) (x,0x 1-x 10) (x,0 t 8x 8- ,0>-<=≤≤+-=≤≤+=>≤≤=∂∂+∂∂或ζζζζζx xt三、上机要求：1.学会对MS-FORTRAN 的基本操作。

2.用Fortran 编写程序计算一维对流方程在A 、B 、C 三种格式下的解。

3.讨论各种格式下不同的t x∆∆值的差分格式解的特点。

四、实验程序以A 格式为例，对微分方程进行离散化， 得出 A 格式的差分解的表达式：B 、C 格式同理可以写出。

由此编写如下的Fortran 程序。

注：除了循环时间层的计算公式略有不同外，三个程序没有区别，因此这里只用一个主程序，并根据!空间节点321，dx=0.05 输出依次为（时间，空间，数值）program mainimplicit nonereal dx_dt !定义的值integer abc,r_t,i,j,k !定义变量，abc 为格式类型，r_t 为时间网格数，其余为循环变量real ,allocatable ::s(:,:) !定义存储矩阵swrite (*,*) "输入dx_dt=0.5,1,2"read (*,*) dx_dt write (*,*) "选择格式，A,B,C 分别输入1,2或3"read (*,*) abc!根据格式选择生成相应的文件if (abc==1) thenopen (unit=8,file='out_a.csv')elseif (abc==2) thenopen (unit=8,file='out_b.csv')elseif (abc==3) then open (unit=8,file='out_c.csv')endifr_t=160/dx_dt !计算时间网格总数allocate(s(r_t+1,321)) !分配存储矩阵的空间!第一层赋初值do i=1,140,1s(1,i)=0write(8,*)"1",",",i,",",s(1,i)end dodo i=141,161,1s(1,i)=1+0.05*(i-161)write(8,*)"1",",",i,",",s(1,i)end dodo i=162,181,1s(1,i)=1-0.05*(i-161)write(8,*)"1",",",i,",",s(1,i)end dodo i=182,321,1s(1,i)=0write(8,*)"1",",",i,",",s(1,i)end do!循环时间层,根据格式的选择来判断执行哪一部分if(abc==1) thendo i=2,r_t,1do j=i,322-i,1s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j-1))/(dx_dt*2)write(8,*)i,',',j,',',s(i,j)end dodo k=1,i-1,1 !余下部分赋值0,下同s(i,k)=0write(8,*)i,',',k,',',s(i,k)end dodo k=322-i,321,1s(i,k)=0write(8,*)i,',',k,',',s(i,k)end doend doelseif(abc==2) thendo i=2,r_t+1,1do j=1,322-i,1s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j))/dx_dtwrite(8,*)i,',',j,',',s(i,j)end dodo k=322-i,321,1s(i,k)=0write(8,*)i,',',k,',',s(i,k)end doend doelseif(abc==3) thendo i=2,r_t+1,1do j=i,321,1s(i,j)=s(i-1,j)-(s(i-1,j)-s(i-1,j-1))/dx_dtwrite(8,*)i,',',j,',',s(i,j)end dodo k=1,i-1,1s(i,k)=0write(8,*)i,',',k,',',s(i,k)end doend doendif!完成提示write(*,*)'数据已输出至源目录'pausestopend program五、实验结果及分析程序运行后在对应目录下生成csv表格文件，根据输入的的值不同生成对应的网格并计算各节点数值。

一维扩散方程差分格式的数值计算一维扩散方程是描述物质在一维空间中扩散过程的方程。

数值计算是一种近似求解微分方程的方法，可以通过离散化空间和时间来求解一维扩散方程。

本文将介绍一维扩散方程差分格式的数值计算方法，并给出一个具体的数值计算实例。

∂u/∂t=D∂²u/∂x²其中，u是扩散物质的浓度，t是时间，x是空间坐标，D是扩散系数。

差分格式的基本思想是将连续的时间和空间变量离散化为一系列有限的点，然后用离散化后的点代替原方程中的连续变量，从而得到一个差分方程。

一维扩散方程的差分格式数值计算方法有很多种，下面介绍两种基本的差分格式：显式差分格式和隐式差分格式。

1.显式差分格式：显式差分格式的基本思路是使用当前时间步的解来计算下一个时间步的解。

通过对一维扩散方程进行差分得到：(u_i)_(n+1)=(u_i)_n+D*(∆t/∆x²)*((u_(i-1))_n-2(u_i)_n+(u_(i+1))_n)其中，(u_i)_(n+1)表示时间步n+1时刻、位置i处的扩散物质浓度。

该公式使用当前时间步n的解来逐点计算下一个时间步n+1的解。

2.隐式差分格式：隐式差分格式的基本思路是使用下一个时间步的解来计算当前时间步的解。

通过对一维扩散方程进行差分得到：((u_i)_(n+1)-(u_i)_n)/∆t=D*(∆x²)*((u_(i-1))_(n+1)-2(u_i)_(n+1)+(u_(i+1))_(n+1))这是一个关于时间步n+1的隐式方程，需要使用迭代方法求解。

数值计算的实例：假设在一根长为L的杆上有一种扩散物质，杆的两端固定浓度为0，即u(0, t) = u(L, t) = 0；初始时刻杆上的浓度分布为一个正弦函数，即u(x, 0) = sin(πx/L)；扩散系数为D。

我们需要计算杆上扩散物质的浓度随时间的变化情况。

首先，选择合适的网格间距∆x和时间步长∆t。

然后将杆上的空间坐标和时间离散化为一系列点，得到网格。

求解一维对流方程的差分格式
李海燕;陈豫眉
【期刊名称】《佳木斯大学学报（自然科学版）》
【年(卷),期】2024(42)3
【摘要】针对一维对流方程,提出了一种具有六阶空间精度的差分格式,形成关于时间的半离散化方程;在时间层上,利用指数函数的Pade近似求解该方程.最后通过数值算例验证其精确性.
【总页数】5页(P169-173)
【作者】李海燕;陈豫眉
【作者单位】西华师范大学数学与信息学院;西华师范大学公共数学学院
【正文语种】中文
【中图分类】O241.82
【相关文献】
1.求解一维对流方程的高精度紧致差分格式
2.一种求解三维非稳态对流扩散反应方程的高精度有限差分格式
3.求解一维对流方程的四阶紧致差分格式
4.求解一维非定常对流扩散方程的紧致差分格式
5.求解对流扩散反应方程的高阶指数型组合紧致差分格式
因版权原因，仅展示原文概要，查看原文内容请购买。

致 谢研究生的学习生活就要结束了，而大学最后的一段学习生活也已宣告尾声。

在本文完成之际，我特别向所有曾经给予我帮助、鼓励和支持的人们表示衷心的感谢。

在此，我首先要向我的导师张佐刚教授表示最衷心的感谢。

从入学之初到毕业论文完成的整个过程中，我都一直得到张老师细致入微的关心和照顾，无论是在学习上还是在生活上都给了我莫大的支持和鼓励。

另外，在论文写作过程中，张老师对我悉心的教诲和严格的要求足以让我受益终生，我深深地被老师严谨的治学精神和一丝不苟的生活态度所折服，他一丝不苟的育人作风更是为我以后的学习工作树立了榜样。

这里我还要特别感谢高雷阜教授在平时一直给予的默默关心，在导师繁忙的时候，给予了我无私的帮助，同时也是他为我们提供了一个方便优越的学习环境，使论文得以顺利完成。

最后还要感谢辽宁工程技术大学理学院所有的老师，是他们在传授我知识的同时培养了我良好的数学思维和数学的思考方式，还有理学院所有的同学，他们在日常的学习生活中给了我无数的帮助和启发。

祝福所有的老师和同学身体健康，工作和学业一切顺利！摘 要对流扩散方程作为流体力学中的基本方程之一，目前求解的数值方法有很多，主要包括有限单元法、有限体积法、有限分析法和有限差分法。

其中，有限差分法作为一种重要的数值离散方法，在科学研究和工程计算中都得到了广泛的应用。

而在有限差分法中，高精度有限差分方法又以其涉及网格点少、边界无需特殊处理、具有较高的计算精度等优点，成为学者们争相研究的热点问题。

本文首先分析了高精度差分格式的研究意义、国内外的研究现状、研究过程中存在的一系列问题和发展的趋势。

然后结合均匀网格上和非均匀网格上差分算子的定义式，具体介绍了一些传统的关于对流扩散方程的差分格式和目前已经得到一定发展的主要的高精度差分格式。

通过对这些差分格式的研究学习，了解了很多高精度差分格式的主要构造方法，从中也得到了一些启发。

并且通过对一些差分方法的优点和局限性的分析，改进了文献中的一些高精度差分格式的构造方法，分别应用待定系数法、降维法等数学思想，给出了均匀网格上简单对流方程的一种高精度差分格式，以及非均匀网格上二维对流扩散反应方程的一种高精度紧致差分格式，并给出了相应的数值算例，证明了该格式较好的计算效果。

一维定常对流扩散反应方程的高精度紧致差分格式祁应楠;武莉莉【摘要】针对一维定常对流扩散反应方程,提出了一种四阶精度的有理型紧致差分格式,其局部截断误差为O(h4);然后通过Richardson外推技术和算子插值法将本文格式的精度提高到六阶.因为格式仅涉及到3个网格基架点,所以对于Dirichlet 边值问题,由差分格式可得三对角线性方程组,可采用追赶法进行求解.最后通过数值算例验证了本文方法的精确性和可靠性.【期刊名称】《华中师范大学学报（自然科学版）》【年(卷),期】2017(051)001【总页数】6页(P1-6)【关键词】对流扩散反应方程;高阶紧致格式;Richardson外推;有限差分法【作者】祁应楠;武莉莉【作者单位】宁夏师范学院数学与计算机科学学院,宁夏固原756000;宁夏师范学院数学与计算机科学学院,宁夏固原756000【正文语种】中文【中图分类】O241.8对流扩散反应问题是流体力学、传热学、传质学等学科以及环境、化工等应用领域中经常遇到的典型问题之一，由于问题的准确解往往很难获得，所以人们经常采用数值方法来寻求问题的近似解.目前，所流行的近似计算方法包括有限差分法、有限元法和边界元法等.其中有限差分方法是一种常用的数值计算方法.目前，国内外已经有许多有关该问题高阶紧致差分格式的研究报道.如：魏剑英[1]针对一维对流扩散方程，提出了一种指数型高阶紧致差分格式.王彩华[2]利用泰勒展开公式和数项级数收敛性给出了一线性对流扩散问题的一类高精度紧致差分格式.田芳和田振夫[3]基于非均匀网格上函数的泰勒级数展开，构造了非均匀网格上的高精度紧致差分格式.Sun和Zhang[4]构造了定常对流扩散反应方程的多项式型四阶紧致差分格式，并用Richardson外推法[5]和算子插值技术将格式的精度提高到了六阶.Tian和Dai [6]构造对流扩散问题的指数型格式，其空间具有四阶精度.文献[7]研究了非定常对流扩散方程的有理型高阶紧致差分格式并得到了很好的计算效果.杨志峰等[8]构造了含源项非定常对流扩散问题的紧致四阶格式.文献[9-11]研究了利用样条插值的方法来构造高精度紧致差分格式.文献[12]通过消除对流项，并利用Pade格式，构造了一维非定常对流扩散反应方程无条件稳定的四阶紧致差分格式.文献[13]针对非定常对流扩散方程，对空间采用三点紧致差分格式，并对时间采用单对角隐式Runge-Kutta方法进行离散，得到了截断误差为O(τ4+h4)的无条件稳定的隐格式.文献[14]通过简单的分裂算法及增加特殊网格点的方法，对时间的处理采用C-N格式与向后欧拉结合的技巧，推导出求解高维非定常对流扩散反应方程的隐式差分格式.本文针对一维定常对流扩散反应方程，基于截断误差余项修正思想，并结合原方程本身，推导得到了求解该方程的一种四阶精度的有理型紧致差分格式.然后采用Richardson外推法和算子插值技术将格式的精度提高到六阶.最后给出了数值算例. 本文讨论的方程模型为两点边值问题：其中，边界条件为：u(0)=q0，u(L)=qL.这里，a，p(x)，b(x)分别为扩散、对流和反应项系数.且a>0，p(x)和b(x)均为关于空间变量x的光滑函数.首先将定义域[0，L]离散为：0=x0<x1<x2<…<xN=L，记为网格数，空间步长定义为.为了方便书写，用ui表示表示，依此类推.将式(1)改写为：由此定义空间一阶和二阶导数的中心差分算子为：将式(1)利用中心差分代替，并利用关于u的一阶数和二阶导数的定义，可得：).如果略去i和i项，及高阶截断误差项O(h4)，即可得到二阶中心差分格式.为了得到高精度的差分格式，将式(5)中的三阶和四阶导数项进行处理，为此对式(2)两边同时关于x求一阶导和二阶导数可得：.将式(6)代入式(7)消去xi化简可得：.将式(6)和式(8)代入式(5)可得：fi+O(h4).将式(3)和式(4)代入式(9)，略去高阶项后化简整理可得：其中式(10)即为多项式型四阶紧致(FOC)差分格式，此格式色散误差和耗散误差较大.为了能精确数值求解此类方程，我们推导一种有理型的四阶紧致差分格式.将式(2)代入式(6)可得：.在式(15)中令，则式(15)可简化为：将式(2)和式(15)代入式(7)，整理可得：.令则式(17)可化简为：将式(16)和式(18)代入式(5)，整理可得：其中，.接着，将式(19)中的一阶导数用中心差分代替加×(5)，利用式(3)和(4)对u的一阶和二阶导数进行离散，整理可得：).将式(15)和式(17)代入式(20)加×(5)，同样利用u的一阶和二阶导数进行离散，整理可得：).令则式(21)可化简为：(γ4f+γ5fx+γ6fxx)i+O(h4).然后，略去式(22)中i和i项，并将式(3)和式(4)代入式(22)，化简整理可得有理型的高阶紧致差分格式：，其中，.此格式的高阶截断误差为O(h4)，即此格式具有四阶精度.本文格式之所以称之为有理型格式，是因为其差分算子的系数为有理型函数，记为RHOC.从推导过程可以看出，FOC格式只是其中的一种特殊情况，有理型格式的推导更具有广泛性.结合原方程可得到具有不同性质的高精度格式，对于不同性质的问题可选用与之相适应的格式进行求解，此类格式均为三个网格基架点，只发生系数的变化.下面使用Richardson外推方法[5]将本文的四阶格式RHOC提高到六阶精度.◆………… ◆——◆0 1 2 ……N/2◆—●—◆…… ◆—●—◆0 1 2 3 4……N-1 N首先，用格式(23)分别在粗网格和细网格上计算一遍可得粗网格和细网格具有四阶精度的近似解.然后利用Richardson外推公式：可得到在粗网格具有六阶精度的解，其中和分别为粗网格和细网格上的解，即通过Richardson外推公式(24)可将细网格上偶数点(菱形点)的精度提高到六阶.为了使细网格上奇数点(圆点)的精度也达到六阶，使用算子插值法.由式(23)可得：).由于细网格上偶数点(菱形点)已经算出，因此只须采用式(26)计算奇数点(圆点)，即可得如下算子插值公式：.通过式(26)利用细网格上具有六阶精度的偶数点来计算奇数点，从而可使得细网格上点的精度均为六阶，整个过程我们将其记为RRHOC，其算法步骤如下：1) 在粗网格上用格式(23)计算一遍，可得；2) 在细网格上用格式(23)计算一遍，可得；3) 通过外推式(24)将细网格上偶数点(菱形点)的精度计算至六阶，4) 通过式(26)将细网格上奇数点(圆点)的精度计算至六阶.为了验证本文格式的精确性和可靠性，分别采用RHOC格式和RRHOC格式对以下两个有精确解的问题进行数值实验，并与中心差分格式、多项式型四阶紧致格式(FOC) [4]和六阶格式(REC) [4]的计算结果进行比较.其中，L∞范数误差和收敛阶(Rate)的定义如下：L∞范数误差，其中，Ui表示点xi处的精确解，ui表示点xi处的数值解，L∞(uh1)和L∞(uh2)分别表示网格步长为h1和h2时对应的L∞范数误差．问题1：该问题的精确解为：u(x)=ex.取：a=1，b(x)=x2+1，f(x)=x2ex.问题2：该问题的精确解为：u(x)=e-4πsin(x).取：a=1，p(x)=1，b(x)=1，f(x)=e-4π(cos x+2sin x).对于问题1和问题2，表1和表3列出了取不同步长h时，采用中心差分格式、FOC格式[4]与本文RHOC格式计算的L∞范数误差和Rate(收敛阶).不难得到，本文所提的四阶精度的有理型格式(RHOC)格式比多项式型格式(FOC) 和中心差分格式均具有更高的准确度.而且，当网格数不断增加时，RHOC格式的L∞范数误差比中心差分格式小四个数量级不等，比同是四阶的FOC格式计算结果更精确.表2和表4列出了取不同网格步长h时，REC格式[4]与本文RHOC格式的最大绝对误差和收敛阶，从表中可以看出经过外推和算子插值之后的REC格式和本文RHOC格式均有六阶精度，但是本文RHOC格式的计算误差明显优于REC格式[4].本文基于中心差分格式的截断误差余项修正，并利用原方程本身，提出了数值求解一维两点边值问题的一种紧致的高精度差分方法，由理论推导可知所提格式为四阶精度.然后采用Richardson外推法和算子插值技术将格式的精度提高到六阶.最后，采用本文两种方法计算了两个数值算例，并与传统的中心差分格式以及文献[4]中的FOC格式和REC格式进行了对比，充分体现了本文方法的精确性和有效性. 【相关文献】[1] 魏剑英. 定常对流扩散反应方程的指数型高阶差分格式[J].宁夏大学学报(自然科学版)， 2012，33(2):140-143．[2] 王彩华. 一维对流扩散方程的一类新型高精度紧致差分格式[J].水动力学研究与进展， 2004，19(5):655-663.[3] 田芳，田振夫. 定常对流扩散反应方程非均匀网格上的高精度紧致差分格式[J].宁夏师范学院学报(自然科学版)， 2009， 26(2):219-225．[4] SUN H， ZHANG J. A High order finite difference discretization strategy based on extrapolation for convection diffusion equations[J].Numer Methods Partial Differential Eq，2004， 20(1):18-32[5] CHENEY W， KINCARD D. Numerical Mathematics and Computing[M]. 4th Ed. CA：Brooks/Cole Publishing，Pacific Grove，CA. 1999．[6] TIAN Z F， DAI S Q. High-order compact exponential finite difference methods for convection-diffusion type problems[J]. J Comput Phys， 2007， 220:952-974．[7] 赵飞，蔡志权，葛永斌. 一维非定常对流扩散方程的有理型高阶紧致差分格式[J].江西师范大学学报(自然科学版)， 2014， 38(4):413-418．[8] 杨志峰，陈国谦. 含源项非定常对流扩散问题紧致四阶差分格式[J].科学通报， 1993，38(2):113-116．[9] LANG F G， XU X P. Quintic B-spline collocation method for second order mixed boundary value problem[J]. Computer Physics Communications. 2012， 183:913-921．[10] GOH J， MAJID A A, ISMAIL A I M. A quartic B-spline for second-order singular boundary value problems[J]. Computers and Mathematics with Applications. 2012，64:115-120．[11] 林建国，许维德，陶尧森.含源项非定常非线性对流扩散方程的三次样条四阶差分格式[J].水动力学研究与进展(A辑)， 1994， 9(2):599-602．[12] 杨录峰，李春光.一种求解对流扩散反应方程的高阶紧致差分格式[J]. 宁夏大学学报(自然科学版)， 2013， 34(2): 101-104．[13] PIAO X， CHOI H J， KIM S D， et al. A fast singly diagonally implicit Runge-Kutta method for solving 1D unsteady convection-diffusion equations[J]. Numercal Methods for Partial Differential Equations，2013(30):788-812．[14] ADEM K. Finite difference approximations of multidimensional unsteady convection-diffusion-reaction equations[J]. Journal of Computational Physics， 2015， 285:331-349．。

一维非定常对流扩散方程非均匀网格上的高精度紧致差分格式黄雪芳;郭锐;葛永斌【摘要】A high accuracy compact finite difference scheme with non-uniform grids is pro-posed to solve unsteady convection diffusion equations, which are used to describe boundary layer problems or locally large gradient problems, etc. The new method starts from the dis-cretization of the steady convection diffusion equation. Firstly, the spatial derivatives are discretized by using the Taylor series expansion on non-uniform grids. Then, the second order backward Eulerian difference formula is used to discretize the temporal derivative term. The three-level full implicit compact difference scheme on non-uniform grids for solving the one-dimensional unsteady convection diffusion equation is derived. The new scheme has the second order accuracy in time and the third to fourth order accuracy in space and is unconditionally stable. Finally, some numerical experiments are conducted to demonstrate the high accuracy and the advantages in solving boundary layer problems or locally large gradient problems.%本文在非均匀网格上给出了求解非定常对流扩散方程的一种高精度紧致差分格式，特别适合边界层和大梯度等问题的求解。

一维扩散方程差分格式的数值计算∂u/∂t=D∂²u/∂x²其中，u(x,t)是在位置x和时间t的扩散现象的浓度，D是扩散系数。

为了对一维扩散方程进行数值计算，可以使用差分格式。

最常用的差分格式是向前差分和中心差分。

1.向前差分格式：使用向前差分格式将时间t和位置x分别离散化，差分步长分别为Δt和Δx。

将扩散方程中的偏导数用有限差分近似替代，可以得到近似方程：(u_i(t+Δt)-u_i(t))/Δt=D(u_i-1(t)-2u_i(t)+u_i+1(t))/Δx²其中，u_i(t)表示在位置x_i和时间t的解，u_i(t+Δt)和u_i(t)是上一时刻和当前时刻的浓度，u_i-1(t)和u_i+1(t)分别是x_i左右两侧位置的解。

这样，一维扩散方程就被转化为一个差分方程。

根据初始条件u(x,0)和边界条件u(0,t)和u(L,t)，L表示空间区域的长度，可以得到差分方程的初始条件。

使用向前差分格式可以得到一个显式迭代公式：u_i(t+Δt)=u_i(t)+DΔt(u_i-1(t)-2u_i(t)+u_i+1(t))/Δx²这个公式可以用来逐步推进时间t的步骤，从而获得扩散过程中的浓度分布。

2.中心差分格式：使用中心差分格式将时间t和位置x分别离散化，差分步长分别为Δt和Δx。

将扩散方程中的偏导数用有限差分近似替代，可以得到近似方程：(u_i(t+Δt)-u_i(t))/Δt=D(u_i-1(t)-2u_i(t)+u_i+1(t))/Δx²与向前差分格式不同的是，在右侧位置x_i+1处使用u_i+1(t)近似。

这个差分方程可以进一步简化为一个稳定的隐式迭代公式：u_i(t+Δt)=u_i(t)+DΔt(u_i-1(t+Δt)-2u_i(t+Δt)+u_i+1(t+Δt))/Δx²这个公式可以通过求解线性方程组来计算下一个时间步长的解。

以上是一维扩散方程差分格式的数值计算的基本原理和方法。

一维对流方程的差分格式∂u/∂t+c∂u/∂x=0其中u是流体密度，t是时间，x是空间位置，c是流体速度。

为了离散化这个方程，我们将时间和空间范围分为小区间，然后在这些离散点上近似连续方程。

我们可以使用中心差分、向前差分或向后差分等方法来近似对流方程。

中心差分格式是一种常用的差分格式，通过使用当前时间步和两个相邻时间步的值来近似时间导数，使用当前空间点和两个相邻空间点的值来近似空间导数。

假设我们在时间方向上将时间分为n个步长Δt，将空间方向上将空间分为m个步长Δx，那么时间和空间步长分别为Δt = t_max / n 和Δx = x_max / m。

在中心差分格式中，时间导数可以使用向前差分或向后差分来计算。

使用向前差分有：∂u/∂t≈(u_i,j-u_i-1,j)/Δt使用向后差分有：∂u/∂t≈(u_i+1,j-u_i,j)/Δt空间导数可以使用中心差分来计算：∂u/∂x≈(u_i+1,j-u_i-1,j)/(2Δx)结合时间导数和空间导数的近似，我们可以得到中心差分格式的一般形式：(u_i,j+1-u_i,j)/Δt+c(u_i+1,j-u_i-1,j)/(2Δx)=0通过对该方程进行变形，可以得到u_i,j+1的计算公式：u_i,j+1=u_i,j-cΔt/(2Δx)(u_i+1,j-u_i-1,j)在空间和时间方向上进行迭代，可以逐步求解一维对流方程。

除了中心差分，还存在其他差分格式，如向前差分和向后差分。

向前差分格式使用当前时间和前一时间步的值来近似时间导数，向后差分格式使用当前时间和后一时间步的值来近似时间导数。

中心差分格式的优点是它具有二阶精度，即误差随着步长的平方而减少。

然而，该格式可能会出现数值耗散或数值扩散的问题，特别是在高梯度区域。

在实际应用中，选择正确的差分格式对于获得准确的数值解至关重要。

一维对流方程的差分格式提供了一种近似求解连续方程的方法，使得我们能够使用计算机来解决复杂的流体力学问题。