1、偏微分方程数值解所在学院： 数学与统计学院 课题名称：抛物形扩散方程的有限差分法及数值实例学生姓名： 向聘 1抛物形扩散方程的有限差分法及数值实例1.1 抛物型扩散方程抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程：(1.1.1)2(),0uafxtTt其中 是常数， 是给定的连续函数。按照初边值条件的不同给法，可将a()fx(1.1.1)的定解分为两类：第一，初值问题（Cauchy 问题）：求足够光滑的函数 ，满足方程txu,(1.1.1)和初始条件：, (1.1.2)xu0, x第二，初边值问题（也称混合问题）：求足够光滑的函数 ，满足方txu,程(1.1.1)和初始条件：,
2、(1.1.3)xu0, l及边值条件， (1.1.4),tlt Tt0假定 和 在相应的区域光滑，并且于 ， 两点满足相容条xf,l件,则上述问题有唯一的充分光滑的解。1.2 抛物线扩散方程的求解下面考虑如下热传导方程(1.2.1)2()(0.),0,(uafxtLtx其中， ， , (常数)是扩散系数。0xlTta2取 为空间步长， 为时间步长，其中 ， 是自然数，用两族NlhMTNM平行直线 , 和 , 将矩形域jxN,10 kt,10 G分割成矩形网格。其中 表示网格节点； 表示网tl0; jkxt h格内点（位于开矩形 中的网格节点）的集合； 表示位于闭矩形 中的网GhG格节点的集合；
3、 表示 - 网格边界点的集合。hh表示定义在网点 处的待求近似解， ， 。kju,jkxt Nj0Mk现在对方程进行差分近似：(一) 向前差分格式(1.2.2)kjju1112()kkjjjjjjuaffxh， = =0 (1.2.3)jjjx0ku0N计算后得：(1.2.4)1 1(2)kkkjjjjjurrf 其中， ， ， 。2arh,0Nj ,0M显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下：(1.2.5)100211322433100121()()NNNurruffurruf若记， ，TkNkku121, TNxx121,fxff则显格式(1.
4、2.4) 可写成向量形式(1.2.6)10,0,1kkMuAf其中3rrr21000210 A而对于向前差分格式，当网比 时稳定，当 时不稳定。这就意味12着给定空间步长 以后，时间步长 必须足够小，才能保证稳定。h1.3 抛物型热传导方程数值算例对于(1.2.1)所描述的扩散方程，取 已知方程的精确解为 1a2sintuex：(1.3.1)2(,0)sin()1,.5utxt设空间步长 ，时间步长为 ，网格比 。/hM0.5/N2/rgh向前格式：1112,.,1,.kkkjjjjjuukh边值条件：，1100112, ,kkkkuuj uh.11112, ,kkkkkMMMj初值条件：(,
5、0)sin,01,.jjuxx对时间和空间进行分割，令 M=40，N=1600 ，通过 Matlab 计算得到该方程的解析解，数值解以及相对误差如下：4图(1)解析解的图像图(2) 数值解的图像5图(3)M=40,N=1000 的相对误差的图像我们取部分精确解和数值解进行比较，结果如表(1)xt数值解 精确解 相对误差0.1 0.1 0.1151 0.1152 41.700.2 0.2 0.0815 0.0816 6580.3 0.3 0.0418 0.0419 4.20.4 0.4 0.0183 0.0184 57100.5 0.45 0.0117 0.0118 .30.6 0.35 0.0
6、300 0.0301 460.7 0.25 0.0684 0.0686 1.700.8 0.15 0.5084 0.5085 5890.9 0.05 0.3654 0.3654 .4表 (1)数值解与精确解的比较6由表(1)我们可以看出，精确解和数值解的绝对误差在 以内，因此可以410得出，在分割 M=40，N=1600 下，该有限差分方法对方程(1.3.1)是收敛和稳定的。下面，我们比较在不同的分割下对有限差分算法精度的影响。在扩散系数 不变的情况下，讲时间和空间进行更加细密的分割，取1a，其中，M 表示空间上的分割， N 表示时间上的分割。观察50,N数值解与精确解在节点 处的绝对误差值，
7、如下图所示：,jkxt图(4)M=50,N=10000 的相对误差的图像由图(3)和图(4)，两者在节点处的误差收敛分别是在 和 以内，因此，4105可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。最后，我们比较网比 以及 时扩散方程的收敛情况。1/2rr当网比 时，此时我们取 M=10，N=50，这时，方程的数值解与解析解1r还有相对误差图如下：7图(5)M=10,N=50 的解析解的图像8图(6)M=10,N=50 的数值解的图像图(7)M=10,N=50 的绝对误差的图像此时，我们观察绝对误差发现，扩散
8、方程(1.3.1)时不收敛不稳定的。而前面我们已经知道，到网格比为 时，方程是收敛稳定的。所以，我们可以验12r证，当网比 时稳定，当 时不稳定。12r9参考文献1李荣华，刘播.微分方程数值解法M.北京.高等教育出版社.2009.1.2王曰朋.偏微分方程数值解OL. http:/ Matlab 解法OL.http:/ 数值分析.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*h2/alfa;%时间步长tm=N*tao;%热传导的总时间 tm%t
9、m=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:Nj=2;while j=MU(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+1X(:,i)=x;endfor j=1:M+1Y(j,:)=t;endmesh(X,Y,U);legend(数值解);xlabel(X);ylabel(T);zlabel(U);11z=zeros(M+1,N+1);for j=1:M+1for k=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);for j=1:M+1for k=1:N+1y(j,k)=abs(z(j,k)-U(j,k);endendmesh(x,t,y);legend(绝对误差);xlabel(X);ylabel(Y);zlabel(error);
