微分方程之液体冷却

在学习微分方程的时候，液体冷却是一个经常使用的例子。利用微分方程对液体冷却的过程进行建模，我们可以求解出液体温度随时间变换的函数。
本文利用一组实际测量的温度数据拟合、验证液体温度和时间之间的关系。

一、模型建立
假设环境温度为$T_e$，时刻$t$液体的温度是$T(t)$，$k$是比例系数，$k>0$。
那么，它们之间的关系可以表示为下面的公式：

$$\frac{dT}{dt}=-k(T-T_e)$$
直观地理解：公式左边 $\frac{dT}{dt}$表示温度降低的速度。如果液体的温度越高于环境温度，即 $(T-T_e)$的数值越大，那么温度下降越快。反之，如果液体的温度与环境温度越接近，温度下降地越慢，最终当液体温度和环境温度一样，即 $T-T_e=0$时，液体的温度不再变化。

二、模型求解
此处的方程比较简单，可以直接推导出解析解。

$$\frac{dT}{dt}=-k(T-T_e)$$
$$\frac{1}{T-T_e}{dT}=-kdt$$
$$\int\frac{1}{T-T_e}{dT}=-k\int dt$$
$$ln \lvert T-T_e\rvert = -kt+C$$
$$\lvert T-T_e\rvert=e^{-kt+C}$$
因为我们这里只考虑液体冷却，即 $T-T_e\ge0$，所以
$$T(t)=e^{-kt+C}+T_e$$
其中 $C$是一个常数。

当$t=0$时， $T(0)=e^C+T_e$，即 $e^C=T(0)-T_e$。
因此，最终的表达式为

$$T(t)=(T(0)-T_e)e^{-kt}+T_e$$

注意：$T(0)$是测量的起始温度，$T_e$是环境温度，都是已知数。公式中只有$k$是变量，需要通过数据拟合来求解。

三、数据拟合
实际测量的数据（表格1）。每1分钟测量一次液体的温度，一共测量了30次，30分钟。

表格1

时间点（分钟）	温度（摄氏度）
0	31.5
1	31.2
2	30.7
3	30.4
4	30.1
5	29.9
6	29.6
7	29.4
8	29.2
9	29.0
10	28.8
11	28.7
12	28.5
13	28.3
14	28.2
15	28.1
16	28.0
17	27.9
18	27.8
19	27.8
20	27.7
21	27.6
22	27.6
23	27.5
24	27.4
25	27.4
26	27.4
27	27.3
28	27.3
29	27.3
30	27.2

温度与时间的散点图。

接着，利用Matlab自带的曲线拟合工具（Curve Fitting）对数据进行拟合。

Step 1 输入数据，运行，读入内存

clear

t = 0:1:30;
T_empirical = [31.5, 31.2, 30.7, 30.4, 30.1, 29.9, 29.6, 29.4, 29.2, 29,   28.8, ...
     28.7, 28.5, 28.3, 28.2, 28.1, 28  , 27.9, 27.8, 27.8 ,27.7, 27.6, ...
     27.6, 27.5, 27.4, 27.4, 27.4, 27.3, 27.3, 27.3, 27.2];

注意：此处要运行程序，将数据读入内存。

Step 2 打开Matlab自带的Curve Fitting工具

Step 3 选择数据

在红色的框中，X data和Y data分别表示横纵坐标的数据，这里我们选择$t$和$T_-empirical$。绿色框中表示用来拟合数据的模型。此时，图中所示的是”Polynomial”，”Degree: 1”，一阶多项式，即线性函数。从拟合图中我们也可以看出。

Step 4 自定义拟合函数

此时，我们需要做的是，改变绿色框中的拟合函数，选择“Custom Equation”（自定义公式）。

因为在冷却模型中，只有$k$是变量。根据实际测量数据，$T(0)=31.5$，$T_e=27$，$T(0)-T_e=31.5-27=4.5$。所以，带入数值的表达式为：

$$T(t)=4.5*e^{-kt}+27$$

我们根据上面的公式，修改绿色框中的表达式。

Step 5 显示结果和参数估计值

现在，我们就可以看到函数拟和的实际效果，及其计算出来的参数$k$的数值。蓝色框中是数据拟合的效果，红色框中是计算出来的$k$的数值。

四、总结
利用微分方程建立的液体冷却模型的确可以很好地拟合实际测量数据。数学建模对于我们了解物理现象背后的规律很有帮助。另外，Matlab提供了方便、全面的数据拟合工具，并不需要我们自己去编写函数来实现。