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

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

2、滑，并且于，两点满足相容条件,则上述问题有唯一的充分光滑的解。1.2抛物线扩散方程的求解下面考虑如下热传导方程 (1.2.1)其中，,(常数)是扩散系数。取为空间步长，为时间步长，其中，是自然数，用两族平行直线, 和, 将矩形域分割成矩形网格。其中 表示网格节点；表示网格内点（位于开矩形中的网格节点）的集合；表示位于闭矩形中的网格节点的集合；表示-网格边界点的集合。表示定义在网点处的待求近似解，。现在对方程进行差分近似：(一) 向前差分格式 (1.2.2)， =0 (1.2.3)计算后得： (1.2.4)其中，。显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如

3、下： (1.2.5)若记，则显格式(1.2.4)可写成向量形式 (1.2.6)其中而对于向前差分格式，当网比时稳定，当时不稳定。这就意味着给定空间步长以后，时间步长必须足够小，才能保证稳定。1.3抛物型热传导方程数值算例对于(1.2.1)所描述的扩散方程，取已知方程的精确解为 ： (1.3.1)设空间步长，时间步长为，网格比。向前格式： 边值条件：，.初值条件： 对时间和空间进行分割，令M=40，N=1600，通过Matlab计算得到该方程的解析解，数值解以及相对误差如下：图(1)解析解的图像图(2)数值解的图像图(3)M=40,N=1000的相对误差的图像我们取部分精确解和数值解进行比较，结

4、果如表(1)数值解精确解相对误差0.10.10.11510.1152 0.20.20.08150.08160.30.30.04180.04190.40.40.01830.01840.50.450.01170.01180.60.350.03000.03010.70.250.06840.06860.80.150.50840.50850.90.050.36540.3654表(1)数值解与精确解的比较由表(1)我们可以看出，精确解和数值解的绝对误差在以内，因此可以得出，在分割M=40，N=1600下，该有限差分方法对方程(1.3.1)是收敛和稳定的。下面，我们比较在不同的分割下对有限差分算法精度的影响

5、。在扩散系数不变的情况下，讲时间和空间进行更加细密的分割，取，其中，M表示空间上的分割，N表示时间上的分割。观察数值解与精确解在节点处的绝对误差值，如下图所示：图(4)M=50,N=10000的相对误差的图像由图(3)和图(4)，两者在节点处的误差收敛分别是在和以内，因此，可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。最后，我们比较网比以及时扩散方程的收敛情况。当网比时，此时我们取M=10，N=50，这时，方程的数值解与解析解还有相对误差图如下：图(5)M=10,N=50的解析解的图像图(6)M=10,N=

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

7、;tao=lambda*h2/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的关系*%for k=1:N j=2; while j<=M U(j,k+1)=lambda*U(j+1,k)+(1-2*lambda)*U(j,k)+lambda*U(j-1,k); j=j+1; endend length(U);%*设置立体网格*%for i=1:N+1 X(:,i

8、)=x;endfor j=1:M+1 Y(j,:)=t;endmesh(X,Y,U);legend('数值解');xlabel('X');ylabel('T');zlabel('U');z=zeros(M+1,N+1);for j=1:M+1 for k=1:N+1 z(j,k)=exp(-pi*pi*t(k)*sin(pi*x(j); endend%mesh(x,t,z')legend('解析解');xlabel('X');ylabel('T');zlabel('Z');for j=1:M+1 for k=1:N+1 y(j,k)=abs