一维热传导方程数值解法及matlab实现分离变量法和有限差分法.docx
《一维热传导方程数值解法及matlab实现分离变量法和有限差分法.docx》由会员分享，可在线阅读，更多相关《一维热传导方程数值解法及matlab实现分离变量法和有限差分法.docx（34页珍藏版）》请在冰点文库上搜索。
一维热传导方程数值解法及matlab实现分离变量法和有限差分法
一维热传导方程的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可计算的最大数相当于无穷
fori=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（[0pi010100]）;
所得到的三维热传导图形为：
有限差分法：
u=zeros（10,25）;%t=1x=pi构造一个1025列的矩阵（初始化为0）用于存放时间t和变量x横坐标为x纵坐标为t
s=（1/25）/（pi/10）^2;
fprintf（'稳定性系数S为:
\n'）;
disp（s）;
fori=2:
9
u（i,1）=100;
end;
forj=1:
25
u（1,j）=0;
u（10,j）=0;
end;
forj=1:
24
fori=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;
fori=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（[0pi010100]）;
所得到的热传导图形为：
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可计算的最大数，相当于无穷
fori=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（[0pi010100]）;
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;
fori=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（[0pi010100]）;
disp（s（:
6））;
当温度低于50度是时间为t=23.5*0.04=0.94
第三种情况（有限差分法）
在原来程序代码的基础上加上disp（u（5,:
））;可得出第五行（即x=pi/2）处温度随时间的变化情况
u=zeros（10,25）;%t=1x=pi10行25列横坐标为x纵坐标为t
s=（1/25）/（pi/10）^2;
fprintf（'稳定性系数S为:
\n'）;
disp（s）;
fori=2:
9
u（i,1）=100;
end;
forj=1:
25
u（1,j）=0;
u（10,j）=0;
end;
forj=1:
24
fori=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（'有限差分法解'）;
disp（u（5,:
））;
得到如下结果
我们知19列为50.350520列是数据为47.8902所以时间t为20*0.04=0.78
结论：
比较一二三种情况，我们得到不同的时间，这是由于：
1、加和不同一种为100，一种为无穷；
2、采用的方法不同：
一种为分离变量法，一种为有限差分法
造成的。
第一题完
解：
根据课本知识：
我们知道稳定性系数S是衡量由差分格式所得的数据是否稳定的因数，如果S>1/2，我们得到的数据就会不稳定，所画出的图形失真也就会很大。
一下就S<1/2和S>1/2，两种情况进行讨论（以第一题为例进行讨论）
当S<1/2时，取步长k=1/25，h=pi/10,S=0.4053;
u=zeros（10,25）;%t=1x=pi10行25列横坐标为x纵坐标为t
s=（1/25）/（pi/10）^2;
fprintf（'稳定性系数S为:
\n'）;
disp（s）;
fori=2:
9
u（i,1）=100;
end;
forj=1:
25
u（1,j）=0;
u（10,j）=0;
end;
forj=1:
24
fori=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（'有限差分法解s<1/2时'）;
获得的数据为
所得到的图形为
当S>1/2,时取步长k=1/25,h=pi/20,S=1.6211
u=zeros（20,25）;%t=1x=pi20行25列横坐标为x纵坐标为t
s=（1/25）/（pi/20）^2;
fprintf（'稳定性系数S为:
\n'）;
disp（s）;
fori=2:
19
u（i,1）=100;
end;
forj=1:
25
u（1,j）=0;
u（10,j）=0;
end;
forj=1:
24
fori=2:
19
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:
20）;
surf（x,t,u）;
xlabel（'t'）,ylabel（'x'）,zlabel（'T'）;
title（'有限差分法解s>1/2时'）;
所得数据为
所得到的图形为
显而易见，我们得到的图形与实际的热传导图形相差太大，这是由于S>1/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可计算的最大数，相当于无穷
fori=0:
m
u=u+8*（-1）^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=1x=pi20行100列横坐标为x纵坐标为t
s=（1/100）/（pi/20）^2;
fprintf（'稳定性系数S为:
\n'）;
disp（s）;
fori=1:
20
u（i,1）=i/20*pi;;
end;
forj=1:
100
u（1,j）=0;
end
forj=1:
99
fori=2:
19
u（i,j+1）=s*u（i+1,j）+（1-2*s）*u（i,j）+s*u（i-1,j）;
end
end
forj=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（'有限差分法解'）;
我们得到如图所示的热传导方程：
结论：
比较可得由以上两种方法作出的三维图形基本相同，符合热传导的热量分布随时间和空间的变化规律
第四题完成
升级会员 