在很多领域中往往需要对实验数据（或者调查数据）通过某种方程进行拟合，得到方程中的关键参数。在统计、物理、化学、生物等领域都是非常常见的方法。 非线性方程是指方程中存在非线性项，例如 e x e^x ex， ln ⁡ x \ln x lnx， x 3 x^3 x3， sin ⁡ x \sin x sinx等项时，求解过程通过人工计算很难得到准确的解。

Python的numpy库中的函数polyfit()能够拟合一元多项式 y = a 1 x m + a 2 x m − 1 + . . . + a m x + a m + 1 y = a_1 x^m + a_2 x^{m-1} +...+ a_m x + a_{m+1} y=a1​xm+a2​xm−1+...+am​x+am+1​numpy.polyfit()专门用于拟合这种形式的多项式，看下面的例子

图1-1多项式拟合 图1-2 多项式拟合方差随多项式次数的变化规律 简单举个例子， y = e x y=e^x y=ex这个函数也可以拟合成多项式的形式，随着多项式次数的增加，拟合结果方差逐渐减小。其原理是连续函数的泰勒公式，高数里证明了连续可导的函数总能表示成多项式的形式。 多项式拟合能在一定程度上获得变化规律，但是也存在着无法忽略的问题，超出原始数据范围的预测效果差，而且随着n的增加，有可能出现过拟合的情况（噪音信号被带入拟合结果）。

curve_fit()是scipy.optimize的函数，用于拟合已知形式的函数，接下来还是看个例子。仍然以 y = exp ⁡ ( x ) y=\exp(x) y=exp(x)为例说明，通过随机数函数增加生成一些噪音，使噪音范围为[-0.2, 0.2），然后分别用curve_fit()和polyfit()拟合并对比。用下面的函数拟合，a和n为待拟合参数。 y = a exp ⁡ ( n x ) y = a\exp(nx) y=aexp(nx) 代码如下：

图1-3 非线性拟合+多项式拟合过拟合 图1-4 非线性拟合 拟合的结果是a=0.9663855 , n=1.01433087，和预期的1, 1很接近，拟合效果还可以。多项式拟合当使用30次方时，就产生了过拟合的情况，偏离了实际的情况，所以多项式拟合还是要慎用，次数取小了误差大，取大了会过拟合。

对于光催化反应甲醛反应，其动力学方程可表示为 r = d c A d t = k 0 exp ⁡ ( − E R T ) c A m r = \frac{dc_A}{dt} = k_0\exp(-\frac{E}{RT})c_A^m r=dtdcA​​=k0​exp(−RTE​)cAm​ 分离变量后得 d c A c A m = k 0 exp ⁡ ( − E R T ) \frac{dc_A}{c_A^m} = k_0\exp(-\frac{E}{RT}) cAm​dcA​​=k0​exp(−RTE​) 根据m取值的不同，浓度与时间的变化规律可分为三种形式的函数 c A ( t ) = c 0 − k 0 exp ⁡ ( − E R T ) t ,   m = 0 ln ⁡ c A ( t ) = ln ⁡ c 0 − k 0 exp ⁡ ( − E R T ) t ,   m = 1 [ c A ( t ) ] 1 − m = c 0 1 − m − ( 1 − m ) k 0 exp ⁡ ( − E R T ) t ,   m ≠ 0 , 1 c_A(t) = c_0 - k_0\exp(-\frac{E}{RT})t,\ m=0\\ \ln c_A(t) = \ln c_0 - k_0\exp(-\frac{E}{RT})t,\ m=1\\ [c_A(t)]^{1-m} = c_0^{1-m} - (1-m)k_0\exp(-\frac{E}{RT})t,\ m\neq 0,1 cA​(t)=c0​−k0​exp(−RTE​)t, m=0lncA​(t)=lnc0​−k0​exp(−RTE​)t, m=1[cA​(t)]1−m=c01−m​−(1−m)k0​exp(−RTE​)t, m​=0,1 通过实验测得浓度随时间的变化规律cA(t)，根据不同的动力学模型对数据进行拟合。 图2-1 实验浓度变化

从实验数据初步判断m≠0，接下来将使用模型2和模型3对实验数据进行拟合。 ln ⁡ r = m ln ⁡ c 0 + ln ⁡ k 0 − E / R T \ln r = m\ln c_0 + \ln k_0 - E/RT lnr=mlnc0​+lnk0​−E/RT 反应速率的表达式更为统一，对反应速率方程取对数后得到的方程形式更简洁，因此为减少在拟合过程中不必要的误差，可使用速率方程进行拟合，计算过程见附录，拟合结果如下。 图2-2 速率与浓度变化规律 拟合得到m=1.75978，k0=24.8696，E=15949.7。拟合曲线与实验数据的对比如下图所示。 图2-3 浓度随时间变化规律 通过非线性回归方法计算得出甲醛光催化反应的反应级数为1.76，活化能为15.9 kJ/mol。

非理想反应器内流体存在返混现象，因此流体微元在反应器中的停留时间并不是完全相同的，而存在与反应器特性有关的概率分布，称为停留时间分布。停留时间分布对于认识反应器传质过程与评估反应器生产能力等都有重要参考价值。 停留时间分布的测试方法分为两类：阶跃法和脉冲法。阶跃法是指在某一时刻突然改变稳定运行的反应器入口的探针物质浓度，通过监测不同时间后出口处探针物质浓度的变化规律，直接得到离散的停留时间分布F(t)。脉冲法是指在某一瞬间，从入口注入一定量探针物质，通过监测出口浓度，直接得到停留时间分布密度E(t)。 图3-1 阶跃法停留时间分布示意图 停留时间分布函数F(t)与停留时间分布密度函数E(t)体现了反应器的某些工作特点，比如返混程度、平均停留时间等等。对于理想平推流反应器无返混，而理想全混釜反应器返混无穷大，通常非理想反应器的返混程度则是在两者之间。可见停留时间分布函数对于评估反应器生产能力、设计反应器体积等极为重要。 图3-2 停留时间分布函数

停留时间分布规律通常呈“S”型曲线，F(t)的最大值为1，因此可使用S型函数表示F(t)，其中b，c为拟合参数。将实验数据通过S函数拟合得到停留时间分布的连续函数。 F ( t ) = 1 1 + b exp ⁡ ( − c t ) F(t) = \frac{1}{1+b\exp(-ct)} F(t)=1+bexp(−ct)1​ 因此，进一步可解算停留时间分布密度函数E(t)，其形式可表示为如下方程： E ( t ) = d F ( t ) d t = b c exp ⁡ ( − c t ) ( 1 + b exp ⁡ ( − c t ) ) 2 E(t) = \frac{dF(t)}{dt} = \frac{bc\exp(-ct)}{(1+b\exp(-ct))^2} E(t)=dtdF(t)​=(1+bexp(−ct))2bcexp(−ct)​ 平均停留时间则可表示为如下： t ˉ = ∫ 0 ∞ t E ( t ) d t = ln ⁡ ( b + 1 ) c \bar t = \int_0^\infty tE(t) dt = \frac{\ln (b+1)}{c} tˉ=∫0∞​tE(t)dt=cln(b+1)​ 停留时间分布的方差表示为： σ t 2 = ∫ 0 ∞ ( t − t ˉ ) 2 E ( t ) d t = ∫ 0 ∞ t 2 E ( t ) d t − t ˉ 2 \sigma^2_t = \int_0^\infty (t-\bar t)^2E(t)dt = \int_0^\infty t^2E(t)dt - \bar t^2 σt2​=∫0∞​(t−tˉ)2E(t)dt=∫0∞​t2E(t)dt−tˉ2

表 1 阶跃法测定停留时间分布实验数据

（1）通过离散方法计算得到=50.58，σt2=341.54。 （2）通过拟合S型函数得到停留时间分布函数： F ( t ) = 1 1 + 63.5564 exp ⁡ ( − 0.0916 t ) F(t) = \frac{1}{1+63.5564\exp(-0.0916t)} F(t)=1+63.5564exp(−0.0916t)1​ 停留时间分布密度函数： E ( t ) = 5.8218 exp ⁡ ( − 0.0916 t ) ( 1 + 63.5564 exp ⁡ ( − 0.0916 t ) ) 2 E(t) = \frac{5.8218\exp(-0.0916t)}{(1+63.5564\exp(-0.0916t))^2} E(t)=(1+63.5564exp(−0.0916t))25.8218exp(−0.0916t)​ 平均停留时间： t ˉ = ln ⁡ ( 63.5564 + 1 ) 0.0916 = 45.50 \bar t = \frac{\ln (63.5564+1)}{0.0916} = 45.50 tˉ=0.0916ln(63.5564+1)​=45.50

停留时间分布方差： σ t 2 = 364.60 \sigma_t^2 = 364.60 σt2​=364.60 图3-3 拟合停留时间分布函数 *图3-4 拟合停留时间分布密度函数 拟合函数的计算结果与离散方法的计算结果相比，数值相差不超过10%，通过图像可以看到，S函数能较好的拟合出停留时间分布函数。经过多次尝试发现，S函数拟合返混程度较大的非理想反应器时效果较好，而对接近平推流的反应器拟合效果差。