下载文档原格式
求解一维扩散反应方程的隐式高精度紧致差分格式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隐式高精度紧致差分格式在上述差分方法中，我们采用了一阶差分法和二阶差分法，这种方法的精度有限。
为了提高求解的精度，可以采用更高阶的差分方法。
对流扩散方程的数值方法流扩散方程是描述物质在流动中同时进行的扩散过程的方程。
在很多科学和工程领域，如物理、化学、生物学等，流扩散方程都具有重要的应用。
为了解决流扩散方程，在数值计算中可以采用不同的数值方法。
本文将介绍几种常用的数值方法，包括有限差分法、有限元法和谱方法等。
有限差分法是一种常用且简单的数值方法，可用于解决流扩散方程。
它将空间和时间离散化，并采用中心差分近似来计算偏导数。
通过将方程离散化为代数方程组，可以使用迭代方法（如雅可比方法、高斯-赛德尔方法等）求解。
有限差分法的主要优点是简单易行，且可以方便地处理复杂的边界条件。
然而，它在处理不规则边界和复杂的时间变化时可能会出现精度问题。
有限元法是一种更加灵活和通用的数值方法，可用于解决流扩散方程。
它将连续的空间和时间域划分为离散的小单元，并利用有限元近似来计算解。
有限元法的优点是适用于各种不规则边界和复杂的几何结构，且能够提供更高的精度。
它通常使用高阶基函数来提高数值解的精确度，但计算复杂度较高，并且需要额外的后处理步骤来获得所需的物理量。
谱方法是一种基于傅里叶级数和函数的展开来计算数值解的方法，也适用于解决流扩散方程。
它使用特殊类的基函数（如傅里叶基函数或Chebyshev基函数）来表示解，并利用傅里叶级数的收敛性和高精度的性质来求解偏微分方程。
谱方法的优点是能够提供非常高的精度，并且适用于各种边界条件和几何结构。
但是，谱方法通常对于非线性问题的数值求解比较困难，且需要合适的扩展性来处理大规模问题。
对于流扩散方程的数值方法，除了上述几种常见的方法外，还有其他一些方法如交替方向隐式方法（ADI方法）和双曲正切方法（双曲正切线性增量法）等。
这些方法在特定情况下可能更适用于一些问题，但在一般情况下，有限差分法、有限元法和谱方法是流扩散方程数值计算的主要选择。
在选择数值方法时，需要综合考虑问题的特点和要求。
有限差分法适用于简单的几何结构和边界条件，有限元法适用于复杂的几何结构和边界条件，谱方法适用于需要高精度和快速收敛的问题。
应用数学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 ，它是描述一种随时间变化的物理过程，即所谓不定常现象。
这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。
本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。
1.扩散方程一维扩散方程为：22u u t xα∂∂=∂∂ （1）式中，u 为因知量，α为扩散系数，x 为坐标，t 为时间。
其定解条件如下： 初始条件： (,0)() 0x u x f x L =≤≤（2）边界条件： 12(0,)() , (,)()u t f t u L t f t ==（3） 一般假定函数()f x ，1()f t ，2()f t 满足连接条件，即1(0)(0) f f =，2()(0) f L f =。
2.有限差分法有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。
其控制方程中的导数用离散点上函数值的差商代替。
差分格式可以分为显格式和隐格式。
所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。
由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。
隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。
因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。
为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。
因此，隐格式总是伴随着求解巨大的代数方程组。
Keller-Segel趋化模型的高精度紧致差分方法专业品质权威编制人：______________审核人：______________审批人：______________编制单位：____________编制时间：____________序言下载提示：该文档是本团队精心编制而成，期望大家下载或复制使用后，能够解决实际问题。
文档全文可编辑，以便您下载后可定制修改，请依据实际需要进行调整和使用，感谢!同时，本团队为大家提供各种类型的经典资料，如办公资料、职场资料、生活资料、进修资料、教室资料、阅读资料、知识资料、党建资料、教育资料、其他资料等等，想进修、参考、使用不同格式和写法的资料，敬请关注!Download tips: This document is carefully compiled by this editor. I hope that after you download it, it can help you solve practical problems. The document can be customized and modified after downloading, please adjust and use it according to actual needs, thank you!And, this store provides various types of classic materials for everyone, such as office materials, workplace materials, lifestyle materials, learning materials, classroom materials, reading materials, knowledge materials, party building materials, educational materials, other materials, etc. If you want to learn about different data formats and writing methods, please pay attention!Keller-Segel趋化模型的高精度紧致差分方法Keller-Segel趋化模型的高精度紧致差分方法引言Keller-Segel模型是一种用于描述生物群体中化学物质扩散和趋化运动的数学模型。
求解变系数对流扩散反应方程的指数型高精度紧致差分方法田芳;葛永斌【摘要】本文给出了一种数值求解变系数对流扩散反应方程的指数型高精度紧致差分方法.我们首先将模型方程变形,借助常系数对流扩散方程的指数型高精度紧致差分格式,采用残量修正法得到变系数对流扩散反应方程的指数型高精度紧致差分格式;并从理论上分析了当Pelect数很大时,本文格式达到四阶计算精度时网格步长的限制条件;离散得到的代数方程组可采用追赶法直接求解.数值实验结果与理论分析完全吻合,表明了本文格式对于边界层问题或大梯度变化的物理量求解问题具有的高精度和鲁棒性的优点.%An exponential high accuracy compact finite difference method is proposed to solve the one-dimension (1D) convection-diffusion-reaction equation with variable coefficients. Fir-stly, the equation is rewritten in the form of convection diffusion equation. Then the exponential high order compact finite difference scheme for the convection diffusion equation with constant coefficients and the remainder term modification approach are utilized to obtain an exponential high accuracy compact finite difference scheme for the 1D convection-diffusion-reaction equation with variable coefficients. Secondly, the necessary condition on grid step length is analyzed theoretically if the scheme in this paper has a fourth-order accuracy when the Peclet number is very high. Lastly, the Thomas approach is applied to deal with the algebraic equations. Numerical examples, mostly with the boundary layer where sharp gradients may appear due to high Peclet number, arepresented to demonstrate the accuracy and robustness of the proposed scheme.【期刊名称】《工程数学学报》【年(卷),期】2017(034)003【总页数】14页(P283-296)【关键词】对流扩散反应方程;指数型有限差分格式;高精度紧致差分格式;对流占优;边界层【作者】田芳;葛永斌【作者单位】宁夏大学数学统计学院,银川 750021;宁夏大学数学统计学院,银川750021【正文语种】中文【中图分类】O241.821 引言对流扩散反应方程是自然界和工程应用中的一类非常重要的方程，常常用来描述大气、海洋、河流等污染中污染物的扩散与分布、细菌的浓度分布、核工业中核反应堆的冷却及工业生产中的化学气相沉积等对流扩散反应现象．因此，对于此类方程的数值求解研究具有重要的理论意义和实际应用价值．目前，求解对流扩散反应方程的方法主要有有限元法、有限体积法、边界元法、特征线方法、有限差分法等[1-18]．文献[11]通过网格加密技术与变分多尺度法的耦合，消除了在对流扩散反应方程中由边界层效应和内部层效应引起的数值伪振荡．文献[12]利用微分方程的通解，提出了数值求解一维定常常系数对流扩散反应方程的一种指数型高精度差分方法；文献[13]分别针对对流扩散方程和对流扩散反应方程采用变量替换消去方程中的对流项，将方程转化为反应扩散方程组进行离散求解；文献[14,15]构造了非均匀网格上的对流扩散反应方程的多项式型的紧致差分格式，该格式较均匀网格系统下的计算无论在精度上还是分辨率上均有很明显的优势；文献[16]基于泰勒级数展开，结合原方程，发展了求解对流扩散方程的多项式型四阶精度的紧致差分格式；文献[17]发展了求解对流扩散反应方程的多项式型四阶精度的紧致差分格式，并将该格式外推得到多项式型六阶精度的紧致差分格式；文献[18]采用差分修正的思想构造了求解对流扩散方程的指数型差分格式，该格式对于对流占优问题在粗网格上能得到很高的计算精度．本文将针对对流扩散反应模型方程，将文献[18]的差分格式进行推广，发展了求解对流扩散反应方程的指数型高精度紧致差分格式．2 对流扩散反应方程的指数型高精度紧致差分格式考虑定常对流扩散反应方程其中a是扩散系数，一般是常数；c(x)和b(x)分别为对流项系数和反应项系数，可以是常数，亦可以是x的函数；f(x)是x的足够光滑的函数；u(x)是待求未知量；当b(x)=0时，模型方程(1)为对流扩散方程．将求解区间[X1,X2]等分为N个子区间在点xi处由泰勒级数展开得到其中2.1 常系数对流扩散反应方程的指数型高精度紧致差分格式考虑常系数对流扩散反应方程将方程(4)改写为对于常系数对流扩散方程文献[18]中给出了形如的O(h4)阶精度的指数型紧致差分格式，其中系数为下面推导常系数对流扩散反应方程(4)的四阶指数型紧致差分格式．对(5)式中的第一个方程，考虑其在点xi处的四阶指数型紧致有限差分格式由(5)式中的第二个方程直接求导得将(10)式代入(9)中，整理得将(2)式代入(11)式右端项，略去高阶项整理得常系数对流扩散反应方程(4)的四阶指数型高精度紧致差分格式其中2.2 变系数对流扩散反应方程的指数型高精度紧致差分格式考虑变系数对流扩散反应方程(1)在点xi处的形如(12)式的差分格式其中差分系数由(13)给出．由泰勒级数展开得差分格式(14)的修正方程为由(15)式知，差分格式(14)是方程(1)的O(h2)阶近似差分，要得到其O(h4)的差分格式，需要在原模型方程(1)的左端添加进行修正，从而得到模型方程(1)的修正微分方程其中然后，对修正微分方程(16)考虑在点xi处的形如(12)式的差分格式，有其中差分格式(18)和(19)即为对流扩散反应方程(1)的四阶指数型紧致差分格式．(17)和(18)式中的fx,fxx可以直接求导得到，或者采用二阶中心差分算子(3)离散得到，都能保证差分格式(18)的O(h4)精度．将差分算子(3)代入到(18)式中，整理得(20)式所对应的代数方程组可以采用追赶法直接求解．3 差分格式的性质当c̸=0时，记则由于当时，故当雷诺数很大时，在粗网格差分格式(7)收敛阶也能达到三阶．若要差分格式(7)达到四阶精度，则需要满足因此，对于差分格式(7)有如下结论：定理1 1) 如果ac>0，则当网格步长满足时，差分格式(7)具有四阶精度；2) 如果ac<0：(i) 当c>−16a时，则当对任意的网格步长h，差分格式(7)具有四阶精度；(ii) 当c<−16a时，则当网格步长满足时，差分格式(7)具有四阶精度．证明不妨假设a>0，则分c>0和c<0两种情况证明．由(22)式得1) 当c>0时，有2ch2+ch−2a<0，由于∆=c2+16ac>0，故2) 当c<0时，有−2ch2+ch+2a>0，由于∆ =c2+16ac>0，故：当c>−16a时，则对任意的网格步长h，(22)式均成立；当c<−16a时，得若要差分格式(12)能达到四阶精度，则除(22)式外还需要满足因此，对于差分格式(12)，我们有如下结论：定理2 假设a>0，那么：1) 如果bc>0，则：(i) 如果b>0,c>0，当网格步长满足时，差分格式(12)具有四阶精度；(ii) 如果b>0,c<0，则当网格步长满足时，差分格式(12)具有四阶精度；2) 如果bc<0，则当网格步长满足时，差分格式(12)具有四阶精度．证明假设a>0，由(24)式得1) 如果bc>0，则：(i) 当b>0,c>0时，有即4ch2+(b+1)ch−2a(b+1)<0，由于∆=(b+1)2c2+32a(b+1)c>0，故(ii) 当b<0,c<0时，有即(b−1)ch+2a(b−1)<0，从而得2) 如果bc<0，则：(i) 当b<0,c>0时，有即(1−b)ch+2a(b−1)<0，从而得(ii) 当b>0,c<0时，有即由于函数是(0,+∞)上关于b的单调递减函数，故有从而得综上(i)和(ii)所述，如果bc<0，则当网格步长满足时，差分格式(12)具有四阶精度．4 数值算例下面，我们将选取典型算例采用本文格式(12)(简称为EHOC)进行计算，并与精确解和文献[17]中的多项式型格式(简称为FOC)的计算结果进行比较，验证本文格式的精确性和可靠性．收敛阶通过公式计算得到，其中err(N1)和err(N2)分别为网格数为N1和N2时最大绝对误差．算例1 −εuxx+ux= επ2sin(πx)+ πcos(πx)精确解为当ε很小时，该精确解在x=1处有一边界层．此算例中，a=ε,c=1，则当时，计算收敛阶为四阶．表1比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．计算数值结果表明，ε=0.1时，FOC格式的计算误差和EHOC格式的计算误差达到同一个数量级，但当ε=0.01,0.001时，EHOC格式的计算精度明显优于FOC格式，即随着ε取值减小，EHOC格式较FOC格式在计算精度方面的优势明显增加．图1和图2比较了当ε=10−3,10−5时采用EHOC格式和FOC格式的计算结果．从图中明显的看到，在粗网格上(N=16)，EHOC格式的计算解和精确解在计算节点上吻合的很好．表1: 算例1最大绝对误差比较图1: 当ε=10−3，网格数=16,128时，采用EHOC、FOC格式的计算解和精确解的比较图2: 当ε=10−5，网格数=16,128,1024,8192时，采用EHOC、FOC格式的计算解和精确解的比较算例2 −εuxx+ux+(1+ε)u=0,0<x<1精确解为此算例中，a=ε,c=1,b=1+ε，满足bc>0且b>0,c>0，则当时，计算收敛阶为四阶．表2比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．从表格中数值结果我们可以看到，对任何的ε取值，EHOC格式的计算精度都高于FOC格式的计算精度．表2: 算例2最大绝对误差比较算例3 −εuxx−ux+εu=sinx精确解为其中此算例中，a=ε,c=−1,b=ε，满足bc<0，则当时，计算收敛阶为四阶．表3比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．从表格中数值结果我们可以看到，对任何的ε取值，EHOC格式的计算精度都高于FOC格式的计算精度．图3和图4给出了当ε=10−3,10−5时，采用EHOC格式和FOC格式的计算数值结果和精确解的比较．图3中，当ε=10−3时，EHOC格式在粗网格(N=16)上能够很好的逼近精确解(此时最大绝对误差为2.94×10−4)，而FOC格式若要获得相当量级的计算精度，则至少需要4200个网格点(此时最大绝对误差为3.18×10−4)．而当ε=10−5时，如图4(d)所示，当网格数增加到9000时，在靠近边界层的地方，FOC格式的计算解和精确解的误差(为1.42)很大，而此时EHOC格式的计算误差为5.76×10−10．表3: 算例3最大绝对误差比较εFOC[17]网格数收敛阶误差误差收敛阶89.58(−2)2.65(−3)16 0.11.22(−2)2.983.54(−4)2.90328.72(−4)3.802.58(−5)3.89 645.22(−5)4.061.56(−6)4.08 1282.91(−2)8.14(−6)0.012562.43(−3)3.587.15(−7)3.51 5121.43(−4)EHOC4.094.25(−8)4.07 10248.93(−6)4.002.68(−9)3.99 5122.92(−1)7.19(−7)0.00110245.78(−2)2.341.57(−7)2.1920485.81(−3)3.311.69(−8)3.22 40963.54(−4)4.041.05(−9)4.00图3: 当ε=10−3，网格数=16,128,4200时，采用EHOC、FOC格式的计算解和精确解的比较算例4精确解为当ε<<1时，该解在x=1处有一边界层．此算例中，满足bc>0且b>0,c>0，则当时，计算收敛阶为四阶．图4: 当ε=10−5，网格数=16,128,9000时，采用EHOC、FOC格式的计算解和精确解的比较表4比较了当ε取不同值时，采用EHOC格式和FOC格式计算的最大绝对误差．计算数值结果表明，EHOC格式较FOC格式在计算精度方面具有明显的优势．表4: 算例4最大绝对误差比较图5和图6给出了当ε=10−3,10−5时，采用EHOC格式和FOC格式的计算数值结果和精确解的比较．从图中明显的看到，在强对流和反应占优的情况下，即使在粗网格上(N=16)，EHOC格式都能很好的逼近精确解．当ε=10−3时，若要获得相当的计算精度，采用FOC格式计算需要在620个网格下的计算，计算量是EHOC格式计算的38.75倍(EHOC的最大误差为4.33×10−4，FOC的最大误差为4.48×10−4)．当ε=10−5时，如图6(c)所示，将网格点增加255倍，采用FOC格式计算在光滑区域能很好的和精确解吻合，但在边界处误差很大，而EHOC格式能够很好的逼近精确解：EHOC的最大误差为4.25×10−8，FOC的最大误差为7.45×10−1．图5: 当ε=10−3，网格数=16,620时，采用EHOC、FOC格式的计算解和精确解的比较图6: 当ε=10−5，网格数=16,4096时，采用EHOC、FOC格式的计算解和精确解的比较5 结论本文基于常系数对流扩散方程的四阶指数型高精度紧致差分格式，结合残量修正法发展了一种数值求解一维变系数对流扩散反应方程的指数型高精度紧致差分方法，并选取了四个典型算例与文献[18]中的四阶多项式型格式的计算结果进行了比较，所有算例计算结果与表5中的理论分析结论一致，对于对流(反应)占优问题、边界层问题和大梯度变化的物理量问题本文格式的计算精度要明显优于文献[17]的格式．表5: 当ε取不同值时，差分格式达到四阶精度的条件参考文献：[1]Radhakrishna Pillai A C.Fourth-order exponential finite difference methods for boundary value problems of convective diffusiontype[J].International Journal for Numerical Methods inFluids,2001,37(1):87-106[2]Chen G Q,Gao Z.A perturbational exponential finite difference schemefor the convection diffusion equation[J].Journal of Computational Physics,1993,104(1):128-138[3]Phongthanapanich S,Dechaumphai bined finite volume and finite element method for convectiondiffusion-reaction equation[J].Journal of Mechanical Science and Technology,2009,23(3):790-801[4]Bause M,Schwegler K.Higher order finite element approximation of systems of convection-diffusion-reaction equations with smalldiffusion[J].Journal of Computational and AppliedMathematics,2013,246:52-64[5]Ayuso B,Marini L D.Discontinuous Galerkin methods for advection-diffusion-reaction problems[J].SIAM Journal on NumericalAnalysis,2009,47(2):1391-1420[6]Shih Y,Kellogg R B,Tsai P.A tailored finite point method for convection-diffusion-reaction problems[J].Journal of ScientificComputing,2010,43(2):239-260[7]Angelini O,Brenner K,Hihorst D.A finite volume method on general meshes for a degenerate parabolic convection-reaction-diffusion equation[J].Numerische Mathematik,2013,123(2):219-257[8]张建松，朱江，郭会，等.对流扩散反应方程的特征分裂最小二乘方法[J].高等学校计算数学学报，2012,34(4):289-299 Zhang J S,Zhu J,Guo H,et al.A characteristic splitting least-squares method for convection-diffusionreaction equations[J].Numerical Mathematics:A Journal of Chinese Universities,2012,34(4):289-299[9]朱国庆，李清善，陈绍春.分层网格上对流扩散反应方程的双二次元逼近[J].河南大学学报(自然科学版)，2007,37(3):221-225 Zhu G Q,Li Q S,Chen SC.Biquadraic element approximation of convection-diffusion-reactionequation under graded meshes[J].Journal of Henan University(Natural Science Edition),2007,37(3):221-225[10]崔翔鹏，贺力平.非线性对流扩散反应方程的预估-校正单调迭代差分方法[J].上海交通大学学报，2007,41(10):1731-1736 Cui X P,He L P.Predictor-corrector monotone iterative difference method for nonlinear convection diffusion reaction equation[J].Journal of Shanghai JiaotongUniversity,2007,41(10):1731-1736[11]朱海涛，欧阳洁.对流-扩散-反应方程的变分多尺度解法[J].工程数学学报，2009,26(6):997-1004 Zhu H T,Ouyang J.Variational multiscale method for the advection-diffusion-reaction equation[J].Chinese Journal of Engineering Mathematics,2009,26(6):997-1004[12]魏剑英.定常对流扩散反应方程的指数型高阶紧致差分格式[J].宁夏大学学报(自然科学版)，2012,33(2):140-143 Wei J Y.High-order exponential finite difference method for 1D convection-diffusion-reaction equation[J].Journal of Ningxia University(Natural Science Edition),2012,33(2):140-143[13]杨录峰，李春光.一种求解对流扩散反应方程的高阶紧致差分格式[J].宁夏大学学报(自然科学版)，2013,34(2):101-104 Yang L F,Li C G.A high-order compact finite difference scheme for solving the convection difference reaction equations[J].Journal of Ningxia University(Natural Science Edition),2013,34(2):101-104[14]田芳，田振夫.定常对流扩散反应方程非均匀网格上高精度紧致差分格式[J].工程数学学报，2009,26(2):219-225 Tian F,Tian Z F.A high accuracy compact difference scheme for convection diffusion reaction equation on non-uniform grid[J].Chinese Journal of EngineeringMathematics,2009,26(2):219-225[15]兰斌，薛文强，葛永斌.对流扩散反应方程基于坐标变换的高阶紧致差分格式[J].青岛科技大学学报(自然科学版)，2014,35(1):100-106 Lan B,Xue W Q,Ge Y B.A High-order compact difference scheme based on the coordinate transformation for the convection diffusion reaction equation[J].Journal of Qingdao University of Science and Technology(Natural Science Edition),2014,35(1):100-106[16]Spotz W F.High-order compact finite difference schemes for computational mechanics[D].Austin:University of Texas at Austin,1995 [17]Sun H W,Zhang J.A high order finite difference discretization strategy based on extrapolation for convection diffusion equations[J].Numerical Methods for Partial Differential Equations,2004,20(1):18-32[18]Tian Z F,Dai S Q.High-order compact exponential finite difference methods for convection-diffusion type problems[J].Journal of Computational Physics,2007,220(2):952-974。
