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

合集下载

一维对流扩散方程的4种差分格式的Jacobi迭代收敛性比较

（．１宁夏大学数学计算机学院，宁夏银川７０２；２复旦大学数学科学学院，海５０１．上
摘
要：对一维常系数对流扩散方程采用不同的差分格式离散后所得到的线性系统，过直接估计Ｊｃｂ迭代针通ａｏｉ
考虑一维对流扩散模型方程：
～
１谱半径
（）１
＋“ 一ｆ（￡，ｘ，）
当系数￡￡）小时，（＞０很该方程常作为奇异摄动问题
定义ｙ／２）网格雷诺数，ｙ１时，＝（ｅ为当＞称
的模型方程．而此时，于该方程，典的中心差分对古格式（称为Ｃ）用Ｊｃｂ迭代方法求解会不简ＤＳ采ａｏｉ
的数值震荡，但其数值精度只有一阶．文献［］１中四阶紧致差分格式（简称为ＨＯｓ通过添加人工黏性ｃ）项抑制了迭代的数值震荡，随着扩散系数的减小但其收敛变得很慢，而文献［］２中指数型高阶紧致差分格式（简称为ＥＣ）０ＨＯＳ￣在集中了以上３格式的高精种度、高稳定性优点的同时克服了数值上的震荡，不但具
为阶三对角矩阵，则存在实的非奇异对角矩阵Ｑ，使得ＱＱ为实对称矩阵，且仅当对每个ｉ１Ａ当（ ≤ ｉ）有ｂＣ ≤ ，ｉ＞０或ｂ一Ｃ一０所得对称矩阵为 … ．
Ａ＝ｔｉ一１７２ｙ一１．＝ｒ＝［ —２，＋２，］定理２对任何网格雷诺数ｙ点Ｊｃｂ迭代矩，ａｏｉ阵的谱半径为

求解一维扩散反应方程的隐式高精度紧致差分格式

求解一维扩散反应方程的隐式高精度紧致差分格式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.。

一维对流方程在三种差分格式下的解

一维对流方程在三种差分格式下的数值解
一、实验目的
用数值方法计算一维对流方程在 A、B、C 三种差分格式下的解。∆x 取为 0.05，∆∆��取值
为 0.5,1,2.并作相关讨论
�� + �� = 0,
IF(I==1)THEN !边界点置为相邻点的值 U(I,J)=U(I+1,J-1)
ELSE IF(I==XN) then U(I,J)=U(I-1,J-1)
ELSE U(I,J)=U(I,J-1)-GAMMA*(U(I+1,J-1)-U(I-1,J-1))/2
END IF WRITE(9,*) X(I),U(I,J) END DO END DO WRITE(*,*)"***********Finish**********" write(*,*)"请在程序文件夹下打开Output_FormA.txt查看结果" return end subroutine
OPEN(UNIT=10,FILE="data.txt",STATUS='OLD',ACTION='READ') READ(10,"(a)")TEMP !标题信息 READ(10,*) TN,DX,GAMMA,LB,RB,FORM !读取6个参数 WRITE(*,*) "Read successfully" ELSE WRITE(*,*)"The document doesn't exist." STOP END IF !计算空间网格数 XN=(RB-LB)/DX+1 !确定U矩阵和X向量的大小 ALLOCATE(U(XN,TN+1)) ALLOCATE(X(XN)) !定义U初值为0 U=0 !计算t=0时，U的函数值U（:,1） DO I=1,XN X(I)=LB+(I-1)*DX IF(X(I)>=-1 .AND. X(I)<=0) THEN

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

应用数学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.PartDiﬀer.Equ.,2014,30:773-787.[10]BOURCHTEIN A,BOURCHTEIN L.Explicitﬁnite 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 Diﬀerence 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 diﬀerence 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 diﬀerence 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 diﬀerence scheme;Unconditional sta-bility;Finite diﬀerence 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三种差分格式Word版

一、上机目的用数值方法计算一维对流方程在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表格文件，根据输入的的值不同生成对应的网格并计算各节点数值。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

程的Ｏ（＋ｈ）格式，但该格式对齐次边界合适．
文献［７］中对空间变量四阶紧致离散，时间应用三
① 收稿日期：２０１３—１２—２５基金项目：新疆大学本科生创新项目（ｘＪｕ— ＳＲＴ一１３０１７）．
步长，网格函数ｔｔ（ｘｉ，ｔ）己作
解，目前对该问题主要差分格式有中心显式差分格式，Ｄｕｆｏｒｔ —Ｆｒａｎｋｅｌ差分格式，指数型差分格式，Ｓａｍａｒｓｋｉ差分格式，Ｃｒａｎｋ— Ｎｉｃｈｏｌｓｏｎ格式等¨ Ｊ，
１本文差分格式的构造
本文考虑如下对流扩散问题的初边界问题：
文［３］中对流扩散方程四中常用差分格式进行比
较和各种差分格式的适用范围进行了详细的讨论．文献［４］中用时间离散方法对对流扩散方程求解
得到了空间方向具有７阶精度的差分格式，但是在时间方向具有２阶精度．文献［５］中利用修正局部Ｃｒａｎｋ—Ｎｉｅｈｏｌｓｏｎ方法对问题进行求解得到对流扩散方程的绝对稳定的显式差分格式并稳定性进行
ＶｏＩ１４
文章编号：１００８—１４０２（２０１４）０１ —０１３５— ０４
求解一维对流扩散方程的高精度紧致差分格式①
开依沙尔・热合曼，阿孜古丽．牙生，祖丽皮耶．如孜
方面，对流扩散方程也有其独特的特点，在自然科学中很多现象是用对流扩散方程或方程组描述的，
次Ｃ样条配置方法得到了对流扩散方程的Ｄ（＋ｈ）格式．
本文利用半离散技术对空间变量四阶紧致格式进行离散，时间变量保持不变，把一维对流扩散方程转化为常微分方程组的初值问题，再利用梯形方法构造对流扩散方程的时间二阶空间四阶精
ｄ％ｏ％
（２）
４，Ｂ是具有如下形式的（ｍ—１）阶三对角矩阵
ａ１ａ２０
口３ｏｌ口２
…
０
０
得到一个求解方程（２）的三点四阶紧致差分
格式如下：
一
Ｏ
… Ｏ
０
（＋１２
）
＋＿
＝（１＋ｈ２（２
（新疆大学数学与系统科学学院。新疆乌鲁木齐８３００．４６）
摘
要：对空间变量四阶紧致格式进行离散，时间变量保持不变，把一维对流扩散方程转化为
常微分方程组的初值问题，再利用梯形方法构造对流扩散方程的时间二阶空间四阶精度的一种差分格式，并稳定性进行分析，数值结果与Ｃｒａｎｋ—Ｎｉｃｈｏｌｓｏｎ格式进行比较，数值结果表明，该方
度的绝对稳定的隐式差分格式．
本文构造的计算方法有较好的数值稳定性和较高的计算精度，用来求解对流占优的扩散问题时克服了通常数值解法所遇到的困难．数值结果表明，该方法可以很好地解决对流扩散问题的数值计算．
方程的准确解很难得到的，因此大多数情况用数值
方法求解．本文考虑采用差分格式对问题进行求
法可以很好地解决对流扩散方程的数值计算
关键词：对流扩散方程；高精度紧致差分格式；梯形公式；Ｃｒａｎｋ—Ｎｉｃｏｌｓｏｎ格式
中图分类号：０２４１．８２文献标识码：Ａ
引言
对流扩散方程是一类在物理上有着重要意义的抛物型方程，通常我们将对流方程和扩散方程的差分方法结合起来就可以得到对流扩散方程的差分方法，当扩散系数极小时，反映对流方程的特征，当对流系数极小时，反映扩散方程的特征，但另一
首先考虑定常对流扩散方程
作者简介：开依沙尔・热合曼（１９７８一）．男，新疆库车人，硕士，讲师．主要从事微分方程数值解法的研究．
１３６
佳木斯大学学报（自然科学版）
２０１４年
一
＋卢＝厂（）
第３２卷第１期
２０１４年０１月
佳木斯大学学报（自然科学版）ＪｏｕｎａｒｌｏｆＪｉａｍｕｓｉＵｎｉｖｅｒｓｉｔｙ（ＮａｔｕｒａｌＳｃｉｅｎｃｅＥｄｉｔｉｏｎ）
（１）将求解区域作剖分，区间［ａ，ｂ］作ｍ等分，将
详细讨论．文献［６］中对空间变量四阶紧致离散，时间应用指数函数的Ｐａｄｅ公式得到了对流扩散方
区间［Ｏ，Ｔ］作ｎ等分，并记＾＝
，７－：，：
ｉｈ，０≤ｉ ≤ｍ，ｔ＝．ｒ．称ｈ和丁为空间步长和时问
ｆ＋卢＝瑶，：［口 ≤ ≤ ６］ × ［０ ≤ ￡ ≤ ７１］
１ｕ（ｘ，０）＝），ｏ≤ ≤ｂ
ｔｕ（ａ，ｔ）：ｇ１（￡），／ｄ，（ｂ，）＝ｇ２（￡），０≤ ｔ ≤Ｔ