数学建模

2. 初等模型

张瑞
中国科学技术大学数学科学学院

初等模型

模型特点：

  • 模型的因素可用数值来表示;
  • 因素之间的联系可用初等函数来表示;
  • 实际问题可转化为函数问题,如函数最值、方程求根等。

椅子摆放问题

例 1. 把椅子往不平的地面上一放,通常只有三只脚着地,放不稳。然而只需 稍挪动几次,就可以使四只脚同时着地,放稳了。用数学语言对此现象给以 描述,并用数学工具来验证此现象是否总是成立。

一定的假设：

  1. 椅子的4条腿一样长。脚与地面的接触看作一个点，4脚的连线是正方形。
  2. 地面高度是连续变化的。地面不出现间断（如台阶）
  3. 相对对于椅脚的长度，地面是平的。也就是，任何时候，椅子可以保证3条脚同时着地

模型： 椅子4个脚记为ABCD。ABCD的中心为O。以O为原点，OA为x轴。椅子变化后，绕O旋转了$\theta$角度。 记$f(\theta)$为A,C两脚到地面的距离和，$g(\theta)$为C,D两脚到地面的距离和。

  • $f(\theta)\geq0$, $g(\theta)\geq0$。它们都是$\theta$的连续函数
  • 根据假设，总有3个脚在地面，则$f(\theta)$与$g(\theta)$总有一个为$0$。即$f(\theta)g(\theta)=0$
  • 若最初时$g(0)=0$, $f(0)>0$。则旋转90度后，$f(\frac{\pi}2)=0$, $g(\frac{\pi}2)>0$

求解问题： $f(\theta)$, $g(\theta)$为$\theta$的连续函数，$f(\theta)g(\theta)=0$，且$f(\frac{\pi}2)=g(0)=0$。证明：存在$\theta_0$，使$g(\theta_0)=f(\theta_0)=0$

证明： 令$h(\theta)=f(\theta)-g(\theta)$，则$h(0)>0$, $h(\frac{\pi}2)<0$。因此，存在$\theta_0$，使得

\[h(\theta_0)=f(\theta_0)-g(\theta_0)=0 \]

又

\[f(\theta_0)g(\theta_0)=0 \]

则有

\[f(\theta_0)=g(\theta_0)=0 \]

利用比例进行建模

两个变量之间呈现比例关系

\[y=kx, y=k x^2, y=k \ln x \]

延伸一下，

\[y=y_0+k x, y=a +b x+c x^2, y=a(1-\frac{x}b)x \]

猜想有比例关系的量，然后获取数据，用最小二乘法确定数据，最后再检验模型

著名的比例性

  1. Hooke定律： $F=kS$， $F$是被拉长或压缩$S$距离的弦的恢复力
  2. Newton定律： $F=ma$, $a$是质量$m$的物理的加速度
  3. Ohm定律： $V=iR$, $i$是电压$V$通过电阻$R$的导线时引起的电流
  4. Einstein相对论： $E=Mc^2$,
  5. Kepler第三定律： $T=cR^{\frac32}$，$T$是周期，$R$是到行星到太阳的平均距离

例 2. (Kepler 第三定律) 得到的行星周期与行星到太阳的平均距离数据为：

行星 周期T 平均距离(万公里)R
水星 88.0 5791 (0.38 天文单位)
金星 224.7 10820(0.72个天文单位)
地球 365.3 14960 (1.00 天文单位)
火星 687.0 22794 (1.52 天文单位)
木星 4331.8 77833 (5.20 天文单位)
土星 10760.0 142940 (9.54 天文单位)
天王星 30684 287099 (19.218 天文单位)
海王星 60188.3 450400 (30.06 天文单位)

确定T与R的关系

Python 3.5.2 (default, Nov 12 2018, 13:43:14) 
[GCC 5.4.0 20160609] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> a=np.array([0.38, 0.72, 1.0, 1.52, 5.2, 9.54, 19.218, 30.06])
>>> b=np.array([88.0, 224.7, 365.3, 687.0,4331.8,10760.0,30684,60188.3])
>>> b/a
array([  231.57894737,   312.08333333,   365.3       ,   451.97368421,
         833.03846154,  1127.88259958,  1596.6281611 ,  2002.27212242])
>>> b/a**(1.5)
array([ 375.67065946,  367.79373549,  365.3       ,  366.59906683,
        365.31153154,  365.16547446,  364.20817565,  365.1981868 ])

mm2-ex1

Python中处理最小二乘法：可以使用 scipy库中的leastsq函数

参考： https://www.cnblogs.com/NanShan2016/p/5493429.html

最小二乘解

有了预想的比例关系，和一些相应的数据，如何确定这个比例系数？

理论上，所有的数据应该契合这个比例关系，但实际上，得到的数据总有误差。

定义 1.
最小二乘法：若有数据集$(x_i, y_i), i=1,\cdots,n$。找参数$\theta$，使得

\[\sqrt{\sum_{i=1}^n(y_i-f(x_i,\theta))^2} \]

达到最小。

线性最小二乘

若$f(x)=a_1 \phi_1(x)+a_2\phi_2(x)+\cdots+a_m\phi_m(x)$，则找$a_1$, $a_2$, $\cdots$, $a_n$，使

\[S=\sum_{i=1}^n[y_i-(a_1\phi_1(x_i)+\cdots+a_m\phi_m(x_i))]^2 \]

达到最小。可得，

\[\begin{pmatrix} (\phi_1, \phi_1) & (\phi_1, \phi_2) \cdots & (\phi_1, \phi_m) \\ (\phi_2, \phi_1) & (\phi_2, \phi_2) \cdots & (\phi_2, \phi_m) \\ \vdots & \vdots & & \vdots \\ (\phi_m, \phi_1) & (\phi_m, \phi_2) \cdots & (\phi_m, \phi_m) \\ \end{pmatrix} \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_m \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^n y_i\phi_1(x_i) \\ \sum_{i=1}^n y_i\phi_2(x_i) \\ \vdots \\ \sum_{i=1}^n y_i\phi_m(x_i) \end{pmatrix} \]

其中$\displaystyle (\phi_k, \phi_l)=\sum_{i=1}^n\phi_k(x_i)\phi_l(x_i)$

由$S$关于$a_k$达到最小，有$\frac{\partial S}{\partial a_k}=0$, $k=1,2,\cdots,m$。即

\[2\sum_{i=1}^n\{[y_i-(a_1\phi_1(x_i)+\cdots+a_m\phi_m(x_i))]\phi_k(x_i)\}=0 \]

$k=1,2,\cdots,m$。

  • 若拟合函数$f(x)=kx$, $f(x)=a+bx$等，可以直接代入前面的公式求解。

非线性最小二乘法

若函数$f(x)$关于要求的参数$\theta$是非线性的，

  1. 经过预处理后，函数可以变为关于参数是线性的
    • 若函数为$f(x)=Ce^{\frac1x}$，则可以对数据作处理，即拟合
      \[\ln y=\ln C+\frac1x , \ln y-\frac1x=C \]
    • 若想算$C$和$\alpha$，得到比例关系$y=Cx^{\alpha}$，则可以最小二乘
      \[\ln y=\ln C+\alpha \ln x \]
  2. 需要采用复杂的优化算法来求解。常用的算法有两类，一类是搜索算法，另一类是迭代算法。

搜索算法的思路是：

  • 按一定的规则选择若干组参数值，分别计算它们的目标函数值并比较大小；
  • 选出使目标函数值最小的参数值，同时舍弃其他的参数值；
  • 然后按规则补充新的参数值，再与原来留下的参数值进行比较，选出使目标函数达到最小的参数值。
  • 如此继续进行，直到选不出更好的参数值为止。

以不同的规则选择参数值，即可构成不同的搜索算法。常用的方法有单纯形搜索法、复合形搜索法、随机搜索法等。

迭代算法是从参数的某一初始猜测值$\theta$出发，然后产生一系列的参数点，如果这个参数序列收敛到使目标函数极小的参数点，那么经过足够步后，就可以认为得到了最小值。

迭代算法的一般步骤是：

  1. 给出初始猜测值$\theta$，并置迭代步数i=1。
  2. 确定一个向量$v$作为第i步的迭代方向。
  3. 用寻优的方法决定一个标量步长$\rho$
  4. 检查停机规则是否满足，如果不满足,则将i加1再从2开始重复；如果满足，则取$θ$为值。

典型的迭代算法有牛顿－拉夫森法、高斯迭代算法、麦夸特算法、变尺度法等。

非线性最小二乘法除可直接用于估计静态非线性模型的参数外，在时间序列建模、连续动态模型的参数估计中，也往往遇到求解非线性最小二乘问题。

Mathematica中的数据拟合

Mathematica中的拟合函数为 

Fit[data, funcs, vars]

其中data是二维数组，funcs是基函数列表，vars是变量名。如 

Fit[data,{1,x},x]  : 拟合y=a+bx
Exp[Fit[Log[data], {1,x}, x]]   : 拟合 y=a*x^b

非线性最小二乘法可以用 

FindFit[data,expr,pars,vars]

如 

FindFit[data, a x Log[b + c x], {a, b, c}, x]

例 3. 利用数据，将行星周期$T$与行星到太阳的平均距离$R$拟合为$T=cR^{\alpha}$

解. 对$\ln T$与$\ln R$作线性拟合，可以得到

\[\ln T=5.9101+1.4957 \ln R \]

mm2-ex2

mathematica: 

data = {{0.38, 88.0}, {0.72, 224.7}, {1.0, 365.3}, {1.52, 
    687.0}, {5.2, 4331.8}, {9.54, 10760.0}, {19.218, 30684.0}, {30.06,
     60188.3}};
f = Fit[Log[data], {1, x}, x]
fd = Plot[f, {x, -1.0, 4.0}];
pd = ListPlot[Log[data], PlotStyle -> Red];
Show[fd, pd]

结果 

5.91013 + 1.49568 x

mm2-ex2-mathematica

Python

Python中可以使用 scipy 库中的模块 optimize 的函数 leastsq() 作非线性的最小二乘法 

import numpy as np
from scipy.optimize import leastsq

Yi=np.array([7.01,2.78,6.47,6.71,4.1,4.23,4.05])
xi=np.array([8.19,2.72,6.39,8.71,4.7,2.66,3.78])
 
def func(p,x):
    k,b=p
    return k*x+b
 
def error(p,x,y,s):
    print(s)
    return func(p,x)-y
 
 #TEST
p0=[100,2]
#print( error(p0,Xi,Yi) )
s="Test the number of iteration"
Para=leastsq(error,p0,args=(xi,Yi,s)) #把error函数中除了p以外的参数打包到args中
k,b=Para[0]
print("k=",k,'\n',"b=",b)

车辆的停止距离

例 4. 3秒法则： 在高速上，跟前车的距离要大于当前车速跑3秒的距离。 这个法则合理吗？

分析：

问题： 车的当前速度为$v$，到车完全停止开了$h$的长度，$h$与$v$的函数关系是什么？

分析： 距离($h$)=司机反应距离($f$)+刹车距离($g$) 。

司机反应距离是司机反应时间与车速的函数，刹车距离是车速与车重的函数。

  • 假定在司机反应时间内，车保持车速行驶，则
    \[f=t_d v \]
    $t_d$是司机的反应时间。这个时间可以取平均偏高的时间。
  • 对于刹车距离来说，刹车时，相当于用最大的力$F$来使车完全停止，则力做的功与车的动能应该一样，有
    \[F g = \frac12 m v^2 \]
    其中$g$是刹车距离，$m$是车重，$v$是车速。可以看到$g$是$v^2$的比例关系。
    \[g=\frac{m}{2F}v^2 \]
    对于车重$m$与刹车力$F$的关系，可以假定它们成比例。这样$\frac{m}{2F}$就维持一个常数。

这样，有关系

\[h=c_0 v+ c_1 v^2 \]

其中$c_0$, $c_1$为待定常数。 需要用实测数据来估计$c_0$, $c_1$。

数据：实际测量。

结论：

3秒原则即要求$h<3v$。 可以看到，$h$与$v^2$相关，因此，3秒原则跟车速相关。也就是说，当车速较低时，可以用2秒原则；当车速较高时，可能需要4秒原则。

研究的对象是按比例放缩的。

例 5. 怎么估计鱼的重量？

假设：同一种类的鱼，它们的比例是一样，可以认为它们的重量与长度有关系。

模型：重量($W$)跟体积($V$)是比例关系，而体积($V$)与长度($l$)的3次方是比例关系。这样，可以有模型

\[W=k l^3 \]

模型二：若还测量鱼的最宽处的腰围($g$)，则体积($V$)与长度($l$)和橫截面($S$)相关，橫截面($S$)与腰围($g$)的平方是比例关系，则

\[W=k g^2 l \]

目录

本节读完

例 6. 谢

6.