一维热传导方程数值解法及matlab实现分离变量法和有限差分法

问题描述

实验原理

分离变量法实验原理

有限差分法

实验目的

利用分离变量法和有限差分法解热传导方程问题利用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(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=zeros(10,25); %t=1 x=pi 10行25列横坐标为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