一维对流扩散方程的数值解法

1一维对流扩散方程的数值解法对流-扩散方程是守恒定律控制方程的一种模型方程，它既是能量方程的表示形式，同时也可以认为是把压力梯度项隐含到了源项中去的动量方程的代表。因此，以对流-扩散方程为例，来研究数值求解偏微分方程的相容性、收敛性和稳定性具有代表性的意义。1数学模型本作业从最简单的模型方程，即一维、稳态、无源项的对流扩散方程出发，方程如下：22,02fffUDxtxx(1)初始条件,0sin(2)fxtAkx(2)解析解224,sin2DktfxteAkxUt(3)式中，1,0.05,0.5,1UDAk函数(3)描述的是一个衰减波的图像，如图1所示t=0t=0.5.t=1图1函数224,sin2DktfxtekxUt的图像（U=1，D=0.05，k=1）2数值解法2.1数值误差分析在网格点,in上差分方程的数值解nif偏离该点上相应的偏微分方程的精确解,fin的值，称为网格节点上的数值误差。当取定网格节点数21N时，观察差分方程的解与微分方程的解在不同时间步长下的趋近程度，其中时间步长分别取值0.05,0.025,0.0125,0.0005t。2（a）21,0.05Nt（b）21,0.025Nt（c）21,0.0125Nt.（d）201,0.0005Nt图2数值误差随步长的变化情况从图2的(a)~(d)可以定性的看出，数值误差与步长的大小有关。在满足稳定性条件的前提下，数值误差随着时间步长的减小而减小，同时，图（d）表示增大网格的分辨率也有助于减小网格误差。为了对数值误差有一个定量的认识，接下来取定时间步长为0.0005t，分别算出11,21,41,61,81,101,121,161N时，指标21NjexactjEhff的数值，如表1所示。表1不同网格节点数下指标E的值网格点数21NjexactjEhff11N0.062121N0.011041N0.001961N6.4657e-0481N3.0793e-04101N1.8223e-04121N1.2965e-04161N9.4418e-053由表1可以看出，指标E的值随着网格节点数的增加呈先大幅减小后平缓变化的趋势，该趋势示于图3中。（a）（b）图3不同网格节点数下指标E的值在图3的（a）和（b）中，随着网格节点数的增加，曲线最后都趋于水平，即此时再增加网格节点数，对于提高数值求解的精确性作用不大。2.2截断误差分析对于方程5fxx(4)其一阶导数和二阶导数的精确解分别为45fxx和23220fxx。一阶导数和二阶导数的中心差分表示式分别为211()2iifffohxh和2211222()iiiffffohxh，2()oh为截断误差。当空间步长分别取值为0.2,0.1,0.05,0.025,0.0125h时，节点1ix取处一阶导数和二阶导数的精确解与离散解的差值见表2。表2截断误差随空间步长的变化h'1fx''1fx0.20.40160.40000.10.10010.10000.050.02500.02500.0250.00630.00620.01250.00160.0016图示如下4图4截断误差随空间步长的变化由图4可以看出，截断误差随着步长的减小而减小，当步长趋近于0时，截断误差趋近于0。3结论（1）在满足稳定性条件的前提下，数值误差随着时间步长的减小而减小。而对于空间步长的减小，数值误差呈先大幅减小后变平缓的趋势，即存在一个网格无关解的网格节点数。（2）对于导数在空间上的离散，二阶精度的中心差分格式能较好的满足数值求解精确性的要求，当选取合适的步长时，其精确性也将显著提高。5附录1%%===============================================================clearclcformatlong%%===============================================================U=1;D=0.05;k=1;A=0.5;T=0.5;dt=0.05;m=T/dt;L=2;N=21;h=L/(N-1);%%===============================================================E0=0;E=0;%误差估计%%===============================================================x=linspace(0,L,N);y=[1:N];fori=1:Nf0(i)=A*sin(2*pi*k*x(i));%初始条件end%%===============================================================foris=1:m,is%%===============================================================fori=2:N-1f(i)=f0(i)-(U*dt/(2*h))*(f0(i+1)-f0(i-1))+(D*dt/(h^2))*(f0(i+1)-2*f0(i)+f0(i-1));end%f(1)=f0(1)-(U*dt/(2*h))*(f0(2)-f0(N-1))+(D*dt/(h^2))*(f0(2)-2*f0(1)+f0(N-1));%边界条件%f(N)=f0(N)-(U*dt/(2*h))*(f0(2)-f0(N-1))+(D*dt/(h^2))*(f0(2)-2*f0(N)+f0(N-1));%边界条件f(1)=exp(-4*D*k^2*pi^2*dt*is)*A*sin(2*pi*k*(0-U*dt*is));%边界条件f(N)=exp(-4*D*k^2*pi^2*dt*is)*A*sin(2*pi*k*(L-U*dt*is));%±边界条件f0=f;end%%===============================================================fori=1:Nfexact(i)=exp(-4*D*k^2*pi^2*T)*A*sin(2*pi*k*(x(i)-U*T));E0=(f(i)-fexact(i))^2+E0;endE=h*sqrt(E0)%%===============================================================f1=exp(-4*D*k^2*pi^2*T)*A*sin(2*pi*k*(x-U*T));plot(y,f1,'r',y,f,'b','linewidt',2),legend('Exact','Numerical'),axis([1N-22]),gridon%plot(y,f1,'r','linewidt',4),legend('Exact'),axis([1N-22]),gridon%holdon%plot(y,f,'b','linewidt',2),legend('Numerical')6附录2%====================================================================i=1;%Çói(x=1)½ÚµãÉϵĵ¼ÊýÖµh=[0.2,0.1,0.05,0.025,0.0125];%====================================================================forn=1:5f1(n)=((i+h(n))^5-(i-h(n))^5)/(2*h(n));%f=x^5;f2(n)=((i+h(n))^5-2*i^5+(i-h(n))^5)/(h(n))^2;E1(n)=f1(n)-5*i^4;E2(n)=f2(n)-20*i^2;end%====================================================================loglog(1./h,E1,'o-',1./h,E2,'*-'),xlabel('1/h'),ylabel('Error'),...legend('first-orderderivative','second-orderderivative')