matlab圆柱内导热分离变量法,一维热传导方程数值解法及matlab实现分离变量法和有限差分法.doc...

一维热传导方程的Matlab解法分离变量法和有限差分法

问题描述

实验原理

分离变量法实验原理

有限差分法

实验目的

利用分离变量法和有限差分法解热传导方程问题 利用matlab进行建模构建图形 研究不同的情况下采用何种方法 从更深层次上理解热量分布与时间、空间分布关系。

模拟与仿真作业

分离变量法(代码)：

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(i*x).*exp(-i^2*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 纵坐标为t

s=(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);

end

end

disp(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(-i^2*t));

end;

surf(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+(200*(1-(-1)^i))/(i*pi)*(sin(i*x).*exp(-i^2*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/(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(-i^2*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=