一维热传导方程数值解法及matlab实现.docx

1、一维热传导方程数值解法及matlab实现一维热传导方程的Matlab解法分离变量法和有限差分法2010/12/20问题描述实验原理分离变量法实验原理有限差分法实验目的利用分离变量法和有限差分法解热传导方程问题 利用matlab进行建模构建图形 研究不同的情况下采用何种方法 从更深层次上理解热量分布与时间、空间分布关系。模拟与仿真作业（1）分离变量法（代码）：x=0:0.1*pi:pi;y=0:0.04:1;x,t=meshgrid(x,y);s=0;m=length(j);%matlab可计算的最大数 相当于无穷for i=1:m s=s+(200*(1-(-1)i)/(i*pi)*(sin(

2、i*x).*exp(-i2*t); end;surf(x,t,s);xlabel(x),ylabel(t),zlabel(T);title( 分离变量法（无穷）);axis(0 pi 0 1 0 100);所得到的三维热传导图形为：有限差分法：u=zeros(10,25); %t=1 x=pi 构造一个1025列的矩阵（初始化为0）用于存放时间t和变量x 横坐标为x 纵坐标为ts=(1/25)/(pi/10)2;fprintf(稳定性系数S为:n);disp(s);for i=2:9 u(i,1)=100;end;for j=1:25 u(1,j)=0; u(10,j)=0;end;for j

3、=1:24 for i=2:9 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endenddisp(u);x,t=meshgrid(1:25,1:10);surf(x,t,u);xlabel(t),ylabel(x),zlabel(T);title( 有限差分法解);所得到的热传导图形为：（2）i分离变量法（取前100项和）x=0:0.1*pi:pi;y=0:0.04:1;x,t=meshgrid(x,y);s=0;for i=1:100 s=s+(200*(1-(-1)i)/(i*pi)*(sin(i*x).*exp(-i2*t); end;s

4、urf(x,t,u);xlabel(x),ylabel(t),zlabel(T);title( 分离变量法);axis(0 pi 0 1 0 100);所得到的热传导图形为：Ii有限差分法根据（1）我们有如下图结论：比较可得这两幅图基本相同，有限差分法和分离变量法对本题都适应（3）第一种情况（取无穷项）：在原来程序代码的基础上加上 disp(s(:,6); 可得出第六列（即x=pi/2）处温度随时间的变化情况x=0:0.1*pi:pi;y=0:0.04:1;x,t=meshgrid(x,y);s=0;m=length(j);%matlab可计算的最大数，相当于无穷for i=1:m s=s+(

5、200*(1-(-1)i)/(i*pi)*(sin(i*x).*exp(-i2*t); end;surf(x,t,s);xlabel(x),ylabel(t),zlabel(T);title( 分离变量法（无穷）);axis(0 pi 0 1 0 100);disp(s(:,6);我们得到如下一组数据： 当温度低于50度是 时间为 t=23.5*0.04=0.94第二种情况（取前100项和）在原来程序代码的基础上加上 disp(s(:,6); 可得出第六列（即x=pi/2）处温度随时间的变化情况x=0:0.1*pi:pi;y=0:0.04:1;x,t=meshgrid(x,y);r=0.04/

6、(0.1*pi)2;fprintf(稳定性系数S为：)disp(r);s=0;for i=1:100 s=s+(200*(1-(-1)i)/(i*pi)*(sin(i*x).*exp(-i2*t); end;surf(x,t,s);xlabel(x),ylabel(t),zlabel(T);title( 分离变量法);axis(0 pi 0 1 0 100);disp(s(:,6); 当温度低于50度是 时间为 t=23.5*0.04=0.94第三种情况（有限差分法）在原来程序代码的基础上加上 disp(u(5,:);可得出第五行（即x=pi/2）处温度随时间的变化情况u=zeros(10,2

7、5); %t=1 x=pi 10行25列 横坐标为x 纵坐标为ts=(1/25)/(pi/10)2;fprintf(稳定性系数S为:n);disp(s);for i=2:9 u(i,1)=100;end;for j=1:25 u(1,j)=0; u(10,j)=0;end;for j=1:24 for i=2:9 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endenddisp(u);x,t=meshgrid(1:25,1:10);surf(x,t,u);xlabel(t),ylabel(x),zlabel(T);title( 有限差分法解);

8、disp(u(5,:);得到如下结果我们知 19列为50.3505 20列是数据为 47.8902 所以时间t为 20*0.04=0.78结论：比较一二三种情况，我们得到不同的时间，这是由于：1、加和不同一种为100，一种为无穷；2、采用的方法不同：一种为分离变量法，一种为有限差分法 造成的。第一题完解：根据课本知识：我们知道 稳定性系数S是衡量由差分格式所得的数据是否稳定的因数，如果S1/2，我们得到的数据就会不稳定，所画出的图形失真也就会很大。一下就S1/2，两种情况进行讨论（以第一题为例进行讨论）当S1/2时，取步长k=1/25，h=pi/10, S= 0.4053;u=zeros(10

9、,25); %t=1 x=pi 10行25列 横坐标为x 纵坐标为ts=(1/25)/(pi/10)2;fprintf(稳定性系数S为:n);disp(s);for i=2:9 u(i,1)=100;end;for j=1:25 u(1,j)=0; u(10,j)=0;end;for j=1:24 for i=2:9 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endenddisp(u);x,t=meshgrid(1:25,1:10);surf(x,t,u);xlabel(t),ylabel(x),zlabel(T);title( 有限差分法解

10、s1/2,时 取步长k=1/25 ,h=pi/20, S=1.6211 u=zeros(20,25); %t=1 x=pi 20行25列 横坐标为x 纵坐标为ts=(1/25)/(pi/20)2;fprintf(稳定性系数S为:n);disp(s);for i=2:19 u(i,1)=100;end;for j=1:25 u(1,j)=0; u(10,j)=0;end;for j=1:24 for i=2:19 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endenddisp(u);x,t=meshgrid(1:25,1:20);surf(x,

11、t,u);xlabel(t),ylabel(x),zlabel(T);title( 有限差分法解s1/2时);所得数据为所得到的图形为显而易见，我们得到的图形与实际的热传导图形相差太大，这是由于S1/2，得到的数据是不稳定的，所以画出的图形也是不准确的。第三题完成根据课上所学知识，我们有如下方程：为便于解释做题，我们令：a=1 l=pi=x;下面开始求解：分离变量法根据课上所讲其中：我们有如下代码：x=0:0.1*pi:pi;y=0:0.4:10;x,t=meshgrid(x,y);u=0;m=length(j);%matlab可计算的最大数，相当于无穷for i=0:m u=u+8*(-1)

12、i/(pi*(2*i+1)2)*(sin(2*i+1)/2*x).*exp(-(2*i+1)2/4*t);end;surf(x,t,u);xlabel(x),ylabel(t),zlabel(T);title( 分离变量法（无穷）);disp(u);得到如图所示的热传导现象：有限差分法u=zeros(20,100); %t=1 x=pi 20行100列 横坐标为x 纵坐标为ts=(1/100)/(pi/20)2;fprintf(稳定性系数S为:n);disp(s);for i=1:20 u(i,1)=i/20*pi;end;for j=1:100 u(1,j)=0;endfor j=1:99 for i=2:19 u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j); endendfor j=1:100 u(20,j)=u(19,j);end;disp(u);x,t=meshgrid(1:100,1:20);surf(x,t,u);xlabel(t),ylabel(x),zlabel(T);title( 有限差分法解);我们得到如图所示的热传导方程：结论： 比较可得由以上两种方法作出的三维图形基本相同，符合热传导的热量分布随时间和空间的变化规律 第四题完成