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

偏微分方程数值解所在学院：数学与统计学院课题名称：抛物形扩散方程的有限差分法及数值实例学生姓名：向聘1抛物形扩散方程的有限差分法及数值实例1.1抛物型扩散方程抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程：22(),0uuafxtTtx(1.1.1)其中a是常数，()fx是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1)的定解分为两类：第一，初值问题（Cauchy问题）：求足够光滑的函数txu,，满足方程(1.1.1)和初始条件：xxu0,,x(1.1.2)第二，初边值问题（也称混合问题）：求足够光滑的函数txu,，满足方程(1.1.1)和初始条件：xxu0,,0xl(1.1.3)及边值条件0,,0tlutu，Tt0(1.1.4)假定xf和x在相应的区域光滑，并且于0,0，0,l两点满足相容条件,则上述问题有唯一的充分光滑的解。1.2抛物线扩散方程的求解下面考虑如下热传导方程22()(0.)(,)0(,0)()uuafxtxutuLtuxx(1.2.1)其中，0xl，Tt0,a(常数)是扩散系数。取Nlh为空间步长，MT为时间步长，其中N，M是自然数，用两族2平行直线jhxxj,Nj,,1,0和kttk,Mk,,1,0将矩形域GTtlx0;0分割成矩形网格。其中,jkxt表示网格节点；hG表示网格内点（位于开矩形G中的网格节点）的集合；hG表示位于闭矩形G中的网格节点的集合；h表示hG-hG网格边界点的集合。kju表示定义在网点,jkxt处的待求近似解，Nj0，Mk0。现在对方程进行差分近似：(一)向前差分格式kjkjuu11122(())kkkjjjjjjuuuafffxh(1.2.2)jjjxu0，ku0=kNu=0(1.2.3)计算后得：111(12)kkkkjjjjjurururuf(1.2.4)其中，2arh，1,,1,0Nj，1,,1,0Mk。显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下：10001210110002321210003432310001121(12)(12)(12)(12)NNNNNurururufurururufurururufurururuf(1.2.5)若记TkNkkkuuu121,,,u，TNxxx121,,,，TNxfxfxf121,,,f则显格式(1.2.4)可写成向量形式10,0,1,,1kkkMuAufu(1.2.6)其中3rrrrrrrrrr21002100210021A而对于向前差分格式，当网比12r时稳定，当12r时不稳定。这就意味着给定空间步长h以后，时间步长必须足够小，才能保证稳定。1.3抛物型热传导方程数值算例对于(1.2.1)所描述的扩散方程，取1a已知方程的精确解为2sintuex：22(,0)sin()(0,)(1,)001,00.5uutxuxxututxt(1.3.1)设空间步长1/hM，时间步长为0.5/N，网格比2/rgh。向前格式：11122,1,...,1,1,...,kkkkkjjjjjuuuuujMkNh边值条件：110010111220,,kkkkkkkuuuuujuuh，11111122,,kkkkkkkMMMMMMMuuuuujMuuh.初值条件：(,0)sin,0,1,...,jjuxxjM对时间和空间进行分割，令M=40，N=1600，通过Matlab计算得到该方程的解析解，数值解以及相对误差如下：4图(1)解析解的图像图(2)数值解的图像5图(3)M=40,N=1000的相对误差的图像我们取部分精确解和数值解进行比较，结果如表(1)xt数值解精确解相对误差0.10.10.11510.115241.1700100.20.20.08150.081641.6580100.30.30.04180.041941.2752100.40.40.01830.018457.4456100.50.450.01170.011855.3755100.60.350.03000.030141.0674100.70.250.06840.068641.7410100.80.150.50840.508557.5897100.90.050.36540.365451.740710表(1)数值解与精确解的比较6由表(1)我们可以看出，精确解和数值解的绝对误差在410以内，因此可以得出，在分割M=40，N=1600下，该有限差分方法对方程(1.3.1)是收敛和稳定的。下面，我们比较在不同的分割下对有限差分算法精度的影响。在扩散系数1a不变的情况下，讲时间和空间进行更加细密的分割，取50,10000MN，其中，M表示空间上的分割，N表示时间上的分割。观察数值解与精确解在节点,jkxt处的绝对误差值，如下图所示：图(4)M=50,N=10000的相对误差的图像由图(3)和图(4)，两者在节点处的误差收敛分别是在410和510以内，因此，可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。最后，我们比较网比1/2r以及1r时扩散方程的收敛情况。当网比1r时，此时我们取M=10，N=50，这时，方程的数值解与解析解还有相对误差图如下：7图(5)M=10,N=50的解析解的图像图(6)M=10,N=50的数值解的图像8图(7)M=10,N=50的绝对误差的图像此时，我们观察绝对误差发现，扩散方程(1.3.1)时不收敛不稳定的。而前面我们已经知道，到网格比为12r时，方程是收敛稳定的。所以，我们可以验证，当网比12r时稳定，当12r时不稳定。9[参考文献][1]李荣华，刘播.微分方程数值解法[M].北京.高等教育出版社.2009.1.[2]王曰朋.偏微分方程数值解[OL].=-UAAXRI8_V547zVS6UwenT8rGKsbwNf9CuDmh2qmsy5K8eQ32PzhYazZ_sWiHfz_Pj3LA7ufHH4O3t8NIlUCnXNUiyYhssqmNBeArsrLwQG[3]未知.偏微分方程的Matlab解法[OL].=WBlR0q9n02YsNg6UOLkCQ4Z5QOCefDKNdEglrNR2TJ8VlxaJ9IkCrLlaEDlJ_OHCLehf19UJU_B7PRe1skQTYaaLcAa1UWcwlv7GYsWNtG7[4]周品，何正风.MATLAB数值分析.[M].北京.机械工业出版社.2009.1.10附录：L=1;M=40;N=1600;alfa=1;lambda=0.5;%网格比%**********************************************%h=L/M;%空间步长x=0:h:L;x=x';tao=lambda*h^2/alfa;%时间步长tm=N*tao;%热传导的总时间tm%tm=0.1;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，与杆长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+111fork=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');forj=1:M+1fork=1:N+1y(j,k)=abs(z(j,k)-U(j,k));endendmesh(x,t,y');legend('绝对误差');xlabel('X');ylabel('Y');zlabel('error');