求解一维磁流体方程的无震荡中心差分格式

基于一维对流方程差分格式的哲学思考【摘要】水力学中对流扩散方程都具有通用微分方程的形式，如N-S方程组中的运动方程等。

通用微分方程通常用数值方法进行求解。

对流项的离散格式对差分方程的稳定性有着重要影响。

一维对流方程是研究对流项特性的模化方程。

对于初值问题，对流项差分格式为一阶精度时，差分方程是有条件稳定的；而对流项差分格式为二阶精度时，差分方程却是不稳定的。

从物理实质上，后者与实际情况也是不符的。

由此可见，局部高精度并不能保证整体的高精度，甚至会导致整体稳定性的破坏。

这表明，做事都要从实际出发，个人利益往往要服从集体利益。

【关键词】哲学；水力学；对流方程；差分格式0 引言同类事物的共同本质及其所表现出来的共同现象，即共性。

由于不同学科的复杂程度或发展的快慢不同，研究成果所达到的水平也不同。

复杂程度是学科研究内容本身所固有的，学科发展的快慢受社会需求的影响。

[1-2]人文科学相对于自然科学是复杂的。

就此而言，人文科学相对于自然科学研究水平相对较低，是必然的也是客观的。

所以，基于同类事物的共性，人文科学可以借鉴自然科学的研究成果。

水力学是自然科学的一个分支。

数值方法是目前水力学常用的研究手段之一。

基于一维对流方程差分格式等进行一些哲学思考是有益的。

1 通用微分方程及一维对流方程[3]1.1 通用微分方程通用微分方程的形式可写为：■+■=■Γ■■+S■（1）上式为张量式。

式中Φ为通用变量，Γ■是扩散系数，S■为源项。

方程各项依次为时变项、对流项、扩散项和源项。

对应于特定意义的Φ，Γ■和S■具有特定的形式。

该通用微分方程可以转化为连续方程、运动方程、雷诺方程、浓度输运方程或热输运方程等。

1.2 一维对流方程运动方程是水力学通用微分方程中最常见的形式。

张量式为：■+■=■ρv■-■+ρfi（2）对于不可压缩流体ρ为常数，并且忽略扩散项和源项（含压强梯度项和质量力项），则上式可简化为：■+■=0（3）一维问题的形式为：■+u■=0（4）对流项线性化之后的形式为：■+α■=0（5）式中α是不为零的常数。

应用数学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。

有限差分求解1维流体方程
有限差分法是一种数值求解偏微分方程的方法，用于离散化空间和时间上的导数，将偏微分方程转化为代数方程组。

对于一维流体方程，我们可以采用有限差分法来求解。

首先，我们需要考虑一维流体方程的形式。

一维流体方程通常包括质量守恒方程和动量守恒方程。

质量守恒方程描述了流体的质量随时间和空间的变化，而动量守恒方程描述了流体的运动状态随时间和空间的变化。

在有限差分法中，我们需要将空间和时间上的导数用差分形式表示。

对于一维空间上的导数，我们可以使用中心差分、向前差分或向后差分等方法进行离散化。

对于时间上的导数，通常使用向前差分或向后差分。

一旦我们将偏微分方程离散化为代数方程组，我们可以使用数值方法（如迭代法、矩阵求解法等）来求解这个代数方程组。

在求解过程中，需要考虑数值稳定性和收敛性等问题。

另外，对于一维流体方程，我们还需要考虑边界条件和初始条
件的处理。

边界条件和初始条件的选择对数值求解的精度和稳定性有很大影响，因此需要仔细考虑和处理。

总之，有限差分法是一种常用的数值求解偏微分方程的方法，对于一维流体方程，我们可以通过离散化空间和时间上的导数，将偏微分方程转化为代数方程组，然后使用数值方法进行求解。

在实际应用中，需要考虑边界条件、初始条件以及数值稳定性等因素。

一、上机目的用数值方法计算一维对流方程在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.周文龙的贡献
正文
1.流体力学简介
流体力学是研究流体在静止和运动状态下的行为和性质的学科，其中流体可以是液体或气体。

流体力学在工程、科学和日常生活中具有广泛的应用，如管道设计、飞机翼的构造、心脏的运作等。

在计算机图形学和数值计算领域，流体力学的应用也日益广泛。

2.差分格式的概念与应用
差分格式是一种数学方法，用于离散化流体力学中的偏微分方程。

通过差分格式，可以将连续的数学模型转化为离散的计算问题，从而便于计算机进行求解。

差分格式在流体力学中的应用包括计算流体动力学（CFD）和有限元方法等。

3.周文龙的贡献
周文龙是一位在计算流体力学领域具有重要影响的学者。

他对差分格式进行了深入研究，提出了一种高效的流体力学差分格式。

周文龙的差分格式在求解流体力学问题时具有较高的精度和稳定性，被广泛应用于各种实际问题中。

第1页共1页。

二维磁流体方程组的差分方法李富智;于佳利;杨慧【摘要】对二维磁流体方程组的初边值问题建立了有限差分格式,并分析了差分格式的截断误差.最后,通过Matlab数值模拟了差分结果,并与精确解进行了误差比较,验证了理论结果.【期刊名称】《云南师范大学学报（自然科学版）》【年(卷),期】2016(036)002【总页数】5页(P22-26)【关键词】磁流体方程;有限差分格式;截断误差;数值模拟【作者】李富智;于佳利;杨慧【作者单位】云南师范大学数学学院,云南昆明650500;云南师范大学数学学院,云南昆明650500;云南师范大学数学学院,云南昆明650500【正文语种】中文【中图分类】O242.1;O246磁流体是一种良好的功能材料，既具有液体的流动性，受到磁场作用后又具有固体材料的磁性.鉴于磁流体的优良性质，其在工业生产、天体物理及核技术等领域内的应用日益广泛[1].目前，国内外研究人员仍然从实验研究、理论分析和数值模拟三个方面对磁流体进行研究，研究结果基本一致.理论分析不仅有助于清晰地认识导电流体的相关性质，而且能使得实验研究和数值模拟得以顺利开展[2].本文主要研究求解二维磁流体力学方程组近似解的数值方法.磁流体力学方程组的数值解法在空间物理、热交换、空间推进和流动控制等众多领域具有广泛的应用[3-4].考察二维不可压缩的磁流体方程组其中磁场B通过一个磁流函数ψ来定义B=▽⊥ψ，且j=Δψ.速度u通过一个速度流函数φ来定义u=▽⊥φ[5].并且常系数μ≥0.本文只考虑μ>0的情形.假设Ω⊂R2是一个有界区域，现对Ω进行网格剖析.记时间步长为τ，即Δt=τ，空间步长为h，为了方便起见使x方向与y方向上的步长相等，即Δx=Δy=h.这样网格线就可以写作用有限差分方法近似求解偏微分方程问题有多种多样的方法，并且也可以用不同的构造方法来建立这些有限差分方法.下面就用Taylor级数展开法来建立差分格式[6]. 首先将方程组(1)转化成如下的形式：假定二维不可压缩的磁流体方程组问题的解φ(x,y,t)是充分光滑的.由Taylor级数展开有同理有ψ(x,y,t)、w(x,y,t)的Taylor级数展开.其中表示括号内的函数在节点(xj,yk,tn)处的取值.现假设为φ(xj,yk,tn)的近似值，用近似值来代替原方程，这样就得到原方程组的有限差分方程，变形后形式如下对方程(1)的解，关于t的向前差分的Taylor级数展开有Δ+tw(x,y,t)=w(x,y,t+τ)-w(x,y,t)对变量x进行Taylor级数展开有同理可以得到ψ(x,y,t)、φ(x,y,t)的Taylor级数展开.故得截断误差为所以有下面考察二维不可压缩的磁流体方程组现在给出一组精确解来进行模拟，为方便起见，给出了一组特殊的精确解得到精确解与数值解模拟如下.综上可以发现，当t=4-6时，差分格式吻合度较好，随着时间的推移，误差图的吻合度较好，变化程度也趋于稳定.以上得到的图像都是采用时间步长为0.05时的情况.我们还发现随着时间步长越来越小，差分解的精确性也越来越高.由此可以说明本文所采用的差分格式是可行的.现在考虑该差分格式的收敛性，将方程组中的三个变量φ、ψ、w分别进行模拟来考察.针对(9)通过Matlab软件，得到图4、图5.图4、图5为在时间与空间步长越来越细的时候，中心差分数值解每相邻5层的差值标准差图像.分析图像可知，随着时间的推移和时间层数的增加，数值解间的距离并不趋于0，所以磁流体方程组(1)的差分格式一般不是收敛的.【相关文献】[1] 刘德明.用差分法解一维磁流体力学方程组[J].数值计算与计算机应用， 1984 (3)：138-150.[2] 陶詹晶，郑华盛.二维磁流体方程的非交错无振荡中心差分格式[J].南昌航空大学学报：自然科学版，2011，25(1)：78-84.[3] 陆金甫，关治. 偏微分方程数值解法[M].2版.北京：清华大学出版社，2003.[4] 张德喜，赵磊生. MATLAB语言程序设计教程[M].2版.北京：中国铁道出版社,2010.[5]CORDOBA DIEGO,MARLIANI CHRISTIANE.Evolution of current sheets and regularity of ide al incompressible magnetic fluids in 2D[J].Commun on pure and Applied Mathematics，2000,53(4):512-524.[6] 郑尚昆，王彩华.非等距网格上奇异扰动问题的有限差分格式[J].天津师范大学学报：自然科学版，2014，34(3)：11-16.。

一维热传导方程的差分格式一维热传导方程的差分格式《微分方程数值解》学生姓名1：许慧卿学号： 20214329 学生姓名2：向裕学号： 20214327 学生姓名3：邱文林学号： 20214349 学生姓名4：高俊学号： 20214305 学生姓名5：赵禹恒学号： 20214359 学生姓名6：刘志刚学号： 20214346 学院：理学院专业： 14级信息与计算科学指导教师：陈红斌2021年6 月25日《一维热传导方程的差分格式》论文一、《微分方程数值解》课程论文的格式 1）引言：介绍研究问题的意义和现状 2）格式：给出数值格式3）截断误差：给出数值格式的截断误差 4）数值例子：按所给数值格式给出数值例子 5）参考文献：论文所涉及的文献和教材二、《微分方程数值解》课程论文的评分标准 1）文献综述：10分；2）课题研究方案可行性：10分； 3）数值格式：20分；4）数值格式的算法、流程图：10分； 5）数值格式的程序：10分；6）论文撰写的条理性和完整性：10分； 7）论文工作量的大小及课题的难度：10分； 8）课程设计态度：10分； 9）独立性和创新性：10分。

一维热传导方程的差分格式考虑如下一维非齐次热传导方程Dirichlet 初边值问题=a 2+f (x , t ), cu (x ,0) =ϕ(x ), c ≤x ≤d , (1.2)u (c , t ) =α(t ), u (d , t ) =β(t ), 0的有限差分方法, 其中a 为正常数, f (x , t ), ϕ(x ), α(t ),β(t ) 为已知常数, ϕ(c ) =α(0),ϕ(d ) =β(0). 称(1.2)为初值条件, (1.3)为边值条件.本文将给出(1.1) (1.3)的向前Euler 格式, 向后Euler 格式和Crank -Nicolson 格式, 并给出其截断误差和数值例子. 经对比发现, Crank -Nicolson 格式误差最小, 向前Euler 格式次之, 向后Euler 格式误差最大.2 差分格式的建立2.1 向前Euler 格式将区间[c , d ]作M 等分, 将[0, T ]作N 等分, 并记 h =(d -c ) /M , τ=T /N ,x j =c +jh , 0≤j ≤M , t k =k τ, 0≤k ≤N . 分别称h 和τ为空间步长和时间步长. 用x =x j , 0≤j ≤M , t =t k , 0≤k ≤N将Ω分割成矩形网格. 记Ωh =x j |0≤j ≤M , Ωτ={t k |0≤k ≤N }, Ωhτ=Ωh ⨯Ωτ. 称x j , t k 为结点.定义Ωh τ上的网格函数Ω=U j |0≤j ≤M ,0≤k ≤N , 其中U j =u x j , t k .在结点x j , t k 处考虑方程(1.1),有∂u (x j , t k )=a∂2u (x j , t k )+f (x j , t k ), 1≤j ≤M -1, 1≤k ≤N -1. (2.1)将u (x j , t k +1)以结点(x j , t k )为中心关于t 运用泰勒级数展开, 有u (x u ''(x j , t k )τ2j , t k +1)=u (x j , t k )+u '(x j , t k )τ++o (τ2).u (x 2j , t k +1)-u (x j , t k )=∂u (x j , t k )τ∂u (x j , t k )τ∂t +2∂t+o (τ). (2.2) 再将(x j -1, t k ), (x j +1, t k )分别以结点(x j , t k )为中心关于x 运用泰勒级数展开, 有u (x j , t k （)-h ）j -1, t k )=u (x j , t k )+u '(x j , t k )(-h )j , t k )(-h ）+u ''(x 2!u '''(x 3u (4)(x j , t k )(-h )+o (h 4) , (x u (x u ''(x j , t k )h 2u '''(x j , t k )h 3u j +1, t k )=j , t k )+u '(x j , t k )h +u (4)(x j , t k )h 4+o (h 4).由上述两式可得u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )∂2u (x j , t k h 2=)∂x 2+h 2∂4(x j , t k )12∂x+o (h 2) . (2.3) 将(2.2), (2.3)两式代入(2.1)中, 得u (x j , t k +1)-u (x j , t k )k )+u (x j +1, t k )u (x j -1, t k )-2u (x j , t h+f (x j , t k )+R k j . (2.4)2其中R k ah ∂4u (x j , t k)τ2∂2(x j , t k )j=-12∂x 4+2∂t 2+o (τ+h 2) 为方程(2.1)的截断误差. 舍去截断误差, 用u kj 代替u (x j , t k ), 得到如下差分方程u k +1-u k j ju k k kj -1-2u j +u j +1+f k j , 1≤j ≤M -1, 0≤k ≤N -1. 结合初边值条件, 可得如下差分格式+1u k -u k j ju k j -1-2u j +u j +1+f j k , 1≤j ≤M -1, 1≤k ≤N -1, (2.6)u 0j =ϕ(x j ), 0≤j ≤M , (2.7)k k u 0=α(t j ), u M =β(t k ), 1≤k ≤N . (2.8)2.2 向后Euler 差分格式在结点x j , t k +1处考虑方程(1.1), 有∂u (x j , t k +1)∂t∂2u (x j , t k +1)∂x 2+f (x j , t k +1), 1≤j ≤M -1, 1≤k ≤N . (2.9)将u x j , t k 以x j , t k +1为中心关于t 运用泰勒级数展开, 有u (x j , t k )=u (x j , t k +1)+u '(x j , t k +1)(-τ) +u ''(x j , t k +1)(-τ) 2+o (τ2).u (x j , t k +1)-u (x j , t k )∂u (x j , t k +1)τ∂2u (x j , t k +1)=-+o (τ). (2.10) 2再将u x j -1, t k +1, u x j +1, t k +1分别以x j , t k +1为中心关于x 运用泰勒级数展开, 有u (x j -1, t k +1)=u (x j , t k +1)+u '(x j , t k +1)(-h ）+u ''(x j , t k +1)(-h )u '''(x j , t k +1（)-h ）u (4)(x j , t k +1)(-h )+o (h 4), u ''(x j , t k +1)h 22! +o (h 4).u '''(x j , t k +1)h 3u (x j +1, t k +1)=u (x j , t k +1)+u '(x j , t k +1)h +u (4)(x j , t k +1)h 4由上述两式可得u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)∂2u (x j , t k +1)h 2∂4(x j , t k +1)=++o (h 2) . (2.11) 224h ∂x 12∂x将(2.10), (2.11)两式代入(2.9)中, 得+f (x j , t k +1)+R k j . (2.12)τ∂(x i , t k +1)ah 2∂u (x i , t k +1)+o (τ+h 2) 为方程(2.9)的截断误差.舍去截断误差, 用u j 代替u x j , t k , 得到如下差分方程+1u k -u k j j+1k +1+1u k +u k j -1-2u j j +1+f j k +1, 1≤j ≤M -1, 1≤k ≤N . (2.13)结合初边值条件, 可得如下差分格式u k -u k j j+1k +1+1u k +u k j -1-2u j j +1+f j k +1, 1≤j ≤M -1, 1≤k ≤N , (2.14)u 0j =ϕ(x j ), 0≤j ≤M , (2.15)k k u 0=α(t k ), u M =β(t k ), 1≤k ≤N . (2.16)2.3 Crank -Nicolson 差分格式在结点x j , t k +1/2处考虑方程(1.1), 有∂u (x j , t k +1/2)∂t∂2u (x j , t k +1/2)∂x+f (x j , t k +1/2), 1≤j ≤M -1, 0≤k ≤N -1. (2.17)将u x j , t k +1, u x j , t k 以x j , t k +1/2为中心关于t 运用泰勒级数展开, 有u ''(x j , t k +1/2)⎛τ⎛2u (x j , t k +1)=u (x j , t k +1/2)+u '(x j , t k +1/2)+ ⎛22! ⎛2⎛u '''(x j , t k +1/2)⎛τ⎛3+ ⎛+o (τ).τ⎛u ''(x j , t k +1/2)⎛τ⎛⎛u (x j , t k )=u (x j , t k +1/2)+u '(x j , t k +1/2) -⎛+ -⎛22! ⎛⎛⎛2⎛u '''(x j , t k +1/2)⎛τ⎛3+ -⎛+o (τ).将上述两式整理得=++o (τ3). (2.18) 3∂t 24∂t再将u x j -1, t k , u x j +1, t k 分别以x j , t k 为中心关于x 运用泰勒级数展开, 有u (x j -1, t k )=u (x j , t k )+u '(x j , t k )(-h ）+u ''(x j , t k )(-h )u '''(x j , t k （) -h ）u (4)(x j , t k )(-h )4! u ''(x j , t k )h 2+o (h 4) , u '''(x j , t k )h 3u (4)(x j , t k )h 4u (x j +1, t k )=u (x j , t k )+u '(x j , t k )h ++o (h 4).u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )∂2u (x j , t k )h 2∂4(xj , t k )2=++o (h ). (2.19) 224h ∂x 12∂x同理, 将u x j -1, t k +1, u x j +1, t k +1分别以x j , t k +1为中心关于x 泰勒级数展开, 整理得u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)∂2u (x j , t k +1)h 2∂4(x j , t k +1)=++o (h 2). (2.20) 224h ∂x 12∂x此时，分别将u x j , t k +1, u x j , t k 以u x j , t k +1/2为中心关于t 泰勒级数展开，有∂2u (x j , t k +1)∂2u (x j , t k +1/2)∂3u (x j , t k +1/2)τ1∂4u (x j , t k +1/2)⎛τ⎛2=+++o τ, () ⎛22222∂x ∂x ∂x ∂t 22! ∂x ∂t ⎛2⎛∂2u (x j , t k )∂2u (x j , t k +1/2)∂3u (x j , t k +1/2)⎛τ⎛=+ -⎛ 222∂x ∂x ∂x ∂t ⎛2⎛1∂u (x j , t k +1/2)⎛τ⎛2+ -⎛+o (τ). 222! ∂x ∂t ⎛2⎛利用上述两式得∂2u (x j , t k +1/2)∂x 21⎛∂u (x j , t k +1)∂u (x j , t k )⎛1⎛∂u (x j , t k +1/2)⎛τ⎛⎛2⎛ ⎛= +-+o τ. () ⎛2222⎛ ⎛2∂x ∂x 2∂x ∂t ⎛2⎛⎛⎛⎛⎛利用(2.19), (2.20)两式, 整理有∂2u (x j , t k +1)∂x∂2u (x j , t k )∂xu (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)h 2∂(x j , t k ) -h 2∂(x j , t k +1)2-+o h . (2.22) ()412∂x结合(2.21),(2.22)两式, 整理得∂2u (x j , t k +1/2)∂xu (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)1⎛∂u (x j , t k +1/2)⎛τ⎛⎛- ⎛⎛⎛2 ∂x 2∂t 22⎛⎛⎛⎛⎛1⎛h 2∂(x j , t k )h 2∂(x j , t k +1)2- ++o (h )⎛. (2.23) 44⎛212∂x 12∂x ⎛⎛将(2.18), (2.23)式代入(2.17)得u (x j , t k +1)-u (x j , t k )u (x j -1, t k )-2u (x j , t k )+u (x j +1, t k )u (x j -1, t k +1)-2u (x j , t k +1)+u (x j +1, t k +1)+f (x j , t k )+R k j . (2.24)⎛h 2∂4(x j , t k )h 2∂4(x j , t k +1)∂4u (x j , t k +1/2)⎛τ⎛2⎛τ2∂3u (x j , t k +1/2)a⎛+R k ++ j =- ⎛44223⎛2 12∂x 12∂x ∂x ∂t 224∂t ⎛⎛⎛⎛+o (τ3+h 2) 为方程(2.17)的截断误差.舍去截断误差, 用u j 代替u (x j , t k ) , 可得如下差分方程u k -u k j ju k j -1-2u j +u j +1+1k +1+1u k +u k j -1-2u j j +1+f j k +1/2,1≤j ≤M -1, 0≤k ≤N -1. (2.24)结合初边值条件, 可得如下差分格式u k -u k j ju k j -1-2u j +u j +1+1k +1+1u k +u k j -1-2u j j +1+f j k +1/2,1≤j ≤M -1, 0≤k ≤N -1, (2.25)u 0j =ϕ(x j ), 0≤j ≤M , (2.26)k k u 0=α(t k ), u M =β(t k ), 1≤k ≤N . (2.27)3 差分格式的求解3.1 向前Euler 格式记r =a τ, 称r 为步长比. 差分格式(2.6) (2.8)中(2.6)可改写为u k +1j =r ⋅u k ) u k k kj -1+(1-2r j +r ⋅u j +1+τ⋅f j ,1≤j ≤M -1 , 0≤k ≤N -1 .将(3.1)写成如下形式U k +1=A ⋅U k +τ⋅F k +r ⋅F 1 .U k +1=⎛u k +1, u k +1, , u k +1⎛TU k =⎛u k , u k k⎛12M -1⎛, ⎛12, , u M -1⎛, F k =⎛f k , f k , , f k⎛T F =⎛u k , 0, k ⎛T⎛12M -1⎛, 1⎛00, , 0, u M ⎛(M -1)*1, ⎛ 1-2r r ⎛r 1-2r r ⎛A = ⎛⎛r 1-2r ⎛⎛(M -1)*(M -1)3.2 向后Euler 格式记r =a τ, 称r 为步长比. 差分格式(2.14) (2.16)中(2.14)可改写为-u k j =r ⋅u k +1k +1k +1k +1j -1-(1+2r ) u j +r ⋅u j +1+τ⋅f j1≤j ≤M -1 , 1≤k ≤N .将(3.2)写成如下形式A ⋅U k +1=-U k -τ⋅F k +1-r ⋅F 1 .k +1k +1k +1k k k k⎛⎛⎛U k +1=⎛u , u , , u U =u , u , , u , 12M -112M -1⎛, ⎛⎛⎛ k +1k +1k +1k +1k +1⎛⎛⎛F k +1=⎛f , f , , f F =u , 0, 0, , 0, u , 12M -110M ⎛⎛⎛⎛T T (M -1)*1r ⎛-(1+2r ) ⎛r -(1+2r ) r ⎛⎛A = . rr -(1+2r ) ⎛⎛⎛(M -1)*(M -1)3.3 Crank -Nicolson 格式记r =a τ, 称r 为步长比. 差分格式(2.25) (2.27)中(2.25)可改写为+1k +1k +1k k k k +1/2, -[r /2⋅u k -(1+r ) u +r /2⋅u ]=r /2⋅u +(1-r ) u +r /2⋅u +τ⋅f j -1j j +1j -1j j +1j将(3.3)写成如下形式1≤j ≤M -1 , 0≤k ≤N -1 .. -A 2⋅U k +1=A 1⋅U k +τ⋅F k +1/2+r (F 1+F 2）k +1k +1k +1k k k k⎛⎛⎛U k +1=⎛u , u , , u U =u , u , , u , 12M -112M -1⎛, ⎛⎛⎛k k k +1k +1F 1=⎛⎛u 0, 0, 0, , 0, u M ⎛⎛(M -1)*1, F 2=⎛⎛u 0, 0, 0, , 0, u M ⎛⎛(M -1)*1,k +1/2k +1/2⎛=⎛, f 2k +1/2, , f M -1⎛, ⎛f 1⎛1-r r /2⎛⎛r /21-r r /2 ⎛⎛A 1= , r /2r /21-r ⎛⎛⎛(M -1)*(M -1)r /2⎛-(1+r ) ⎛ ⎛r /2-(1+r ) r /2 ⎛⎛A 2= . r /2r /2-(1+r ) ⎛⎛⎛(M -1)*(M -1)在本文中, 将应用向前Euler 格式(2.6) (2.8), 向后Euler 格式(2.14) (2.16)和Crank -Nicolson 格式(2.25) (2.27)分别来求解如下算例.算例现有如下定解问题∂u ∂2u x +2t=2+e , 0u (x ,0) =e x , 0≤x ≤1,u (0,t ) =e 2t , u (1,t ) =e 1+2t , 0上述定解问题的精确解为u (x , t ) =e 4.1 向前Euler 格式在表4.1.1中给出了取步长h =1/10和τ=1/200(步长比r =1/2) 时计算得到的部分数值结果. 表4.1.2给出了取步长h =1/10和τ=1/100(步长比r =1) 时计算得到的部分数值结果, 随着计算层数的增加, 误差越来越大, 数值结果是发散的. 经步长比r 的多次取值发现, 当r ≤1/2时, 其数值结果是收敛的, 当r >1/2时, 其数值结果是发散的. 表4.1.3给出了r =1/2时, 取不同步长, 数值解的最大误差E ∞(h , τ) =max |u (x j , t k ) -u j k |.1≤j ≤M -11≤k ≤N从表4.1.3可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/4时, 最大误差约缩小到原来的1/4.图4.1.1给出了t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线. 图4.1.2给出了t =1时不同步长数值解的误差曲线图(r =1/2) , 由图4.1.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.1.3给出了不同步长的误差曲面图.表4.1.1 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/200)(x , t )(0.5,0.1)数值解 1.802810精确解 1.803988|精确解-数值解|1.178e-360 80 100 120 140 160 180 200(0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)2.6886513.2838564.010886 4.8988975.983523 7.308290 8.926366 10.9026872.6912343.2870814.014850 4.9037495.989452 7.315534 8.935213 10.9134942.583e-33.225e-3 3.965e-34.852e-35.929e-3 7.243e-3 8.847e-3 1.081e-3表4.1.2 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/100)1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17(x , t )(0.5,0.01) (0.5,0.02) (0.5,0.03) (0.5,0.04) (0.5,0.05) (0.5,0.06)(0.5,0.07) (0.5,0.08) (0.5,0.09) (0.5,0.10) (0.5,0.11) (0.5,0.12) (0.5,0.13) (0.5,0.14) (0.5,0.15) (0.5,0.16) (0.5,0.17)数值解 1.521674 1.552123 1.583184 1.614870 1.647386 1.679785 1.716057 1.741446 1.807089 1,744143 2.092546 1.164923 4.148393 -4.709415 21.998093 -57.422399 178.294623精确解 1.521962 1.552707 1.584074 1.616074 1.648721 1.682028 1.7160071.750673 1.786038 1.822119 1.858928 1.896481 1.934792 1.9738782.0137532.054433 2.095936|精确解-数值解|2.879e-4 5.845e-4 8.901e-4 1.205e-3 1.336e-3 2.242e-3 5.051e-5 9.227e-3 2.105e-2 7.798e-2 2.336e-1 7.316e-1 2.213e+0 6.684e+0 1.998e+1 5.948e+11.762e+2表4.1.3 算例取不同步长时数值解的最大误差(r =1/2)1/10 1/20 1/40 1/801/200 1/800 1/3200 1/12800E ∞(h , τ)1.178e-22.973e-3 7.431e-4 1.858e-4E ∞(2h ,4τ) /E ∞(h , τ)* 3.964 4.000 4.000图4.1.1 算例取t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线图4.1.2 算例取t =1时不同步长数值解的误差曲线(r =1/2)图4.1.3 算例取不同步长数值解的误差曲面(r =1/2)4.2 向后Euler 格式在表4.2.1和表4.2.2中分别给出了h =1/10,τ=1/200(步长比r =1/2) 和h =1/10, τ=1/100(步长比r =1) 时计算得到的部分数值结果, 发现两表的数值结果都是收敛的, 经步长比r 的多次取值发现, 无论r 取任何值, 其数值结果都是收敛的. 表4.2.3 给出了r =1/2时, 取不同步长, 数值解的最大误差E ∞(h , τ) =max |u (x j , t k ) -u j k |.1≤j ≤M -11≤k ≤N从表4.2.3可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/4时, 最大误差缩小到原来的1/4. 表4.2.4给出了r =1时的类似结论.图4.2.1给出了t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线. 图4.2.2给出了t =1时, 取不同步长数值解的误差曲线图(r =1/2) , 由图4.2.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.2.3和图4.2.4分别当给出了步长比r =1/2和r =1下不同步长的误差曲面图.表4.2.1 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/200)20 40 60 80 100 120 140 160 180 200(x , t )(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)数值解 1.805343 2.205674 2.694258 3.290867 4.019509 4.909453 5.9964257.324052 10.926203 10.926203精确解 1.803988 2.203396 2.691234 3.287081 4.014850 4.903749 5.989452 7.315534 8.935213 10.913494|精确解-数值解|1.355e-32.278e-33.023e-3 3.785e-34.659e-35.704e-36.973e-3 8.518e-3 1.041e-2 1.271e-2表4.2.2 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/100)10 20 30 40 50 60 70 80 90 100(x , t )(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)数值解 1.788500 2.185742 2.670171 3.261551 3.983745 4.865788 5.943098 7.258921 8.866068 10.829041精确解 1.786038 2.181472 2.664456 3.254374 3.974902 4.854956 5.929856 7.242743 8.846306 10.804903|精确解-数值解|2.461e-3 4.269e-3 5.715e-3 7.177e-3 8.843e-3 1.083e-2 1.324e-2 1.618e-2 1.976e-2 2.414e-2表4.2.3 算例取不同步长时数值解的最大误差(r =1/2)1/10 1/20 1/40 1/801/200 1/800 1/3200 1/12800E ∞(h , τ)1.386e-2 3.509e-3 8.779e-42.195e-4E ∞(2h ,4τ) /E ∞(h , τ)* 3.950 3.997 3.999表4.2.4 算例取不同步长时数值解的最大误差(r =1)1/10 1/20 1/401/100 1/400 1/1600E ∞(h , τ)2.659e-2 6.744e-3 1.688e-3E ∞(2h ,4τ) /E ∞(h , τ)* 3.943 3.9551/6400 4.222e-4 3.999图4.2.1 算例取t =1时的数值解曲线(h =1/10,τ=1/200) 与精确解曲线图4.2.2 算例取t =1 时不同步长数值解的误差曲线(r =1/2)图 4.2.3 算例取不同步长数值解的误差曲面(r =图 4.2.4 算例取不同步长数值解的误差曲面(r =1)4.3 Crank -Nicolson 格式在表4.3.1中给出了取步长h =1/10,τ=1/10时计算得到的部分数值结果, 表4.3.2给出了取步长h =1/100,τ=1/100时计算得到的部分数值结果, 发现两表的数值结果都是收敛的, 经步长比r 的多次取值发现, 无论r 取任何值, 其数值结果都是收敛的. 表4.3.3给出了取不同步长时, 所得数值解的最大误差E ∞(h , τ) =max |u (x j , t k ) -u j k |.1≤j ≤M -11≤k ≤N从表可以看出当空间步长缩小到原来的1/2, 时间步长缩小到原来的1/2时, 最大误差缩小到原来的1/4.图4.3.1给出了t =1时的数值解曲线(h =1/10,τ=1/10)与精确解曲线. 图4.3.2给出了t =1时, 取不同步长数值解的误差曲线图, 由图4.3.2可知, 当步长比r 保持不变, 随着时间与空间方向上的加密, 误差越来越小. 图4.3.3给出了取不同步长时所得数值解的误差曲面图. 图4.3.4给出了向前Euler , 向后Euler 和Crank -Nicolson 三种方法数值解的误差对比图(h =1/10,τ=1/200) , 由图4.3.4分析可得, 当时间方向t 与空间方向x 都取相同步长时, Crank -Nicolson 方法误差最小, 精度最高, 向前Euler 方法次之,向后Euler 方法所得误差最大, 精度最低.表4.3.1 算例在部分结点处数值解、精确解和绝对误差(h =1/10,τ=1/10)1 2 3 4 5 6 7 8 9 10(x , t )(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)数值解 1.823114 2.227196 2.720414 3.322776 4.058460 4.957021 6.054520 7.395008 9.032284 11.032056精确解 1.822119 2.225541 2.718282 3.320217 4.055200 4.953032 6.049647 7.389056 9.025013 11.023176|精确解-数值解|9.949e-4 1.656e-3 2.132e-3 2.659e-3 3.260e-3 3.988e-3 4.873e-3 5.952e-3 7.271e-3 8.880e-3表4.3.2 算例在部分结点处数值解、精确解和绝对误差(h =1/100,τ=1/100)10 20 30 40 50 60 70 80 90 100(x , t )(0.5,0.1) (0.5,0.2) (0.5,0.3) (0.5,0.4) (0.5,0.5) (0.5,0.6) (0.5,0.7) (0.5,0.8) (0.5,0.9) (0.5,1.0)数值解 1.993726 2.435147 2.974297 3.632815 4.437131 5.919524 6.619421 8.084979 9.875016 12.061372精确解 1.993715 2.435130 2.974297 3.632786 4.437096 5.419481 6.619369 8.084915 9.874938 12.061276|精确解-数值解|1.081e-5 1.749e-52.296e-5 2.864e-53.521e-54.308e-55.266e-56.433e-57.857e-5 9.597e-5表4.3.3 算例取不同步长时数值解得最大误差1/10 1/20 1/40 1/80 1/160 1/320 1/6401/10 1/20 1/40 1/80 1/160 1/320 1/640E ∞(h , τ)9.588e-3 2.429e-3 6.078e-4 1.520e-4 3.799e-5 9.499e-6 2.375e-6E ∞(2h ,2τ) /E ∞(h , τ)* 3.9457 3.9961 3.9990 3.9998 3.9999 4.0000图4.3.1 算例取t =1时数值解曲线(h =1/10,τ=1/10)与精确解曲线图4.3.2 算例取t 1时不同步长数值解的误差曲线图4.3.3 算例取不同步长数值解的误差曲面图4.3.4 算例取t =1时三种方法数值解的误差对比曲线(h =1/10,τ=1/200)在本文中, 给出了向前Euler , 向后Euler 和Crank -Nicolson 三种差分格式, 分别对抛物型方程进行求解, 其中向前Euler 格式是条件稳定的, 后两者是无条件稳定的. 经对比发现, Crank -Nicolson 格式的精度最高, 误差最小, 向前Euler 格式次之, 向后Euler 格式精度最低, 误差最大.[1] 孙志忠. 偏微分方程数值解法[M]. 第二版. 北京: 科学出版社, 2021.[2] 李荣华. 偏微分方程数值解法[M]. 第二版. 北京: 高等教育出版社, 2021.[3] 刘卫国. MATLAB程序设计教程[M]. 第二版. 北京: 中国水利水电出版社, 2021.◆ 程序流程图1 向前Euler 法2 向后Euler 法3 Crank Nicolson 法(向前Euler 法）建立函数文件ForwardEuler.m ，如下：% 求解一维非齐次热传导方程 ut-uxx=e^(x+2t), 0% u(x,0)=e^x, 0% u(0,t)=e^(2t), u(1,t)=e^(1+2t), 0% 该问题的定解为 u(x,t)=e^(x+2t).% 其中取a=1, h1表示时间t 方向上的步长, h2表示空间x 方向上的步长.function [uxy,u,eps]=ForwardEuler(h1,h2,M,N,T)% 求解一维非齐次热传导方程 ut-uxx=f(x,t),(a>0)% 其中M,N 分别为x 与t 方向上的网格数。

对流方程的四阶中心差分格式为了求解流体力学问题，我们经常需要离散化偏微分方程以在计算机上进行数值求解。

流动问题通常可以通过连续性方程和动量方程来描述，在这些方程中，对流项是非常重要的。

在数值求解中，对流项的离散化是一个复杂的问题，因为它与流体流动的特性密切相关。

在本文中，我将讨论一种用于解决流体动力学问题的四阶中心差分格式，该格式对流项进行了有效的离散化。

首先，让我们考虑一维线性对流方程：\[\frac{\partial u}{\partial t} + c\frac{\partialu}{\partial x} = 0\]其中，\(u\) 是速度或浓度的函数，\(c\) 是对流速度。

为了离散化这个方程，我们需要在时间和空间上对其进行网格化。

假设我们使用均匀网格，网格间的空间步长为 \(\Delta x\)，时间步长为 \(\Delta t\)。

现在，我们可以使用四阶中心差分格式来离散化对流方程。

我们可以使用前三个时间步和三个空间步来定义一个差分方程，如下所示：\[\frac{u_i^{n+1} - u_i^{n-2}}{3\Delta t} +c\frac{u_{i+2}^{n} - u_{i-2}^{n}}{4\Delta x} = 0\]其中，\(u_i^{n}\)表示时间步\(n\)的网格点\(i\)上的速度或浓度。

这个差分方程可以进一步重写为显式形式：\[u_i^{n+1} = u_i^{n-2} - \frac{3c\Delta t}{4\Deltax}(u_{i+2}^{n} - u_{i-2}^{n})\]这个差分方程可以通过迭代的方式求解，但注意到它的时间步长是\(\Delta t\) ，而其他差分项的时间步长只有 \(\Delta x\)。

为了确保数值稳定性，我们需要满足一个条件：Courant-Friedrichs-Lewy (CFL) 稳定性条件。

CFL条件告诉我们，稳定的数值解需要满足以下条件：\[\frac{c\Delta t}{\Delta x} \leq \frac{1}{2}\]这意味着对流速度与网格尺寸之比需要小于等于1/2上述对流方程的四阶中心差分格式可以很容易地扩展到二维和三维情况。