应用物理软件训练
前言
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