一维热传导方程的前向 紧差分格式

1、中南林业科技大学本科课程论文学 院： 理学院 专业年级： 09信息与计算科学一班 课 程： 偏微分方程数值解法 论文题目：一维热传导方程的前向Euler和紧差分格式指导教师： 陈红斌 2012年7月学生姓名： 唐 黎 学 号： 20093936 分 工： 程序编写，数值例子 学生姓名： 何雄飞 学 号： 20093925 分 工： 格式建立，资料收集 学生姓名： 汪 霄 学 号： 20093938 分 工： 文档编辑，资料整理 学生姓名： 毛博伟 学 号： 20093931 分 工： 公式编辑，查找资料 学生姓名： 倪新东 学 号： 20093932 分 工： 数据分析，查找资料 学生姓名：

2、何凯明 学 号： 20093924 分 工： 数据分析，查找资料 目录1引言 .12物理背景.1 3网格剖分 .2 41.1向前格式建立 .24.1.2差分格式的求解 .44.1.3收敛性与稳定性.4 4.1.4 数值例子 .7 4.2.1紧差分格式建立 .10 4.2.2差分格式求解 .124.2.3数值例子 .13总结.17 参考文献 .18 附录 .191 引言本文考虑的一维非齐次热传导方程的定解问题：其中为正常数，为已知函数，目前常用的求解热传导方程的差分格式有前向Euler差分格式、向后Euler差分格式、Crank-Nicolson格式、Richardson格式本文将给出前向Eul

3、er格式和紧差分格式，并给出其截断误差和数值例子2 物理背景热传导是由于物体内部温度分布不均匀,热量要从物体内温度较高的点流向温度较低的点处.以函数表示物体在时刻,处的温度,并假设关于具有二阶连续偏导数,关于具有一阶连续偏导数.是物体在处的热传导系数,取正值.设物体的比热容为,密度为.根据热传导定律,热量守恒定律以及公式得如果物体是均匀的,此时以及均为常数.令,上式方程化为若考虑物体内有热源,其热源密度函数为,则有热源的热传导方程为其中.3 网格剖分取空间步长和时间步长，其中都是正整数用两族平行直线和将矩形域分割成矩形网格，网格节点为记以表示网格内点集合，即位于开矩形的网点集合；表示所有位于闭

4、矩形的网点集合；是网格界点集合. 引进如下记号:，分别称为无穷范式（一直范式）2范数（范数，平均范数），差商的2范数（差商的范式）和范式.向前格式建立定义上的网格函数其中 在结点处考虑微分方程(3.1-1),有 (3.2)将 和 代入(3.2),得到 (3.3-1)注意到初边值条件(3.1-2)和(3.1-3),有 (3.3-2) (3.3-3)在(3.3-1)(3.3-3)中略去小量项 (3.4)并用代替，得到如下差分格式 (3.5-1) (3.5-2) (3.5-3)称为差分格式的局部截断误差。记 (3.6-1)则有 (3.6-2)4.1.2 差分格式的求解记，称为步长比。差分格式(3.5

5、-1)可写为 上式表明第k+1层上的值由第k层上的值显示表示出来。若已知第k层的值，则由上式就可直接得到第k+1层上的值。有时也称(3.5-1)为古典显格式。可把古典显格式写成矩阵形式4.1.3 收敛性与稳定性收敛性 设为定解问题（3.1-1）（3.1-3）的解，为差分格式（）（）的解，则当时，有，其中由（3.6-1）定义证明记将（3.3-1）（3.6-1）分别于（3.5-1）（3.5-3）相减，得到误差方程当时，有证明完毕.稳定性 如果在应用差分格式(3.5-1)(3.5-3)时，计算右端函数有误差，计算初值有误差，则实际得到的是如下差分方程的解。 (3.8) 令 将(3.5-1)(3.5-

6、3)与(3.8)相减，可得摄动方程组 (3.9) 当时，有 (3.10)上式说明当和很小时，误差也很小。 摄动方程组(3.9)和差分方程(3.5-1)(3.5-3)的形式完全一样。上述结果可叙述如下。 当时，差分格式(3.5-1)(3.5-3)关于初值和右端项在下述意义下是稳定的：设为差分方程组 的解，则有 下面来考虑的情况。此时必存在当时， 于是 设 则（3.9）为 可以验证其解为 由此易知 由于当时， 所以不论初始误差多么小，均会使解有较大的误差。我们得到如下结论。定理3.1.4 当时，差分格式(3.5-1)(3.5-3)是不稳定的。由定理3.1.3和定理3.1.4知，当步长比时，差分格式

7、(3.5-1)(3.5-3)是不稳定的；当步长比时，差分格式(3.5-1)(3.5-3)是不稳定的。我们把这种稳定性称为条件稳定性。实际计算时选取步长比必须满足，即4.1.4 数值例子应用向前Euler 格式(3.5-1)-(3.5-3)计算定解问题上述定解问题的精确解为.部分节点处数值解、精确解和误差的绝对值数值解精确解|精确解-数值解|20(0.5,0.1)1.8219e+0001.8221e+0002.3008e-00440(0.5,0.2) 2.0134e+000 2.0138e+0003.4361e-00460(0.5,0.3) 2.2251e+000 2.2255e+000 4.1

8、249e-00480(0.5,0.4) 2.4591e+000 2.4596e+000 4.6788e-004100(0.5,0.5) 2.7178e+000 2.7183e+000 5.2148e-004120(0.5,0.6) 3.0036e+000 3.0042e+000 5.7794e-004140(0.5,0.7) 3.3195e+000 3.3201e+000 6.3932e-004160(0.5,0.8) 3.6686e+000 3.6693e+000 7.0677e-004180(0.5,0.9) 4.0544e+000 4.0552e+000 7.8118e-004200(0

9、.5,0.10) 4.4808e+000 4.4817e+000 8.6337e-004取不同不长时数值解的最大误差1/101/2008.6337e-004*1/201/8002.1748e-0043.9699e+0001/301/32005.4366e-0054.0003e+0001/401/128001.3591e-0054.0001e+000的误差曲面图（）的误差曲面图（）的误差曲线图的误差曲线紧差分格式建立令 则有定义网格函数 由Taylor展开，有 其中将上式中上标为k和k+1的两个等式相加并除以2，可得利用(3.42),得到 其中则有注意到初值条件，忽略小项，并用代替，得到如下差分

10、格式4.2.2 差分格式求解差分格式可写成将其写出矩阵形式其中在每一时间层上只要解一个三对角线性方程组. 数值例子应用紧差分格式(3.47-1)-(3.47-3)计算定解问题上述定解问题的精确解为.最大误差 从表中可以看出，当空间步长缩小到原来的1/2，时间步长缩小到原来的1/4时，最大误差约缩小到原来的1/16.部分节点处数值解、精确解和误差的绝对值（）数值解精确解|精确解-数值解|10(0.5,0.1)1.8221e+0001.8221e+0001.0836e-00620(0.5,0.2) 2.0138e+000 2.0138e+000 1.6247e-00630(0.5,0.3) 2.2

11、255e+000 2.2255e+000 1.9547e-00640(0.5,0.4) 2.4596e+000 2.4596e+000 2.2195e-00650(0.5,0.5) 2.7183e+000 2.7183e+000 2.4750e-00660(0.5,0.6) 3.0042e+000 3.0042e+000 2.7435e-00670(0.5,0.7) 3.3201e+000 3.3201e+000 3.0351e-00680(0.5,0.8) 3.6693e+000 3.6693e+000 3.3554e-00690(0.5,0.9) 4.0552e+000 4.0552e+0

12、00 3.7087e-006100(0.5,1.0) 4.4817e+000 4.4817e+000 4.0990e-006取不同步长的数值解的最大误差1/101/1004.0990e-006*1/201/4002.5823e-0071.5873e+0011/401/16001.6139e-0081.6000e+0011/8001/48001.0090e-0091.5995e+0011/1601/192006.5710e-0111.5355e+001的误差曲面图的误差曲面图的误差曲线图的误差曲线图的误差曲线图总结：本文采用差分格式和紧差分格式来求解抛物型方程.差分格式采用二层共三个点,条件稳定

13、显格式,当时误差随着无限增长其稳定条件为网比,因此我们给出了当时的最大误差.而紧差分格式却采用二层共六个点,无条件隐格式,在每一时间层上均只需解三对角方程组.参 考 文 献1 孙志忠偏微分方程数值解法北京：科学出版社，2005.2 李荣华偏微分方程数值解法北京：高等教育出版社，2005.3 附录A程序代码流程图：clear;n=input('n=');%空间剖分数m=input('m=');%时间剖分数x=(0:1/n:1);t=(0:1/m:1);r=n*n/m;%网比l=2;A=diag(ones(1,n-1)*(1-2*r)+diag(ones(1,n-2

14、)*r,-1)+diag(ones(1,n-2)*r,1);%-精确解求解-for i=1:n+1 for j=1:m+1u(i,j)=exp(x(i)+t(j); endend%-%-初值条件求解-for i=1:n+1 u1(i,1)=exp(x(i);endfor j=1:m+1 u1(1,j)=exp(t(j);endfor j=1:m+1 u1(n+1,j)=exp(x(n+1)+t(j);end%-数值解求解-for j=1:m u1(2:n,j+1)=A*u1(2:n,j)+r*u1(1,j),zeros(1,n-3),r*u1(n+1,j)'end%-%for i=1:

15、10% M(i)=(u(6,i*20+1);% N(i)=(u1(6,i*20+1);%end%abs(M-N)'%-误差图-u-u1;x,t=meshgrid(x,t);surf(x,t,u'-u1');grid on;xlabel('x');ylabel('t');zlabel('|u(x,t)-u1(x,t)|');%-误差曲线-%E=abs(u(:,2)-u1(:,2);%plot(x,E,'-o');%u(:,m+1)'-u1(:,m+1)'%plot(x,u(:,m+1)'

16、;-u1(:,m+1)','-o');%xlabel('x');%ylabel('|u(x,1)-u1(x,1)|');%-最大误差求解-%for i=2:n% for j=2:m+1% K(i)=max(u(i,j)-u1(i,j);% end%end%K=max(K)%-附录B紧差分格式程序代码clear;n=input('n=');%空间剖分m=input('m=');%时间剖分H=1/n;%空间步长T=1/m;%时间步长x=(0:H:1);t=(0:T:1);r=n*n/m;%网比l=2;A=dia

17、g(ones(1,n-1)*(5/6+r)+diag(ones(1,n-2)*(1/12-r/2),-1)+diag(ones(1,n-2)*(1/12-r/2),1);%k+1层的系数矩阵B=diag(ones(1,n-1)*(5/6-r)+diag(ones(1,n-2)*(1/12+r/2),-1)+diag(ones(1,n-2)*(1/12+r/2),1);%k层的系数矩阵%-精确解求解-for i=1:n+1 for j=1:m+1u(i,j)=exp(x(i)+t(j); endend%-%-初值条件求解-for i=1:n+1 u1(i,1)=exp(x(i);%条件u(x,0

18、)=exp(x)endfor j=1:m+1 u1(1,j)=exp(t(j);%条件u(0,t)=exp(t)endfor j=1:m+1 u1(n+1,j)=exp(x(n+1)+t(j);%条件u(1,t)=exp(1+t)end%-数值解的计算-for j=1:m u1(2:n,j+1)=A(B*u1(2:n,j)+(1/12+r/2)*u(1,j)-(1/12-r/2)*u(1,j+1),zeros(1,n-3),(1/12+r/2)*u(n+1,j)-(1/12-r/2)*u(n+1,j+1)');end%-%for i=1:n% N(i)=(u(6,i*10+1);% M(i)=(u1(6,i*10+1);%end%abs(N-M)'%-误差曲面-%u-u1;%t,x=meshgrid(t,x);%surf(t,x,u1-u);%grid on;%xlabel('x');%ylabel('t');%zlabel('|u(x,t)-u1