抛物形扩散方程的有限差分法及数值实例

1、偏微分方程数值解所在学院:数学与统计学院课题名称:抛物形扩散方程的有限差分法及数值实例学生姓名：向聘抛物形扩散方程的有限差分法及数值实例抛物型扩散方程抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程：2auf(x),0tT(1.1.1)tx其中a是常数，f(x)是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1)的定解分为两类：第一,初值问题(Cauchy问题)：求足够光滑的函数ux,t,满足方程(1.1.1)和初始条件：ux,0x,x(1.1.(2)第二，初边值问题(也称混合问题)：求足够光滑的函数ux,t,满足方程(1.1.1)和初始条件：ux,0x,0xl(1.1.(

2、3)及边值条件u0,tul,t0,0tT(1.1.(4)假定fx和x在相应的区域光滑，并且于0,0,l,0两点满足相容条件，则上述问题有唯一的充分光滑的解。抛物线扩散方程的求解卜面考虑如下热传导方程取h2为空间步长,N平行直线XXjjh,G0xl;0tT网格内点（位于开矩形uatu(0.t)u(x,0)2uf(X)Xu(L,t)0(x)T,a（常数）是扩散系数工为时间步长,M0,1,N和ttk分割成矩形网格。其中其中M是自然数，用两族0,1,M将矩形域Xj,tk表示网格节点；Gh表示G中的网格节点）的集合;Gh表示位于闭矩形G中的网格节点的集合；h表示Gh-Gh网格边界点的集合。u：表示定义在

3、网点Xj,tk处的待求近似解，0jN,现在对方程进行差分近似:(一)向前差分格式k1kUjUjkuj1a2u：u：1h2fj(fjf(Xj)(1.2.2)计算后得:其中，r0UjXjuk=uN=0(1.2.3)k1uj01,kruj1(12r)ukkruj1fj(1.2.4),N1,k0,1,显然，这是一个四点显示格式,每一层各个节点上的值是通过一个方程组求解到的。方程组如下:1ui1u21u3rururu020304(1(1(12r)u°2r)u02r)u:rururu000102f1(1.2.5)若记1un1ruN(12r)u；10ruNfN1kkkkTuUi,U2,Un1,而对

4、于向前差分格式，当网比r1时稳定，当r21一.-时不稳定。这就意味着2给定空间步长h以后，时间步长必须足够小，才能保证稳定ffXi,fX2,fXn1T则显格式（1.2.4）可写成向量形式(1.2.6)uk1Aukf,k0,1,JU，M10u其中12rr00r12rrA00r12rr00r12r抛物型热传导方程数值算例对于（1.2.1）所描述的扩散方程，取a1已知方程的精确解为-2tuesinx2uutX2u(x,0)sin(x)u(0,t)u(1,t)00x1,0t0.5设空间步长h1/M,时间步长为向前格式：k1kk1kkujujuj2ujuj12,j1,.,Mh(1.3.1)0.5/N,网

5、格比rg/h21,k1,.,N边值条件:k1kk1八kku°u0u12u°u1kkj0,2,u15,hk1kk1kkUMUMUM12UMUM1kkJM,2,UM1UM1.h初值条件：u(xj,0)sinxJ,J0,1,.,M对时间和空间进行分割，令M=40N=160Q通过Matlab计算得到该方程的解析解，数值解以及相对误差如下：3.一一II解析解:.一:*':，：IJ.-"-!"'"1一“:r'08T00X图(1)解析解的图像图(2)数值解的图像图(3)M=40,N=1000的相对误差的图像我们取部分精确解和数值解进行

6、比较，结果如表Xt数值解精确解相对误差41.1700101.6580104_41.2752107.4456105一一55.375510-41.06741041.7410107.58971051.7407105表(1)数值解与精确解的比较由表(1)我们可以看出，精确解和数值解的绝对误差在104以内，因此可以得出，在分割M=40N=1600下，该有限差分方法对方程(1.3.1)是收敛和稳定向下面，我们比较在不同的分割下对有限差分算法精度的影响。在扩散系数a1不变的情况下，讲时间和空间进行更加细密的分割，取M50,N10000,其中，M表示空间上的分割，N表示时间上的分割。观察数值解与精确解在节点X

7、j,tk处的绝对误差值，如下图所示：图(4)M=50,N=10000的相对误差的图像维时误差由图(3)和图(4),两者在节点处的误差收敛分别是在104和105以内，因此,可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。最后，我们比较网比r1/2以及r1时扩散方程的收敛情况。当网比r1时，此时我们取M=10N=50,这时，方程的数值解与解析解还有相对误差图如下:图(5)M=10,N=50的解析解的图像数值解图(6)M=10,N=50的数值解的图像绝对误差x10Y00图(7)M=10,N=50的绝对误差的图像此时

8、，我们观察绝对误差发现，扩散方程(1.3.1)时不收敛不稳定的。而前1面我们已经知道，到网格比为r1证，当网比r一时稳定，当r时不稳定。2时，方程是收敛稳定的。所以，我们可以验2参考文献1李荣华，刘播.微分方程数值解法M.北京.高等教育出版社2王曰朋.偏微分方程数值解OL.未知.偏微分方程的Matlab解法OL.周品，何正风.MATLABa值分析.M.北京.机械工业出版社.附录:L=1;M=40;N=1600;alfa=1;lambda=;%网格比%*%h=L/M;%空间步长x=0:h:L;x=x'tao=lambda*hA2/alfa;%时间步长tm=N*tao;%热传导的总时间tm

9、%tm=;t=0:tao:tm;t=t'%十算初值和边值U=zeros(M+1,N+1);U(:,1)=sin(pi*x);U(1,:)=0;U(M+1,:)=0;%*用差分法求出温度U,与卞f长L,时间t的关系*%fork=1:Nj=2;whilej<=MU(j,k+1)=lambda*U(j+1,k)+(1-2*lambda)*U(j,k)+lambda*U(j-1,k);j=j+1;endendlength(U);%*设置立体网格*%fori=1:N+1X(:,i)=x;endforj=1:M+1Y（j,：）二t；endmesh(X,Y,U);legend('数值解');xlabel('X');ylabel('T');zlabel('U');z=zeros(M+1,N+1);forj=1:M+1fork=1:N+1z(j,k)=exp(-pi*pi*t(k)*sin(pi*x(j);endend%mesh(x,t,z')legend('解析解');xlabel('X');ylabel('T');zlabel('Z');f