一维常系数对流方程的步长定律和固有差分格式

基于一维对流方程差分格式的哲学思考

基于一维对流方程差分格式的哲学思考【摘要】水力学中对流扩散方程都具有通用微分方程的形式，如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）式中α是不为零的常数。

计算程序_计算流体力学_对流方程_有限差分法_Lax格式_迎风格式_FTCS格式

% 一维对流方程迎风格式、Lax格式、FTCS格式差分法计算% 潭花林清华大学航天航空学院% FTCS格式对于一维对流方程不稳定，最好不用clcclear all% 1.参数定义dx=1;x1=-18;x2=18;x=x1:dx:x2;L1=length(x);% dt=0.5*dx; % 收敛dt=2*dx; % 不收敛t1=0;t2=t1+80*dt;t=t1:dt:t2;L2=length(t);alpha=1;lambda=alpha*dt/dx;geshi=1; % 迎风格式% geshi=2; % Lax格式% geshi=3; % FTCS格式% 2.显式求解zeta=zeros(L1,L2);for kk=1:3geshi=kk;for ii=1:L1if x(ii)>0zeta(ii,1)=1;elseif x(ii)==0zeta(ii,1)=1/2;elseif x(ii)<0zeta(ii,1)=0;endendendendif geshi==1for ii=2:L1for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda*(zeta(ii,jj)-z eta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);endzeta1=zeta;else if geshi==2for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=(zeta(ii+1,jj)+zeta(ii-1,jj))/2-. ..lambda/2*(zeta(ii+1,jj)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta2=zeta;else if geshi==3for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda/2*(zeta(ii+1,j j)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta3=zeta;endendendend% 3.绘图% 3.1 t=0figure(1)n=1;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=0时刻的计算结果') %%% 标题axis([-18,18,-0.2,1.2])% 3.2 t=10figure(2)n=(10-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=10s时刻的计算结果') %%% 标题% 3.3 t=20figure(3)n=(20-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') %%% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=20s时刻的计算结果') %%% 标题% 3.4 t=40figure(4)n=(40-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=40s时刻的计算结果') %%% 标题。

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

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

对流扩散方程有限差分方法

对流扩散方程有限差分方法流扩散方程是描述流体内部物质的扩散过程的方程，它可以用于描述溶质的扩散、热量的传导以及动量的传递。

在许多工程和科学领域中，比如地球科学、生物医学和工程学等，流扩散方程都有着广泛的应用。

在数值计算中，有限差分方法是一种常用的数值解法，可以非常有效地解决流扩散方程。

下面将详细介绍对流扩散方程有限差分方法的原理和步骤。

首先，考虑一维流扩散方程的一般形式：∂C/∂t=D∂²C/∂x²-V∂C/∂x其中，C是扩散物质的浓度，t是时间，x是空间位置，D是扩散系数，V是对流速度。

为了使用有限差分方法求解上述方程，我们需要将时间和空间分布离散化，得到方程在网格点上的近似表示。

首先，将时间轴分为n个等间隔的时间步长Δt，空间轴分为m个等间隔的网格点，网格点之间的间距为Δx。

然后，我们使用数值方法来逼近方程中的各个导数项，采用中心差分公式：∂C/∂t≈(C_i^(n+1)-C_i^n)/Δt∂²C/∂x²≈(C_i+1^n-2C_i^n+C_i-1^n)/Δx²∂C/∂x≈(C_i+1^n-C_i-1^n)/(2Δx)将上述近似代入流扩散方程，可以得到：(C_i^(n+1)-C_i^n)/Δt=D(C_i+1^n-2C_i^n+C_i-1^n)/Δx²-V(C_i+1^n-C_i-1^n)/(2Δx)整理上式，可以得到对流扩散方程的有限差分方程：C_i^(n+1)=C_i^n+(DΔt/Δx²)(C_i+1^n-2C_i^n+C_i-1^n)-(VΔt/2Δx)(C_i+1^n-C_i-1^n)上述方程给出了方程在时刻n+1时刻网格点i的值，即C_i^(n+1)，它的值通过已知时刻n时刻各个网格点的值C_i^n来计算。

最后，我们可以使用迭代的方法，从初始条件C_i^0开始，依次计算下一个时刻的网格点C_i^(n+1)，直到达到所需的计算精度或者计算到需要的时间步长。

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

(x, 0) 1 x -1 x 0 这里初始时刻的方程为 (x, 0) x 1 0 x 1 ，即为三角波，其中传播速度为 1，我们 (x, 0) 0 x 1 或 x 1
可以写如下 Fortran90 程序：
!主程序，根据输入选择差分格式，输入1为A格式，2为B格式，3为C格式 !其中输入DX_DT为差分比值，根据输入的比值划分时间网格，其中空间网格为0.05，区间为-8到8 !在A格式中，根据依赖区域，不能计算的左右上三角强行赋值为0 !在B格式中，根据依赖区域右上三角赋值为0 !在C格式中，根据依赖区域左上三角赋值为0 !最后输出为csv格式表格，其中第一列为时间，第二列为空间，第三列为波形数值 !画图时，需要选定同一时间点的二三两列，即可画图 PROGRAM MAIN IMPLICIT NONE !声明差分比值 REAL:: DX_DT !声明差分类型 INTEGER::TYPE_ABC WRITE(*,*)'请输入Dx/Dt，0.5，1或2' READ(*,*)DX_DT WRITE(*,*)'请选择差分格式，A格式输入1，B格式输入2，C格式输入3' READ(*,*)TYPE_ABC !判断类型调用不同的子程序 IF(TYPE_ABC==1) THEN CALL FORMAT_A(DX_DT) ELSE IF(TYPE_ABC==2) THEN CALL FORMAT_B(DX_DT) ELSE IF(TYPE_ABC==3) THEN CALL FORMAT_C(DX_DT) END IF WRITE(*,*)'请在目录下查看结果' PAUSE STOP END
x =1 数值的 t
x =1 为例进行 A，B，C 三种格式的计算结果说明。 t A 格式，根据所的结果，我们绘制不同时刻的波形图像，选取时刻为初始时刻，第 35 个时间区间时刻（3.5s），和最后的 80 区间时刻（8s），图像如下：

一维定常对流扩散反应方程的高精度紧致差分格式

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

一类对流扩散方程的特征—修正差分格式及最大模估计

界处的值，到了ｏ（＋＾阶的最大模误差估计．为了提高空间的逼近阶，个自然的想得）一法是采用高次插值．１中利用二次插值得到了关于ｚ模ｏ（＋ｈ）收敛阶，没有给出［］。。的却
ＤｕＮｉｇｎ
（ｐ．ｏａｈＤｅｔｆＭｔ．，ＳａｇｏｇＵｎｖｒｉｙ，Ｊｎｎ２００）ｈｎｄｎｉｅｓｔｉａ５１０
ＡｂｓｒｃＡｗｅｈｎｑｕｅｉｎｔｏｃｄｔｔａｔｎｅｔｃｉｓｉｒｄｕｅｏａｐｐｒｘｉａｅｔａｌｓｏｆｔｎｔｒｏｍｔｈｅｖｕｅｈｅｉｅ — ｐｏａｔｄｐｏｎｔｎｓｅｄｏｉｏａｌｅｉｓｉｔａｆｕｓｎｇｌｃｌｑｕａａｉｎｔｒｄｒｔｃｉｅｐｏｌｔｏａｉｎ．Ａｏｆｅｃａｒｃｍｄｉｉｄｈａ — ｔｒｓｉｆｅｒｎｃｃｍｅｂｓｄｏｎｔｅｔｃｅｉｔｃｄｉｆｅｅｓｈｅａｅｈｅｎｗｅｈｎｉｅｉｏｒｕｌｔｄｔｒａｔａｃｑｕｓｆｍａｅｏｔｅｏｎ—
作效率．特征法与传统的数值计算方法（分法，限元法等）结合，解决一类广泛的科差有相可
学与工程计算中的实际问题口．已有文献大多分析了数值解关于Ｌ。ｚ）Ｈ＾）的收敛性，得到了很好的收敛性结（。，（模并果．但是，于最大模的相应分析却并不多见．１中采用线性插值来逼近特征线与网格交关［］

计算程序_计算流体力学_对流方程_有限差分法_Lax格式_迎风格式_FTCS格式

% 一维对流方程迎风格式、Lax格式、FTCS格式差分法计算% 潭花林清华大学航天航空学院% FTCS格式对于一维对流方程不稳定，最好不用clcclear all% 1.参数定义dx=1;x1=-18;x2=18;x=x1:dx:x2;L1=length(x);% dt=0.5*dx; % 收敛dt=2*dx; % 不收敛t1=0;t2=t1+80*dt;t=t1:dt:t2;L2=length(t);alpha=1;lambda=alpha*dt/dx;geshi=1; % 迎风格式% geshi=2; % Lax格式% geshi=3; % FTCS格式% 2.显式求解zeta=zeros(L1,L2);for kk=1:3geshi=kk;for ii=1:L1if x(ii)>0zeta(ii,1)=1;elseif x(ii)==0zeta(ii,1)=1/2;elseif x(ii)<0zeta(ii,1)=0;endendendendif geshi==1for ii=2:L1for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda*(zeta(ii,jj)-z eta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);endzeta1=zeta;else if geshi==2for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=(zeta(ii+1,jj)+zeta(ii-1,jj))/2-. ..lambda/2*(zeta(ii+1,jj)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta2=zeta;else if geshi==3for ii=2:(L1-1)for jj=1:(L2-1)zeta(ii,jj+1)=zeta(ii,jj)-lambda/2*(zeta(ii+1,j j)-zeta(ii-1,jj));endzeta(1,jj+1)=zeta(2,jj+1);zeta(L1,jj+1)=zeta(L1,jj)-lambda*(zeta(L1,jj)-z eta(L1-1,jj));endzeta3=zeta;endendendend% 3.绘图% 3.1 t=0figure(1)n=1;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=0时刻的计算结果') %%% 标题axis([-18,18,-0.2,1.2])% 3.2 t=10figure(2)n=(10-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=10s时刻的计算结果') %%% 标题% 3.3 t=20figure(3)n=(20-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') %%% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=20s时刻的计算结果') %%% 标题% 3.4 t=40figure(4)n=(40-t(1))/dt;plot(x,zeta1(1:L1,n),'-k',x,zeta2(1:L1,n),'-.k' ,x,zeta3(1:L1,n),'--k')% 作图% axis equal %%% 是否要求x、y坐标间距相等% grid on %%% 是否要求画网格xlabel('x/m'),ylabel('t/s') % %% x,y轴表示的变量含义%text(1,2,'f(x)') %%% 图中文字标识legend('迎风格式','Lax格式','FTCS格式') %%% 不同曲线的线型区分title('t=40s时刻的计算结果') %%% 标题聚乙烯（PE）简介1.1聚乙烯化学名称：聚乙烯英文名称：polyethylene，简称PE结构式：聚乙烯是乙烯经聚合制得的一种热塑性树脂，也包括乙烯与少量α-烯烃的共聚物。

一维对流方程的迎风格式

一维对流方程的迎风格式标题：一维对流方程的迎风格式简介：本文将介绍一维对流方程的迎风格式，探讨其数值求解方法，并分析其在实际问题中的应用。

正文：一维对流方程是描述流体运动中物质传输的基本方程之一，其数值求解在科学计算和工程应用中具有重要意义。

迎风格式是一种常用的求解一维对流方程的数值方法，其基本原理是根据流体速度来调整网格上的近似解。

在迎风格式中，我们首先将一维对流方程离散化，得到差分格式。

然后，根据网格上的迎风方向速度来选择合适的近似解，以确保数值解的稳定性和精确性。

具体而言，对于一维对流方程u_t+au_x=0，我们假设网格点(xi,tn)处的解为ui^n。

利用向前差分和向后差分的组合，可以得到迎风格式的数值离散方程：ui^{n+1}=ui^n-a(dt/dx)(ui^n-ui^{n-1})其中，a为对流速度，dt为时间步长，dx为空间步长。

通过迎风格式的数值求解，我们可以得到在下一个时间步的解ui^{n+1}。

迎风格式的特点是它能够处理对流方程中的激波和间断等现象。

由于迎风格式中的近似解取决于流体速度，它能够更好地捕捉到流体中物质传输的特点。

因此，迎风格式在流体动力学、气象学和地质学等领域中得到广泛应用。

需要注意的是，迎风格式的稳定性和精确性与时间步长和空间步长的选择有关。

较小的时间步长和空间步长能够提高数值解的精确性，但也会增加计算量。

因此，在实际应用中，需要综合考虑计算效率和数值解的准确性。

总结起来，一维对流方程的迎风格式是一种常用的数值求解方法，能够有效地处理流体中的物质传输问题。

在进行数值计算时，我们需要合理选择时间步长和空间步长，以确保数值解的稳定性和精确性。

迎风格式在科学计算和工程应用中具有广泛的应用前景，对于解决一维对流方程相关问题具有重要意义。

注意：本文旨在提供关于一维对流方程的迎风格式的基本介绍，并未涉及具体的数值实现和应用案例。

如需深入研究和应用，请参考相关专业文献和资料。

一维对流方程在A、B、C三种差分格式

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

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

并作相关讨论。

二、实验原理1x1x0)(x,1x10)(x,x1-x10)(x,t8x8-,0>-<=≤≤+-=≤≤+=>≤≤=∂∂+∂∂或ζζζζζxxt三、上机要求：1.学会对MS-FORTRAN的基本操作。

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

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

四、实验程序以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_dtwrite(*,*) "选择格式，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) thenopen(unit=8,file='out_c.csv')endifr_t=160/dx_dt !计算时间网格总数allocate(s(r_t+1,321)) !分配存储矩阵的空间!第一层赋初值do i=1,140,1s(1,i)=0write(8,*)"1",",",i,",",s(1,i)end dodo i=141,161,1s(1,i)=1+0.05*(i-161)write(8,*)"1",",",i,",",s(1,i)end dodo i=162,181,1s(1,i)=1-0.05*(i-161)write(8,*)"1",",",i,",",s(1,i)end dodo i=182,321,1s(1,i)=0write(8,*)"1",",",i,",",s(1,i)end do!循环时间层,根据格式的选择来判断执行哪一部分if(abc==1) thendo i=2,r_t,1do j=i,322-i,1s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j-1))/(dx_dt*2)write(8,*)i,',',j,',',s(i,j)end dodo k=1,i-1,1 !余下部分赋值0,下同s(i,k)=0write(8,*)i,',',k,',',s(i,k)end dodo k=322-i,321,1s(i,k)=0write(8,*)i,',',k,',',s(i,k)end doend doelseif(abc==2) thendo i=2,r_t+1,1do j=1,322-i,1s(i,j)=s(i-1,j)-(s(i-1,j+1)-s(i-1,j))/dx_dtwrite(8,*)i,',',j,',',s(i,j)end dodo k=322-i,321,1s(i,k)=0write(8,*)i,',',k,',',s(i,k)end doend doelseif(abc==3) thendo i=2,r_t+1,1do j=i,321,1s(i,j)=s(i-1,j)-(s(i-1,j)-s(i-1,j-1))/dx_dtwrite(8,*)i,',',j,',',s(i,j)end dodo k=1,i-1,1s(i,k)=0write(8,*)i,',',k,',',s(i,k)end doend doendif!完成提示write(*,*)'数据已输出至源目录'pausestopend program五、实验结果及分析程序运行后在对应目录下生成csv表格文件，根据输入的的值不同生成对应的网格并计算各节点数值。

一维对流方程的差分格式

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

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

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

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

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

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

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

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

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

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

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

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

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

对流方程及算法介绍

1 引言2 对流方程及算法介绍2.1对流方程的概述对流：是指由于流体的宏观运动，从而使流体各部分之间发生相对位移，冷热流体相互掺混所引起的热量传递过程。

对流仅发生在流体中，对流的同时必伴随有导热现象。

人们研究对流扩散方程，主要的研究对象是流体在流动过程中，流体所携带的某种物质的物理量的变化规律，例如传热过程中温度的变化规律或者溶解于流体中溶质的物质浓度等物理量的变化规律。

这些变化通常包括对流、扩散以及由于某种物理或者化学的因素而引起的物理量的自身衰减或增长。

最简单的一维对流扩散方程形如(2-1)式: （2-1）其中C 是常数，它属于双曲型方程，可以被用来描述流体的运动等物理现象。

2.2水对流现象的简易演示2.2.1 基本步骤用两只相同的小烧杯，各装上冷水，再如图1所示插入长短两根吸管，虹吸管由普通化学实验用玻璃管在酒精灯上加热弯成，一根查到被子底部，一根只插入水的表面，再在右杯中滴入几滴墨水并搅拌均匀，现在开始用酒精灯加热左边的烧杯，一段时间后就可以明显的看到染了颜色的水从右杯源源不断的流入左杯，左杯的水源源不断的流入右杯，最后两杯水都变成了墨水的颜色，与此同时用手摸右边的杯子，右边的水也热了起来，这就是冷热水发生了对流的缘故。

2.2.2 实验注意事项短虹吸管只插入水的表面，不能过深。

玻璃管宜选壁较厚一些的，这样绝热性好一些，效果也好一些。

0=∂∂+∂∂xuC t u2.2.3 实验原理分析对左边的水杯用酒精灯加热，水受热密度变小开始上升，右边水杯的冷水从下边的吸管流向左边的水杯进行补充，左边水杯的热水从上边的吸管流向右边的水杯，这样一会儿两杯水都变成墨水的颜色了[1]。

在冷水里面掺热水也是一样的道理，在不搅拌的情况下，最后水温基本都是一个温度，这就是水的对流，除了水的对流还有刮风是空气的对流，气压高的一方向气压低的一方补充空气，这就形成了对流，就会产生风；还有冬天在家里开空调，形成空气对流，最后整个房间的温度都升了起来。

基于一维对流方程差分格式的哲学思考

基于一维对流方程差分格式的哲学思考
李占松
【期刊名称】《科技视界》
【年(卷),期】2013(000)033
【摘要】水力学中对流扩散方程都具有通用微分方程的形式，如N-S方程组中的运动方程等。

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

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

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

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

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

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

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

【总页数】2页(P198-199)
【作者】李占松
【作者单位】郑州大学水利与环境学院，河南郑州 450001
【正文语种】中文
【相关文献】
1.一维非定常对流扩散方程非均匀网格上的高精度紧致差分格式 [J], 黄雪芳;郭锐;葛永斌
2.一维非定常对流扩散反应方程的高精度紧致差分格式 [J], 杨晓佳;田芳
3.一维定常对流扩散反应方程的高精度紧致差分格式 [J], 祁应楠;武莉莉
4.一维对流扩散方程Crank-Nicolson特征紧致差分格式 [J], 鞠萍;杨青
5.求解一维对流方程的高精度紧致差分格式 [J], 侯波;葛永斌
因版权原因，仅展示原文概要，查看原文内容请购买。

有限差分求解1维流体方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

差分方程通用

程序如下：a=input('差分方程左端的系数向量a=[a(1),...a(na)] a=');b=input('差分方程右端的系数向量b=[b(1),...b(na)] b=');u=input('输入信号序列u=');na=length(a);nb=length(b);nu=length(u);s=['起算点前',int2str(na-1),'点y的值=[y(',int2str(na-2),'),...y(0)]=']; ym=zeros(1,na+nu-1);ym(1:na-1)=input(s);um=[zeros(1,na-2),u];for n=na:na+nu-1ys=ym(n-1:-1:n-na+1);us=um(n:-1:n-nb+1);ym(n)=(b*us'-a(2:na)*ys')/a(1);endy=ym(na:na+nu-1);sten(y),grid on;line([0,nu],[0,0]);运行程序输入差分方程左端的系数向量a=[a(1),...a(na)] a=[1,0.1,0.15,-0.225] 差分方程右端的系数向量b=[b(1),...b(na)] b=[3,7,1]输入信号序列u=exp(0.1*[1:20])起算点前3点y的值=[y(2),...y(0)]=[0,0,0]??? Index exceeds matrix dimensions.Error in ==> ca at 12us=um(n:-1:n-nb+1);1 一维对流方程1.1 一维对流方程的形式：，其中，u代表物质的量，a代表物质的运动速度。

此一维对流方程仅仅表示物质的运动情况，而与边界条件或是约束条件无关。

当a为常数时，此一维对流方程为一维常系数对流方程，当a不为常数时，方程为一维变系数对流方程。

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