此文档下载后即可编辑偏微分方程数值解所在学院：数学与统计学院课题名称：抛物形扩散方程的有限差分法及数值实例学生姓名：向聘抛物形扩散方程的有限差分法及数值实例1.1抛物型扩散方程抛物型偏微分方程是一类重要的偏微分方程。

考虑一维热传导方程：22(),0u ua f x t T t x∂∂=+<≤∂∂(1.1.1)其中a 是常数，()f x 是给定的连续函数。

按照初边值条件的不同给法，可将(1.1.1)的定解分为两类：第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ϕ=0,, ∞<<∞-x(1.1.2)第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ϕ=0,, 0x l <<(1.1.3) 及边值条件()()0,,0==t l u t u ， T t ≤≤0(1.1.4)假定()x f 和()x ϕ在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。

1.2抛物线扩散方程的求解下面考虑如下热传导方程22()(0.)(,)0(,0)()u ua f x t x u t u L t u x x ϕ⎧∂∂=+⎪∂∂⎪⎪==⎨⎪=⎪⎪⎩(1.2.1)其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。

取Nlh =为空间步长，MT =τ为时间步长，其中N ，M 是自然数，用两族平行直线jhx x j ==, ()N j ,,1,0Λ=和k t t k τ==,()M k ,,1,0Λ=将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。

其中(),jkx t 表示网格节点；hG表示网格内点（位于开矩形G 中的网格节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。

k j u 表示定义在网点(),j k x t 处的待求近似解，N j ≤≤0，M k ≤≤0。

现在对方程进行差分近似：(一) 向前差分格式=-+τk jk j u u 11122(())k k kj j j j j j u u u af f f x h +--++=(1.2.2)()j j j x u ϕϕ==0，ku 0=kNu =0(1.2.3) 计算后得：111(12)k k k k j j j j ju ru r u ru f τ++-=+-++(1.2.4) 其中，2a r h τ=，1,,1,0-=N j Λ，1,,1,0-=M k Λ。

显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。

方程组如下：10001210110002321210003432310001121(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----⎧=+-++⎪=+-++⎪⎪=+-++⎨⎪⎪⎪=+-++⎩M(1.2.5)若记()TkN k k k u u u 121,,,-=Λu ，()()()()T N x x x 121,,,-=ϕϕϕϕΛ，()()()()T N x f x f x f 121,,,-=τττΛf则显格式(1.2.4)可写成向量形式10,0,1,,1k k k M φ+⎧=+=-⎨=⎩L u Au f u(1.2.6) 其中⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----=r r r r r r r rr r21002100210021ΛO MO O O M O ΛA 而对于向前差分格式，当网比12r ≤时稳定，当12r >时不稳定。

这就意味着给定空间步长h 以后，时间步长τ必须足够小，才能保证稳定。

1.3抛物型热传导方程数值算例对于(1.2.1)所描述的扩散方程，取1a =已知方程的精确解为2sin t u e x ππ-=：22(,0)sin()(0,)(1,)001,00.5u ut x u x x u t u t x t π⎧∂∂=⎪∂∂⎪⎪=⎨⎪==⎪≤≤≤≤⎪⎩(1.3.1)设空间步长1/h M =，时间步长为0.5/N τ=，网格比2/r g h =。

向前格式：11122,1,...,1,1,...,k k k k kj jj j j u u u u u j M k Nh τ++---+==-=边值条件：110010111220,,k kk k kk k u u u u u j u u h τ++----+===， 11111122,,k k k k k k k M MM M M M M u u u u u j M u u hτ+++-+---+===. 初值条件：(,0)sin ,0,1,...,j j u x x j Mπ==对时间和空间进行分割，令M=40，N=1600，通过Matlab 计算得到该方程的解析解，数值解以及相对误差如下：图(1)解析解的图像图(2)数值解的图像图(3)M=40,N=1000的相对误差的图像我们取部分精确解和数值解进行比较，结果如表(1)数值解精确解相对误差x t0.1 0.1 0.1151 0.1152 4⨯1.170010-0.2 0.2 0.0815 0.0816 4⨯1.658010-0.3 0.3 0.0418 0.0419 4⨯1.275210-0.4 0.4 0.0183 0.0184 57.445610-⨯0.5 0.45 0.0117 0.0118 5⨯5.375510-0.6 0.35 0.0300 0.0301 4⨯1.067410-0.7 0.25 0.0684 0.0686 4⨯1.741010-0.8 0.15 0.5084 0.5085 5⨯7.589710-0.9 0.05 0.3654 0.3654 5⨯1.740710-表(1)数值解与精确解的比较由表(1)我们可以看出，精确解和数值解的绝对误差在410-以内，因此可以得出，在分割M=40，N=1600下，该有限差分方法对方程(1.3.1)是收敛和稳定的。

下面，我们比较在不同的分割下对有限差分算法精度的影响。

在扩散系数1a=不变的情况下，讲时间和空间进行更加细密的分割，取50,10000==，其中，M表示空间上的分割，N表M N示时间上的分割。

观察数值解与精确解在节点,x t处的绝对误差j k值，如下图所示：图(4)M=50,N=10000的相对误差的图像由图(3)和图(4)，两者在节点处的误差收敛分别是在410-10-和5以内，因此，可以得出的结论是：在收敛范围内，随着时间和空间的分割越细，节点数越多，精确解和解析解之间的绝对误差也越小，有限差分法的算法精度也越高。

最后，我们比较网比1/2r=时扩散方程的收敛情况。

r=以及1当网比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].北京.高等教育出版社.2009.1.[2]王曰朋.偏微分方程数值解[OL]. /link?url=-UAAXRI8_V547zVS6UwenT8rG KsbwNf9CuDmh2qmsy5K8eQ32PzhYazZ_sWiHfz_Pj3LA7ufHH4O 3t8NIlUCnXNUiyYhssqmNBeArsrLwQG[3]未知.偏微分方程的Matlab 解法[OL]./link?url=WBlR0q9n02YsNg6UOLkCQ4Z5 QOCefDKNdEglrNR2TJ8VlxaJ9IkCrLlaEDlJ_OHCLehf19UJU_B7P Re1skQTYaaLcAa1UWcwlv7GYsWNtG7[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';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的关系****************%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;endendlength(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');z=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');。