应用物理软件训练
前言
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 程序流程图
四、 源程序
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
un=0;
for k=1:N
unn=2*(U0-u0)*(-1)^k.*sin(k.*pi.*(X.^2+Y.^2).^0.5).*... exp(-k^2*pi^2*a2*t(tt))./(pi.*(X.^2+Y.^2).^0.5);
un=unn+un;
end
mesh(X,Y,un);
axis([-1 1 -1 1 -0.4 0]);
pause(0.1)
end
figure(2)
wn=0;
for k=1:N
wnn=2*(U0-u0)*(-1)^k.*sin(k.*pi.*RR).*...
exp(-k^2*pi^2*a2*TT)./(pi*k.*RR);
wn=wnn+wn;
end
waterfall(RR,TT,wn)
xlabel('r')
ylabel('t')
五、程序调试情况
1、有限长细杆的热传导
开始时
一段时间后2、（1）非齐次方程的解析解
（2）非齐次方程的数值解（差分法）
3、二维热传导问题
4、三维热传导问题
解析解的动画图
解析解的瀑布图
六、仿真中遇到的问题
几乎所有的工程问题都能转化成数学模型来解，而且借助MATLAB，大多数的模型的数值解的精确度均能满足要求。但是，存在的问题也不少。首先，数值解法存在许多局限性，一个解只能适用于一个或几个模型，或者一个或几个方程。而解析解的得到能使我们得出所有同类问题的通解，并且精确度高于数值解。这是由于数学的发展程度还不足以满足自然科学的发展要求，数值解法只是一个权宜之计。其次，MATLAB虽然能处理大量的数学问题，但其命令繁多，再加上各种工具箱，要完全学会和很好的使用MATLAB不是一件容易的事情，在编辑和阅读程序时通常要借助工具书查询相关命令，这样就增加了使用难度，使得MATLAB 不能广泛的普及。再者，要合理的使用MATLAB来解决数学问题，必需是建立在良好的数学基础之上的，这就势必要求MATLAB的使用者有扎实的数学功底，这又给MATLAB的普及带来了挑战。最后，由于工程中的导热问题的数学模型并不一都能很顺利的建立，这就给使用MATLAB解决导热问题增加了难度。
七、结束语
在这短短的一周内从开始的一头雾水，到自己看书学习，到同学讨论，再进行整个题目的理论分析和计算，参考课程上的代码，写出自己的代码。
我们也明白了学无止尽的道理，在我们所查的很多参考书中，很多知识是我们从没有接触过的，我们对它的了解还仅限于皮毛，对它的很多功能以及函数还不是很了解，所以在这个学习的过程中我们穿越在知识的海洋中，一点一点吸取着它的知识。在MATLAB编程中需要很多的参考书，要尽量多的熟悉matlab自带的函数及其作用，因为matlab的自带函数特别多，基本上能够满足一般的数据和矩阵的计算，所以基本上不用你自己编函数。这一点对程序非常有帮助，可以使程序简单，运行效率高，可以节省很多时间。本次课设中用了很多MATLAB自带的函数，使程序变得很简单。
把基本的知识看过之后，就需要找一个实际的程序来动手编一下，不要等所有的知识都学好之后再去编程，你要在编程的过程中学习，程序需要什么知识再去补充，编程是一点一点积累的，所以你要需做一些随手笔记什么的。
在编写程序代码时，需要什么函数，需要什么模块就应该去着重看那个知识点，不要一步登天，一步一步学，如果太急于把所有东西都学到，也是不好的，
更是实现不了的。所以那时一天一天积累的，慢慢地学通这个软件。
八、参考文献
《数学物理方程的MATLAB解法与可视化》彭芳麟著清华大学出版社
《量子物理学中的常用算法与程序》井孝功赵永芳蒿凤有编著哈尔滨工业大学出版社
《计算物理基础》彭芳麟著高等教育出版社
§2热传导方程的初值问题 一维热传导方程的初值问题(或Cauchy 问题) ?? ???+∞<<∞-=>+∞<<∞-=??-??x x x u t x t x f x u a t u ),()0,(0 ,),,(2 2 2? （） 偏导数的多种记号xx x t u x u u x u u t u =??=??=??22,,. 问题也可记为 ?? ?+∞ <<∞-=>+∞<<∞-=-x x x u t x t x f u a u xx t ),()0,(0 ,,),(2?. Fourier 变换 我们将用Fourier 变换法求解热传导方程的柯西问题.为此我们将着重介绍Fourier 变换的基本知识.Fourier 变换在许多学科中是重要使用工具. 可积函数,设)(x f f =是定义在),(+∞-∞上的函数, 且对任意A B <，()f x 在[,]A B 上 可积，若积分 ? +∞ ∞ -dx x f )(收敛,则称)(x f 在),(+∞-∞上绝对可积。 将),(+∞-∞上绝对可积函数形成的集合记为),(1 +∞-∞L 或),(+∞-∞L ， 即{ } ∞<=+∞-∞=+∞-∞? +∞ ∞ -dx x f f L L )(| ),(),(1 ，称为可积函数空间. 连续函数空间: ),(+∞-∞上全体连续函数构成的集合,记为),(+∞-∞C , {}上连续在),(|),(+∞-∞=+∞-∞f f C ， {}上连续在),(,|),(1+∞-∞'=+∞-∞f f f C 。 定义 若),(+∞-∞∈L f ,那么积分 ),(?)(21 λπ λf dx e x f x i =? +∞ ∞ -- 有意义,称为Fourier 变换, )(? λf 称为)(x f 的Fourier 变式(或Fourier 变换的象). ? +∞ ∞ --= =dx e x f f Ff x i λπ λλ)(21)(?)( 定理 (Fourier 积分定理)若),(),(1 +∞-∞?+∞-∞∈C L f ,那么我们有
使用差分方法求解下面的热传导方程 2 (,)4(,) 0100.2t xx T x t T x t x t =<<<< 初值条件：2(,0)44T x x x =- 边值条件：(0,)0(1,)0 T t T t == 使用差分公式 1,,1,2 2 2 (,)2(,)(,) 2(,)()i j i j i j i j i j i j xx i j T x h t T x t T x h t T T T T x t O h h h -+--++-+= +≈ ,1,(,)(,) (,)()i j i j i j i j t i j T x t k T x t T T T x t O k k k ++--= +≈ 上面两式带入原热传导方程 ,1,1,,1,2 2i j i j i j i j i j T T T T T k h +-+--+= 令2 24k r h =，化简上式的 ,1,1,1,(12)()i j i j i j i j T r T r T T +-+=-++ i x j t j
编程MA TLAB 程序，运行结果如下 1 x t T function mypdesolution c=1; xspan=[0 1]; tspan=[0 0.2]; ngrid=[100 10]; f=@(x)4*x-4*x.^2; g1=@(t)0; g2=@(t)0; [T,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid); [x,t]=meshgrid(x,t); mesh(x,t,T); xlabel('x') ylabel('t') zlabel('T') function [U,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid) % 热传导方程：
前言 本文只是针对小白而写，可以使新手对热传导理论由很浅到不浅的认识，如想更深学习热传导知识，请转其它文档。 一、概念与常量 1、温度场： 指某一时刻下，物体内各点的温度分布状态。 在直角坐标系中：； 在柱坐标系中：； 在球坐标系中：。 补充：根据温度场表达式，可分析出导热过程是几维、稳态或非稳态的现象，温度场是几维的、稳态的或非稳态的。 2、等温面与等温线： 三维物体内同一时刻所有温度相同的点的集合称为等温面； 一个平面与三维物体等温面相交所得的的曲线线条即为平面温度场中的等温线。 3、温度梯度： 在具有连续温度场的物体内，过任意一点P温度变化率最大的方向位于等温线的法线方向上。称过点P的最大温度变化率为温度梯度(temperature gradient)。用grad t表示。 定义为： 补充：温度梯度表明了温度在空间上的最大变化率及其方向，是向量，其正向与热流方向恰好相反。对于连续可导的温度场同样存在连续的温度梯度场。
在直角坐标系中： 3、导热系数 定义式：单位 导热系数在数值上等于单位温度降度(即1)下，在垂直于热流密度的单位面积上所传导的热流量。导热系数是表征物质导热能力强弱的一个物性参数。 补充：由物质的种类、性质、温度、压力、密度以及湿度影响。 二、热量传递的三种基本方式 热量传递共有三种基本方式：热传导；热对流；热辐射 三、导热微分方程式（统一形式：） 直角坐标系： 圆柱坐标系： 球坐标系： 其中，称为热扩散系数，单位，为物质密度，为物体比热容，为物体导热系数，为热源的发热率密度，为物体与外界的对流交换系数。 补充： 1处研究的对象为各向同性的、连续的、有内热源、物性参数已知的导热物体。 2稳态温度场，即则有：，此式称为泊松方程。 3无内热源的稳态温度场，则有：，此式称为拉普拉斯方程。 四、单值条件 导热问题的单值条件通常包括以下四项： 1几何条件：表示导热物体的几何形状与大小（一维、二维或三维）
第一章 热传导方程 本章介绍最典型的抛物型方程—热传导方程,在研究热传导,扩散等物理现象时都会遇 到这类方程. §1 热传导方程及其定解问题的导出 1.1热传导方程的导出 物理模型 在三维空间中,考虑一均匀,各向同性的物体Ω,假定它内部有热源,并且与周围介质有热交换,需要来研究物体内部温度的分布和变化. 以函数),,,(t z y x u 表示物体Ω在位置),,(z y x 及时刻t 的温度.物体内部由于各部分温度不同,产生热量的传递,它们遵循能量守恒定律. 能量守恒定律 物体内部的热量的增加等于通过物体的边界流入的热量与由物体内部的热源所生成的热量的总和 . 在物体Ω内任意截取一块D .现在时段],[21t t 上对D 使用能量守恒定律. 设),,,(t z y x u u =是温度(度),c 是比热(焦耳∕度·千克),ρ是密度(千克/米3), q 是热流密度(焦耳/秒·米2),0f 是热源强度(焦耳/千克·秒). 注意到在dt 时段内通过D 的边界D ?上小块dS 进入区域D 的热量为dSdt n q ?-(n 是 D ?的外法向),从而由能量守恒律,我们有 ,)||(21 21 120??????????+?-=-?==t t D t t D D t t t t dxdydz f dt ds n q dt dxdydz u u c ρρ （1.1） 大家知道,热量流动的原因是因为在物体内部存在温差.依据传热学中的傅立叶实验定律,在一定条件下,热流向量与温度梯度成正比 ,u k q ?-= （梯度? ?? ? ????????==?z u y u x u gradu u ,,） （1.2） 这里负号表明热量是由高温向低温流动,k 是物体的导热系数.
应用物理软件训练 前言 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）有限长细杆的热传导定解问题
热传导方程傅里解
————————————————————————————————作者：————————————————————————————————日期：
热传导在三维的等方向均匀介质里的传播可用以下方程表达： 其中： ?u =u(t, x, y, z) 表温度，它是时间变量t 与空间变量(x,y,z) 的函数。 ?/是空间中一点的温度对时间的变化率。 ?, 与温度对三个空间座标轴的二次导数。 ?k决定于材料的热传导率、密度与热容。 热方程是傅里叶冷却律的一个推论（详见条目热传导）。 如果考虑的介质不是整个空间，则为了得到方程的唯一解，必须指定u 的边界条件。如果介质是整个空间，为了得到唯一性，必须假定解的增长速度有个指数型的上界，此假定吻合实验结果。 热方程的解具有将初始温度平滑化的特质，这代表热从高温处向低温处传播。一般而言，许多不同的初始状态会趋向同一个稳态（热平衡）。因此我们很难从现存的热分布反解初始状态，即使对极短的时间间隔也一样。 热方程也是抛物线偏微分方程最简单的例子。 利用拉普拉斯算子，热方程可推广为下述形式
其中的是对空间变量的拉普拉斯算子。 热方程支配热传导及其它扩散过程，诸如粒子扩散或神经细胞的动作电位。热方程也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与 Ornstein-Uhlenbeck 过程。热方程及其非线性的推广型式也被应用于影像分析。量子力学中的薛定谔方程虽然有类似热方程的数学式（但时间参数为纯虚数），本质却不是扩散问题，解的定性行为也完全不同。 就技术上来说，热方程违背狭义相对论，因为它的解表达了一个扰动可以在瞬间传播至空间各处。扰动在前方光锥外的影响通常可忽略不计，但是若要为热传导推出一个合理的速度，则须转而考虑一个双曲线型偏微分方程。 以傅里叶级数解热方程[编辑] 以下解法首先由约瑟夫·傅里叶在他于1822年出版的著作Théorie analytique de la chaleur（中译：解析热学）给出。先考虑只有一个空间变量的热方程，这可以当作棍子的热传导之模型。方程如下： 其中u = u(t, x) 是t和x的双变量函数。 ?x是空间变量，所以x∈[0,L]，其中L表示棍子长度。
例4 周期初始温度分布 求解热传导方程t xx u u =,(,0)x t -∞<<+∞>给定初始温度分布 (,0)1cos 2,()u x x x =+-∞<<+∞。 解 4(,)1cos2t u x t e x -=+. 初始高斯温度分布 例 5求解定解问题22 22 0,(,0) (,0),()kx u u a x t t x u x e x -???-=-∞<<+∞>?????=-∞<<+∞? ， 其中常数0k >. 解 22()4(,)()x s a t u x t s e ds ?-- +∞ -∞ = ? 22 2()4x s ks a t e e ds -- +∞ --∞ = ? 222 2(41)24ka t s xs x a t e ds +-+- +∞ -∞ = ? 222 22224(41)()41414x ka t ka t s x ka t ka t a t e ds +- +++- +∞ -∞ = ? 22 2 222(41)()41 441 k ka t x x s ka t a t ka t e e ds +---+∞ ++-∞ = ? 2241 k x ka t e - += 2241 k x ka t - += . §3初边值问题 设长度为l ,侧表面绝热的均匀细杆,初始温度与细杆两端的温度已知,则杆上的温度分布 ),(t x u 满足以下初边值问题 ?? ? ??≤<==≤≤=<<<<=-T t t g t l u t g t u l x x x u T t l x t x f u a u xx t 0),(),(),(),0(,0), ()0,(0,0),,(212? 对于这样的问题,可以用分离变量法来求解. 将边值齐次化
第八章 热传导方程的傅里叶解 第一节 热传导方程和扩散方程的建立 8.1.1 热传导方程的建立 推导热传导方程和前面弦振动所用的数学方法完全相用，不同之处在于具体的物理规律不同。这里用到的是热学方面的两个基本规律，即能量守恒和热传导的傅里叶实验定律。 热传导的傅里叶实验定律：设有一块连续的介质，选定一定的坐标系，并用(,,,)u x y z t 表示介质内空间坐标为的一点在t 时刻的温度。若沿x 方向有一定的温度差，在x 方向也就一定有热量的传递。从宏观上看，单位时间内通过垂直x 方向的单位面积的热量q 与温度的沿x 方向的空间变化率成正比，即 x u q k x ?=-? （8-1.1） q 称为热流密度，k 称为导热系数。公式中的负号表示热流的方向和温度变化的方向正好相 反，即热量由高温流向低温。 研究三维各向同性介质中的热传导，在介质中三个方向上存在温度差，则有 x u q k x ?=-?，y u q k y ?=-?，z u q k z ?=-? 或 q k u =-?r 即热流密度矢量q r 与温度梯度u ?成正比。 下面以一维均匀细杆为例，根据傅里叶实验定律和能量守恒定律推导介质中的热传导方程。 第一步，定变量。研究介质x 位置处在t 时刻的温度(,)u x t 。 第二步，取局部。在介质内部隔离出从x 到x x +?一段微元长度，在t 到t t +?时间内温度的变化(,)(,)u u x t t u x t ?=+?-。 第三步，立假设。假设均匀介质的横截面积为A ，质量密度为ρ，比热为c ，热传导系数为k 。
第四步，找规律。隔离出来的微元长度在t 到t t +?时间内吸收的热量为： Q c m u c A x u ρ=????=???? （8-1.2） 在t 到t t +?时间内，同过x 位置处的横截面的热量为： 1x x x Q q A t k u A t =???=-?? （8-1.3） 在t 到t t +?时间内，同过x x +?位置处的横截面的热量为： 2x x x x x Q q A t k u A t +?+?=???=-?? （8-1.4） 如果在微元段内有其他的热源，假设在单位时间单位体积内产生的热量为(,)F x t ，则该热源在微元内产生的热量为： (,)3Q F x t t A x =???? （8-1.5） 第五步，列方程。根据能量守恒定律，净流入的热量应该等于介质在此时间内温度升高所需要的热量。 123Q Q Q Q =-+ 即 (,)x x x x x c A x u k u A t k u A t F x t t A x ρ+?????=-??+??+???? (,)x x x x x u u u c k F x t t x ρ+?-??? =+?? 得到： (,)t xx k F x t u u c c ρρ = + 令 a = (,)(,)F x t f x t c ρ= 则得到热传导方程为 (,)2t xx u a u f x t =+ （8-1.6） 当介质内部无其他热源时，热传导方程是齐次的，为 2t xx u a u = （8-1.7） 8.1.2 扩散方程的建立