热方程

热方程（heat equation）描述温度因热传导而变化的规律。均匀细杆的一维无热源模型为 $u_t=\kappa u_{xx}$，其中 $u(x,t)$ 是温度，$\kappa >0$ 是热扩散率。它把固定位置的升降温速度，与附近温度剖面的弯曲程度联系起来。

中间热、两端冷的一根杆

考虑长 $L=1\mathrm{m}$ 的杆，侧面绝热，两端始终接触20°C的恒温装置。初始温度选为 $$u(x,0)=20+12\sin \frac{\pi x}{L}({}^{\circ }\mathrm{C}).$$ 中点最热，为32°C，两端正好是20°C。取热扩散率 $\kappa =0.01{\mathrm{m}}^{\mathrm{2}}/\mathrm{s}$。这是便于复算的合成案例，不把该数值对应到某一种实际材料。

我们希望回答：中点多久降到26°C？整根杆是否“每处都以同样速度降温”？这需要同时求各处温度的变化，单个平均温度的常微分方程不能保留全部空间信息。

从小段杆的能量收支开始

设杆截面积为 $A$，材料密度为 $\rho$，比热容为 $c_p$，三者均取常数。单位分别是 ${\mathrm{m}}^{\mathrm{2}}$、$\mathrm{k}\mathrm{g}/{\mathrm{m}}^{\mathrm{3}}$、$\mathrm{J}/\mathrm{(}\mathrm{k}\mathrm{g}\mathrm{K}\mathrm{)}$。温升一摄氏度与温升一开尔文大小相同，因此下面的温差、导数可以用任一种温标的增量计算。

令 $q(x,t)$ 是朝右为正的热流密度，单位 $\mathrm{W}/{\mathrm{m}}^{\mathrm{2}}$。区间 $[x,x+\Delta x]$ 中，热流从左进入、从右离开。若没有体热源，精确的积分收支为 $$\rho c_{p}A\frac{d}{dt}{\int }_x^{x+\Delta x}u(\xi ,t)d\xi =Aq(x,t)-Aq(x+\Delta x,t).$$

两支箭头采用同一个正方向；净流入要用左端值减右端值，而不是把两端相加。

除以 $A\Delta x$，令段长趋于零，在足够光滑的条件下得到 $$\rho c_p{u}_t=-q_x.$$ 负号有直接含义：若右端流出的热多于左端流入的热，杆段储能减少。

傅里叶导热定律在本模型中写为 $$q=-k{u}_x,$$ 其中导热系数 $k>0$ 的单位是 $\mathrm{W}/\mathrm{(}\mathrm{m}\mathrm{K}\mathrm{)}$。若右边温度较高，$u_x>0$，热流便向左，故需要负号。这是一条材料输运关系，与能量守恒共同构成模型。将它代入收支，并假设 $k$ 常数，得到 $$u_t=\kappa u_{xx},\kappa =\frac{k}{\rho c_p}.$$ 右侧单位为 $({\mathrm{m}}^{\mathrm{2}}/\mathrm{s})(\mathrm{K}/{\mathrm{m}}^{\mathrm{2}})=\mathrm{K}/\mathrm{s}$。若导热系数随位置变化，应保留 $\rho c_p{u}_t=(k{u}_x{)}_x$；若每单位体积每秒另生热 $Q$，还要加上 $Q$。MIT 18.03 讲义，Lecture 29：由能量收支推导热方程

为什么出现二阶空间导数

在光滑剖面中，中心差分 $$u_{xx}(x,t)\approx \frac{u(x-h,t)-2u(x,t)+u(x+h,t)}{h^2}$$ 比较的是中间点与两侧的平均值。若两侧平均比中间暖，分子为正，热方程预言中间升温；若中间形成向下弯的热峰，分子为负，热峰下降。

温度的一阶导数决定热流，热流的一阶空间变化又决定积累，因此温度出现二阶导数。均匀斜坡虽有热流，但各处流入与流出相等，内部温度可以暂时不变；“有热流”和“在升温”是不同判断。

把方程、初值和边界一起代入

令 $\theta =u-20$ 表示相对恒温端的温差。它满足 $${\theta }_t=\kappa {\theta }_{xx},\theta (0,t)=\theta (L,t)=0,\theta (x,0)=12\sin \frac{\pi x}{L}.$$ 初始剖面的形状是一个正弦拱。先尝试让它只改变高度：$\theta (x,t)=B(t)\sin (\pi x/L)$。于是 $${\theta }_t=B^{\prime }(t)\sin \frac{\pi x}{L},{\theta }_{xx}=-\frac{{\pi }^2}{L^2}B(t)\sin \frac{\pi x}{L}.$$ 在内部正弦不为零，约去它得到 $B^{\prime }=-(\kappa {\pi }^2/L^2)B$。由 $B(0)=12$ 解得 $$u(x,t)=20+12e^{-\kappa {\pi }^{2}t/L^2}\sin \frac{\pi x}{L}.$$ 它在两端恒为20，在零时刻恢复给定剖面，而且刚才已经逐项核验了 PDE。指数中的 $\kappa t/L^2$ 无量纲。

比较同一个位置在不同曲线上的高度。各点相对20°C的温差按相同比例衰减，但绝对降温量不相同。

中点的温差为 $12e^{-0.01{\pi }^{2}t}$。降到26°C即温差减半，因此 $$12e^{-0.01{\pi }^{2}t}=6,t=\frac{\log 2}{0.01{\pi }^2}\approx 7.0230\mathrm{s}.$$ 十秒时中点约为24.4725°C。四分之一处还要乘 $\sin (\pi /4)=1/\sqrt{2}$，所以十秒时约为23.1625°C；两处都趋向20°C，但温度值并不相同。

一种形状怎样推广到更多初始温度

前面用的是一个分离变量解：空间形状乘时间幅度。对齐次端点，假设 $\theta =X(x)B(t)$，代入后在乘积非零处可分离为 $$\frac{B^{\prime }}{\kappa B}=\frac{X^{″}}{X}=-\lambda .$$ 左边只随时间变、右边只随位置变，要在整个区域相等，必须是同一个常数。由 $X^{″}+\lambda X=0$ 和 $X(0)=X(L)=0$，非零解要求 $${\lambda }_n={\left(\frac{n\pi }{L}\right)}^2,X_n(x)=\sin \frac{n\pi x}{L},n=1,2,\dots .$$ 这里负的或零的 $\lambda$ 只产生满足两端零值的零解；正值时，从 $X(0)=0$ 去掉余弦项，再由 $\sin (\sqrt{\lambda }L)=0$ 选出上述数列。

每个形状的时间因子是 $e^{-\kappa (n\pi /L{)}^{2}t}$。例如合法初值 $20+10\sin (\pi x/L)+3\sin (2\pi x/L)$ 的解就是两种模式按各自指数衰减后相加。第二模式的衰减率是第一模式的四倍，细密起伏因此更快消失。更一般初值可用正弦级数展开；级数及其导数是否收敛须按初值的正则性讨论，不应仅凭形式相加就声称在所有端点都逐项可微。Lebl：分离变量与端点绝热

端点绝热时，消失的是温差，不是平均温度

现在换一个实验：两端也绝热，不再由恒温器带走能量。热流为零意味着 $u_x(0,t)=u_x(L,t)=0$，称为齐次 Neumann 边界条件；前面的指定端点温度称为 Dirichlet 边界条件。

取新的初值 $u(x,0)=30+6\cos (\pi x/L)$，解为 $$u(x,t)=30+6e^{-\kappa {\pi }^{2}t/L^2}\cos \frac{\pi x}{L}.$$ 其空间导数在两端为零。左端由36°C下降，右端由24°C上升，中点一直为30°C。

这幅图的边界已换成零热流。与恒温端图比较，不能把两种实验的最终温度混为一谈。

积分方程给出平均温度守恒的直接证明： $$\frac{d}{dt}{\int }_0^{L}udx=\kappa [u_x{]}_0^L=0.$$ 余弦项的积分为零，所以平均一直为30°C。前一个两端恒温的实验则把能量传给外界，平均温度会下降。是否守恒由边界通量决定，不能从“内部没有热源”单独推出。

为什么这个解唯一，怎样检查数值结果

以恒温端问题为例，设两个足够光滑的解具有相同初值和端点值，它们之差 $w$ 满足零初值、零边界的热方程。考虑温差的平方积分： $$\frac{d}{dt}\frac{1}{2}{\int }_0^L{w}^{2}dx=\kappa {\int }_0^{L}w{w}_{xx}dx=\kappa [w{w}_x{]}_0^L-\kappa {\int }_0^L{w}_x^{2}dx\le 0.$$ 边界项因为 $w=0$ 而消失。这个非负量起初为零且不能增加，只能一直为零，故 $w=0$。这里的平方积分用于衡量两个解的差，不是温度本身的物理内能。

数值上，把空间二阶导数用中心差分、时间导数用向前差分，可得到 $$U_j^{n+1}=r{U}_{j-1}^n+(1-2r)U_j^n+r{U}_{j+1}^n,r=\frac{\kappa \Delta t}{(\Delta x{)}^2}.$$ 当 $0\le r\le 1/2$ 时，新值是三个旧值的非负加权平均，不会凭空超过它们的最大值。对本例取 $\Delta x=0.25\mathrm{m}$、$\Delta t=1\mathrm{s}$，有 $r=0.16$；起初中间三个值约为28.4853、32、28.4853°C，中点下一步为 $$32+0.16(28.4853-64+28.4853)\approx 30.8753({}^{\circ }\mathrm{C}).$$ 解析解在一秒时中点约30.8722°C，可直接比较误差。若只把时间步长增到4秒，$r=0.64$，中间权重变成负数，前面的最大值论证失效；例如让两端误差为零、三个内部点的误差成比例于 $(1,-\sqrt{2},1)$；把它代入更新式，每个内部误差都会乘上 $1-2r-\sqrt{2}r\approx -1.1851$，绝对值大于一，反复更新会放大这个模式。连续热方程的平滑性质，并不能挽救一个不稳定的离散格式。

模型还省略了相变、温度依赖材料参数和辐射等机制。在适用的宏观尺度上，热方程解释温差扩散；它的无限空间解可在任意正时间出现遍及全域的非零尾部，这一数学性质不能解释为真实物质或信息可以任意快传播。

参考资料

• MIT，18.03 Lecture Notes，Lecture 29，局部能量收支、傅里叶定律与热方程。
• Jiří Lebl，Notes on Diffy Qs，§4.6，正弦模式、分离变量与绝热边界。
• Victor Ivrii，Partial Differential Equations，§3.2，热方程的最大值原理与唯一性背景。
• 先修：偏微分方程、积分、一阶线性微分方程；相关：波动方程、微分方程数值解、量纲分析。