偏微分方程数值解 所在学院： 数学与统计学院 课题名称：抛物形扩散方程的有限差分法及数值实例 学生姓名： 向聘

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

抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程： 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两点满足相容条件,则上述问题有唯一的充分光滑的解。

抛物线扩散方程的求解 下面考虑如下热传导方程 22()(0.)(,)0(,0)()uuafxtxutuLtuxx











(1.2.1)

其中，0xl，Tt0,a(常数)是扩散系数。 取Nlh为空间步长，MT为时间步长，其中N，M是自然数，用两族平行直线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。 显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下： 100012101100023212100034323

10001121(12)(12)(12)(12)NNNNNurururufurururufurururufurururuf













(1.2.5)

若记 TkNkkkuuu121,,,u，TNxxx121,,,，TNxfxfxf121,,,f

则显格式(1.2.4)可写成向量形式 10,0,1,,1kkkM

uAuf

u (1.2.6)

其中





rrrrrrrrrr21002100210021



A

而对于向前差分格式，当网比12r时稳定，当12r时不稳定。这就意味着给定空间步长h以后，时间步长必须足够小，才能保证稳定。

抛物型热传导方程数值算例 对于(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计算得到该方程的解析解，数值解以及相对误差如下：

图(1)解析解的图像 图(2)数值解的图像

图(3)M=40,N=1000的相对误差的图像 我们取部分精确解和数值解进行比较，结果如表(1) x t

数值解 精确解 相对误差

41.170010

41.658010

41.275210

57.445610

55.375510

41.067410

41.741010

57.589710

51.740710

表(1)数值解与精确解的比较 由表(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，这时，方程的数值解与解析解还有相对误差图如下： 图(5)M=10,N=50的解析解的图像

图(6)M=10,N=50的数值解的图像 图(7)M=10,N=50的绝对误差的图像 此时，我们观察绝对误差发现，扩散方程(1.3.1)时不收敛不稳定的。而前面我们已经知道，到网格比为12r时，方程是收敛稳定的。所以，我们可以验证，当网比12r时稳定，当12r时不稳定。 [参考文献] [1]李荣华，刘播.微分方程数值解法[M].北京.高等教育出版社.. [2]王曰朋.偏微分方程数值解[OL]. 未知.偏微分方程的Matlab解法[OL]. 周品，何正风.MATLAB数值分析.[M].北京.机械工业出版社..

附录： L=1; M=40; N=1600; alfa=1; lambda=;%网格比 %**********************************************% h=L/M;%空间步长 x=0:h:L; x=x'; tao=lambda*h^2/alfa;%时间步长 tm=N*tao;%热传导的总时间tm %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，与杆长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; end end length(U); %***************设置立体网格****************% for i=1:N+1 X(:,i)=x; end

for j=1:M+1 Y(j,:)=t; end mesh(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)); end end %mesh(x,t,z') legend('解析解'); xlabel('X'); ylabel('T'); zlabel('Z'); for j=1:M+1 for k=1:N+1