热传导方程的求解

应用物理软件训练

前言

MATLAB 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。

MATLAB是矩阵实验室（Matrix Laboratory）的简称，和Mathematica、Maple 并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其

他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。

本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。

题目：热传导方程的求解

目录

一、参数说明 (1)

二、基本原理 (1)

三、MATLAB程序流程图 (3)

四、源程序 (3)

五、程序调试情况 (6)

六、仿真中遇到的问题 (9)

七、结束语 (9)

八、参考文献 (10)

一、参数说明

U=zeros(21,101) 返回一个21*101的零矩阵

x=linspace(0,1,100);将变量设成列向量

meshz(u)绘制矩阵打的三维图

axis([0 21 0 1]);横坐标从0到21，纵坐标从0到1

eps是MATLAB默认的最小浮点数精度

[X,Y]=pol2cart(R,TH);效果和上一句相同

waterfall(RR,TT,wn)瀑布图

二、基本原理

1、一维热传导问题

（1）无限长细杆的热传导定解问题

利用傅里叶变换求得问题的解是：

取得初始温度分布如下

这是在区间0到1之间的高度为1的一个矩形脉冲，于是得

（2）有限长细杆的热传导定解问题

其中20x 0≤≤，即L=20，取a=10且

得的解是

（3）非齐次方程定解问题是

解析解是

其中

2、二维热传导问题 定解问题

Ut=k^2(Uxx+Uyy) (b y a ≤≤≤≤0,

x 0) U(x=0,y,t)=0, u(x=a,y,t)=b

3sin

y

πμ U （x,y=0,t ）=0, u(x,y=b,t)=a

x

x ππμcos a 3sin

U （x,y,t=0）=0

3、三维热传导问题

球体内的热传导

令u=w+Uo,则w 的定解问题是 Wt=w ∆w W （r=ro ）=0 W(t=to)=uo-Uo

解为

ro

r

n e

n

r

uo Uo w o

r t a n n

n

πππsin

)1()

(22

222/1

-∞

=∑--=

r 为空间变量，并用x ，y 表示。

三、 MATLAB 程序流程图 开始

初始化定义

预设矩阵

初始条件

用for 语言

绘制动态图

四、 源程序

1、一维有限长细杆的热传导

x=0:20;t=0:0.01:1;a2=10; r=a2*0.01; u=zeros(21,101);

u(10:11,1)=1; 是把上述矩阵中的第10行，11行的第一列全部设成1

for j=1:100

u(2:20,j+1)=(1-2*r)*u(2:20,j)+r*(u(1:19,j)+u(3:21,j));

plot(u(:,j));

axis([0 21 0 1]);横坐标0到21，纵坐标0到1

pause(0.1)暂停0.1秒

end

meshz(u)

2、非齐次方程的定解问题

a2=50;b=5;L=1;

[x,t]=meshgrid(0:0.01:1,0:0.000001:0.0005);

Anfun=inline('2/L*(x-L/2).^2.*exp(-b*x/2/a2).*sin(n*pi*x/L)','x ','n','L','b','a2');%定义内联函数

u=0;

for n=1:30

An=quad(Anfun,0,1,[],[],n,L,b,a2);%inline函数中定义x为向量，其它为标量

un=An*exp(-(n*n*pi*pi*a2/L/L+b*b/4/a2/a2).*t).*exp(b/2/a2.*x).* sin(n*pi*x/L);

u=u+un;

size(u);

mesh(x,t,u);%x,t，u都为501行101列的矩阵

figure

subplot(2,1,1)

plot(u(1,:))

subplot(2,1,2)

plot(u(end,:))

end

差分法

dx=0.01;dt=0.000001;a2=50;b=5;c=a2*dt/dx/dx;

x=linspace(0,1,100);%将变量设成列向量

uu(1:100,1)=(x-0.5).^2;%初温度为零

figure

subplot(1,2,1)%初始状态

plot(x,uu(:,1),'linewidth',1);

axis([0,1,0,0.25]);

subplot(1,2,2)%演化图

h=plot(x,uu(:,1),'linewidth',1);

set(h,'EraseMode','xor')

for j=2:200

uu(2:99,2)=(1-2*c)*uu(2:99,1)+c*(uu(1:98,1)+ uu(3:100,1))-... b*dt/dx*(uu(3:100,1)-uu(2:99,1));

uu(1,2)==0;uu(100,2)==0;%边界条件

uu(:,1)=uu(:,2);

uu(:,1)

set(h,'YData',uu(:,1));

drawnow;

pause(0.01)

end

三维热传导问题

U0=2; u0=0; a2=2; N=10;

r=eps:0.05:1; theta=linspace(0,2*pi,100);

t=0.1:0.001:0.2;

[RR,TT]=meshgrid(r,t);

figure(1)

[R,TH]=meshgrid(theta,r);

[X,Y]=pol2cart(R,TH);

for tt=1:100