TASK1_高等数学

高等数学与数值计算

(一)最优化问题

一元函数

❓ GitModel 公司发明了一种专用于数值计算的平板电脑 GitNum$,$ 考虑到员工生产效率的问题$,$ GitNum 的生产成本是与生产台数相关的. 函数关系为

$$C(n)=1250n(2-e^{-(n-5000{)}^2}).$$

请问如何确定这批 GitNum 的生产台数使得每台的平均生产成本最小使得经济效益最大呢?

分析： 求平均生产成本最小$,$ 本质上就是求平均成本函数 $\overline{C}(n)=\frac{C(n)}{n}$ 的最小值. 事实上$,$ 初中我们就经常接触类似的最值问题$,$ 只不过所谓的目标函数往往是最简单的二次函数. 而在高中之后我们学习了求导工具之后发现对任意可导的函数我们都可以探究其导数异号零点来确定其极值点.

诚然$,$ 仅仅是导数零点并不足以说明其就是极值点$,$ 如下图的函数 $y=x^3,$ 我们发现 $x=0$ 是其导数零点但在图像上仅仅只是函数的一个$"$拐点$"$.

*极值vs最值

这是因为导函数 $y^{\prime }=3x^2$ 在导数零点 $x=0$ 处仅仅是$"$切过$"$ $x$ 轴而非$"$穿过$"$ $x$ 轴! 而分析学告诉我们$,$ 导函数$f$ 在零点 $x_0$ 处 $"$穿过$"$ $x$ 轴的图像性质本质上就是 $f$ 在零点处附近是单调的$,$ 亦即二阶导函数 $($导数的导数$)$ $f^{\prime \prime }(x_0)>0$ 或 $f^{\prime \prime }(x_0)<0$!

回到问题上来$,$ 要求 $\overline{C}(n)$ 的极值点$,$ 我们首先需要求出导函数 ${\overline{C}}^{\prime }(n)$ 的零点 $n_0,$ 如果 ${\overline{C}}^{\prime \prime }(n_0)>0,$ 那么 $n_0$ 即为 $\overline{C}(n)$ 的极小值点.

代码如下 :

from sympy import *
n = symbols('n')
y = 1250*(2-exp(-(n-5000)**2))

func1 = diff(y,n) # 求一阶导数
func1

$-1250\cdot \left(10000-2n\right)e^{-{\left(n-5000\right)}^2}$

stagnation = solve(diff(y,n),n)
print("该函数驻点为",stagnation) # 计算驻点

该函数驻点为 [5000]

func2 = diff(y, n, 2) # # 计算二阶导数
func2

$2500\cdot \left(1-2{\left(n-5000\right)}^2\right)e^{-{\left(n-5000\right)}^2}$

# 计算驻点的二阶导数值，验证正负
print(func2.evalf(subs = {n:5000}))
func2.evalf(subs = {n:5000}) > 0

2500.00000000000
True

# 函数的极小值
y.evalf(subs = {n:5000})

1250.00000000000

经验证$,$ $n_0=5000$ 确实为 $\overline{C}(n)$ 的极小值点$,$ 即第一批 GitNum 生产台数计划为 $5000$ 台时平均成本最小$,$ 为 $1250$ 元/台.

抽象提炼

数学建模中$,$ 优化问题是渗透到各个方面的$,$ 小到最优参数的确定$,$ 大到最优策略的规划. 每一个优化问题都可以以如下标准形式给出 :
$$\begin{gathered} maxf(x) \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{l}{g}_i(x)⩾0,i=1,2,\cdots ,n\\ h_j(x)=0,j=1,2,\cdots ,m\end{array} \end{gathered}$$
其中 $f(x)$ 即是我们需要求最值的函数$,$ $x$ 是我们用来 $"$决策$"$ 的变量. $g_i(x)⩾0,h_j(x)=0$ 分别是决策变量 $x$ 满足的条件$,$ 即决策变量 $x$ 的约束. 当决策变量 $x$ 没有约束时$,$ 称为无约束优化问题.

上文提到的GitNum生产台数问题即为一个最简单的无约束优化问题——目标函数有明确表达式，二阶以上可导，定义域离散程度不高. 往往生产一批产品的数量范围是足够大的以至于我们可以将离散的整数视为连续的. 对于这样的简单优化问题$,$ 利用数学分析的知识可以做到精确的分析. 我们给出优化领域最经典也是最实用的判别条件.

命题1:[一阶最优化必要条件]

设 $f(x)$ 是 $\mathbf{R}$ 上的可导函数. 若 $x^{\ast }$ 是 $f$ 的极值点则
$$f^{\prime }(x^{\ast })=0$$
其中 $f^{\prime }$ 是 $f$ 的一阶导数.

命题2:[二阶最优化必要条件]

命题3:[二阶最优化充分条件]
设 $f(x)$ 是 $\mathbf{R}$ 上的二阶可导函数. 若 $x^{\ast }$ 满足
$$f^{\prime }(x^{\ast })=0,f^{\prime \prime }(x^{\ast })<0(f^{\prime \prime }(x^{\ast })>0),$$
则 $x^{\ast }$ 是 $f$ 的极大$($小$)$值点.

多元函数

❓ GitModel 公司通过 GitNum 的发售赚到第一桶金后马不停蹄地投入到 GitNum 的改良研制中$,$ 经市场调查后 GitModel 公司发现根据用户需求可推出轻便版以及高性能版的 GitNum. 三种 GitNum:

在这一批生产中$,$ GitModel 公司应该制定怎么样的生产计划即三种 GitNum 的生产台数才能使得利润与成本的比值最大化?

我们把三种 GitNum 的生产台数按上表顺序标记序号 $n_i,i=1,2,3,$ 分别用 $R,Q,C$ 来代表收入$,$ 利润以及成本. 问题旨在求出利润/成本 $r=\frac{Q}{C}$ 的最大值$,$ 而我们需要同时考虑三种 GitNum 的生产台数$,$ 也就是说$,$ 需要求极值的函数 $r$ 是 $n_1,n_2,n_3$ 的多元函数 !

总收入函数为 $$R(n_1,n_2,n_3)=2000n_1+3000n_2+4500n_3,$$
总成本函数为
$$C(n_1,n_2,n_3)=750n_1(2-e^{-(n_1-6000{)}^2})+1250n_2(2-e^{-(n_2-5000{)}^2})+2000n_3(2-e^{-(n_3-3000{)}^2}),$$
总利润函数为
$$Q(n_1,n_2,n_3)=R(n_1,n_2,n_3)-C(n_1,n_2,n_3),$$
所以总利润与总成本比值为

\begin{align}
r(n_1,n_2,n_3)&=\dfrac{R(n_1,n_2,n_3)-C(n_1,n_2,n_3)}{C(n_1,n_2,n_3)}\
&=\dfrac{R(n_1,n_2,n_3)}{C(n_1,n_2,n_3)}-1\
&=\dfrac{2000n_1+3000n_2+4500n_3}{750n_1(2-e^{-(n_1-6000)2})+1250n_2(2-e^{-(n_2-5000)2})+2000n_3(2-e^{-(n_3-3000)2})}-1.
\end{align}

优化问题即为
$$\underset{n_1,n_2,n_3⩾1}{\max }r(n_1,n_2,n_3).$$

依样画葫芦$,$ 该问题不过是第一个问题的变量 Plus 版$,$ 我们仍然可以采用研究导数以及二阶导数性态的方式去求出极大值组合 $(n_1^{\ast },n_2^{\ast },n_3^{\ast }).$

那有的同学可能就会问了 : 一元函数分析跟多元函数分析有什么区别呢?

在多元函数 $f$ 中$,$ 我们只能对单个变量 $x$ 求偏导数 $f_x,$ 这远远不足以刻画导数的全部信息. 所以$,$ 梯度 $($全导数$)$ 以及 Hesse 矩阵的概念替代了一元函数中导数以及二阶导数.

定义1:[Hesse 矩阵] 设 $f:\mathbf{R}^n\to R$ 是 $n$ 元函数$,$ $x_i,i=1,2,\cdots,n$ 是变量.

记 $f$ 对变量 $x_i$ 方向的偏导数为 $f_i^{\prime },$ 称所有偏导数组成的列向量 $[f_1^{\prime },f_2^{\prime },\cdots ,f_n^{\prime }{]}^T:=\nabla f,$ 为函数 $f$ 的全导数$,$ 亦称梯度.

记偏导数 $f_i^{\prime }$ 对变量 $x_j$ 方向的二阶偏导数为 $f_{ij}^{\prime \prime }$ 称 $n\times n$ 方阵$$\left[\begin{array}{cccc}{f}_{11}^{\prime \prime }& f_{12}^{\prime \prime }& \cdots & f_{1n}^{\prime \prime }\\ f_{21}^{\prime \prime }& f_{22}^{\prime \prime }& \cdots & f_{1n}^{\prime \prime }\\ ⋮& ⋮& \ddots & ⋮\\ f_{n1}^{\prime \prime }& f_{n2}^{\prime \prime }& \cdots & f_{nn}^{\prime \prime }\end{array}\right]:={\nabla }^{2}f,$$
为函数 $f$ 的 Hesse 矩阵.

事实上$,$ 通过 Taylor 展开我们即可验证猜想的正确性$,$ 理论证明不难而繁$,$ 我们只保留最精华实用的命题 : 最优化条件对应的多元函数版本则如下 :

命题4:[多元无约束问题一阶最优化必要条件]

命题5:[多元无约束问题二阶最优化必要条件]
命题6:[多元无约束问题二阶最优化充分条件]
设 $f(x)$ 是 ${\mathbf{R}}^n$ 上的 $n$ 元二阶可导函数. 若 $x^{\ast }$ 满足
$$\nabla f(x^{\ast })=0,{\nabla }^{2}f(x^{\ast })\text{负定}({\nabla }^{2}f(x^{\ast })\text{正定}),$$
则 $x^{\ast }$ 是 $f$ 的极大$($小$)$值点.

而代数学的知识告诉我们,验证一个对称矩阵是否正(负)定只需要 check 其所有特征值是否大(小)于0.

所以$,$ 要制定最佳的生产计划 $(n_1^{\ast },n_2^{\ast },n_3^{\ast })$ 使得利润与成本比值最大$,$ 即求解问题(2)$,$ 只需要求出$3$元方程组
$$\nabla r=0$$
的解集$,$ 再挑出解集中使得 ${\nabla }^{2}r$ 负定的解$,$ 这就是我们的最佳生产计划 $(n_1^{\ast },n_2^{\ast },n_3^{\ast })$!

最优化条件给出了这样一个基本事实 : 不管函数 $f$ 表达式复杂与否$,$ 只要二阶以上可导$,$ 我们只需要找出满足 Hesse 矩阵 ${\nabla }^{2}f$ 半负定$($半正定$)$的梯度零点$,$ 极大$($小$)$值点必定存在其中!

现实生活中我们遇到的许多优化问题需要解的方程都是超越方程$,$ 早已超出人力计算的范畴$,$ 而我们衍生的许多数值计算方法就是为了求出这类方程的近似解的. 优化问题中一个最经典的求极值算法便是针对优化条件发明的梯度下降法——通过迭代的方式找到一列梯度递减的点,当点的**梯度下降的足够接近$0$**时便可认为该点即是极值的候选之一.

scipy求解多元函数的极值

import scipy.optimize as opt
from scipy import fmin
import numpy as np

def func0(cost, x, a):
    return cost*x*(2 - exp(-(x - a)**2))
func = lambda x: (2000*x[0] + 3000*x[1] + 4500*x[2]) / (func0(750, x[0], 6000) + func0(1250, x[1], 5000) + func0(2000, x[2], 3000)) - 1 
bnds = ((1, 10000), (1, 10000), (1, 10000))
res = opt.minimize(fun=func, x0=np.array([2, 1, 1]), bounds=bnds)
res

      fun: 0.126003456985265
 hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
      jac: array([ 6.24167384e-04,  3.71258579e-04, -8.06021916e-06])
  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 64
      nit: 15
   status: 0
  success: True
        x: array([  1.        ,   1.        , 123.56936554])

(二) : 以插值知识解决数据处理问题

❓ GitModel 公司工作室刚刚建完，准备安装宽带以及路由，基于传输数据要求，GitModel 需要的宽带运营商网速越快越好，而由于网速是一个随时间变化的非固定量，简单衡量平均网速也并非一个好的评价手段, GitModel 准备收集在同一个地点各个运营商的宽带一天内的网速数据 以做接下来的建模分析, 出于人力考虑,网速监视每小时汇报一次数据. A运营商的宽带24小时网速如下表:

然而$,$ 应该注意到离散的数据点无法完美体现网速的整体性质$,$ 我们需要通过仅有的数据点还原宽带一天内的实时网速. 以 $A$ 运营商为例$,$ 请问 GitModel 公司应如何$"$模拟$"$出网速曲线呢?

该问题说数学话版$,$ 就是需要找出一条经过上表 $24$ 个点 $(t_i,s_i),i=1,2,\cdots ,24$ 的曲线 $s=s(t),$ 我们自然会问$,$ 那怎么构造 $s(t)$ 合适呢? 又应该怎么说明构造的 $s(t)$ 能有较好的模拟效果呢?

我们接着往下分析 : 由于宽带用户千千万万$,$ 大基数导致了使用宽带的用户数关于时间 $t$ 的变化是近似连续的$,$ 从而网速 $s(t)$ 可以假设为 $t$ 的连续函数$,$ 更甚之$,$ 可以假设为斜率以及曲率 $s^{\prime }(t),s^{\prime \prime }(t)$ 连续! 我们将 $s(t)$ 的满足条件罗列如下 :

• $s(t)$ 过 $24$ 个数据点 $s_i=s(t_i),i=1,2,\cdots ,24.$
• $s(t)$ 连续.
• $s^{\prime }(t)$ 连续.
• $s^{\prime \prime }(t)$ 连续.

$s(t)$ 应该满足的条件暂且如此$,$ 我们再进行思考 : 有什么构造 $s(t)$ 方式是比较简单的呢?

在高中我们就已经学过了给定一组数据点可通过最小二乘法来拟合出一条回归直线$,$ 诚然$,$ 两点一直线的数学直觉告诉我们能过 $24$ 个数据点的直线几乎不存在$,$ 那么我们可否放宽条件$,$ 构造出过 $24$ 个数据点的折线呢? 这显然是可以的! 过数据点的条件可由分段函数来解决.

(三) : 以积分知识解决工程问题

❓ 鉴于 GitModel 公司刚上市$,$ 承接了许多数学建模相关业务. 有一家工程企业给出了这样一个任务 : 该企业准备制造一批抽水机$,$ 需要得到其各种数据来决定抽水机的性能参数. 企业首先对圆柱形的储水桶进行了抽水测试$,$ 需要求出将桶内的水全部抽出需要做多少功？

from scipy import integrate # 已知函数表达式积分
from scipy import pi
def f(h):
    return 88.2 * pi * (5 - h)
v, err = integrate.quad(f, 0, 5) # 被积函数与积分区间
v

3463.605900582747

📋 小结：
本质上 , 微元法其实就是分析学中积分的应用 : 将连续变化的函数 $($变力做功$,$ 恒力变位移做功$,$ 连续型随机变量期望值等等连续物理量$)$ 先微分再积分得到相关量$,$ 这在处理非离散或离散程度小的建模问题中尤为常见.

❓ 紧接着 GitModel 公司又接到另一家工程企业的任务 : 由于吊桥的年久失修导致一些铁链生锈而质量发生了改变$,$ 为保证安全必须重新测量铁链的总体质量$,$ 而体积巨大难以拆卸重组的铁链无法直接测量其质量$,$ 仅能靠检测每块铁环的质量来估计整体变化的质量$,$ 然而逐块检测是一项巨大的工程$,$ 需要耗费巨额的时间与人力.

据专家分析$,$ 若以桥梁中心为原点建立空间坐标系$,$ 铁链的 $y$ 轴坐标可以近似相同. 铁链上每块铁环的密度仅与其所处的横向距离以及海拔 $(x,z)$ 有关$,$ 经数值拟合可得到密度 $(\mathrm{k}\mathrm{g}/{\mathrm{m}}^3)$ 与位置 $(\mathrm{m})$ 的函数为
$$\rho (x,z)=7860\left(1+1.5^{2-\frac{z}{10}-0.1{\left(\frac{x}{50}\right)}^2}\right)$$ 及铁链的垂直面曲线方程为
$$z=30\cosh \frac{x}{300},$$ 铁环是圆柱形的$,$ 半径为 $r=0.15\mathrm{m}$. GitModel 公司如何通过不直接检测的方式来估计铁链的质量呢?

我们发现$,$ 这其实也是类似于上文问题的离散程度小的问题 $($因为铁链的铁环数足够多使得我们可以近似于连续问题$),$ 每块铁环可以当作圆柱形状的微元$,$ 其具有密度与位置的函数 $\rho (x,z),$ 长度微分 $ds$ 以及体积微分 $dV$ 那么我们同样可以应用积分来求出质量.

那么该问题是否可以照搬上文的方法呢? 我们接着往下分析两个问题的不同. 应该注意到 : 第一个问题是位移在垂直直线上的分布, 而该问题是密度 $($质量$)$ 在曲线上的分布 $,$ 这代表着我们无法通过常规积分的方式来求$,$ 因为我们的坐标轴不是一条直线!

当然$,$ 分析学同样提供了在曲线上的积分方式$,$ 我们称为曲线积分. 本问题就是一个经典的第一型曲线积分问题——知道曲线上每个点 $(x,z)$ 上的分布函数 $\rho (x,z),$ 求曲线的质量 $M={\int }_{L}m(x,z)ds,$ 其中$,$ 类似于我们对坐标 $x$ 做微分 $dx,$ 曲线上的坐标可以用长度来表示$,$ 所以我们对曲线的长度做了微分 $ds$ 后进行积分处理. 与直线积分有所不同的是$,$ 我们无法以曲线为轴$,$ 但仍然可以将曲线放在坐标系中$,$ 此时由勾股定理我们可以求出长度微分的表达式$,$ 如下

下面给出第一型曲线积分python代码，铁链质量的估算值为 $98635\mathrm{t}$

# 第一型曲线积分
from scipy import sinh, cosh, sqrt
def f(x):
    index1 = (2 - 3 * cosh(x/300) - 0.1 * ((x/50)**2))
    return (1 + 1.5 ** index1) * sqrt(100 + (sinh(x/300))**2)

v, err = integrate.quad(f, -300, 300)
v * 11.79

98635.09908278256

❓ 有了直线积分$,$ 曲线积分$,$ 我们自然会问 : 有没有平面积分以及曲面积分呢? 当然! 事实上$,$ 由于我们还会接触变力曲线做功$,$ 这代表在每个点上的分布函数不再是标量函数而是向量函数! 分析学还给出了第二类曲线积分以及曲面积分.

有一位爱钻牛角尖的同学发出了提问——老师$,$ 这样的方法还是太$"$数学$"$了$,$ 并且理想很丰满$,$ 现实很骨感$,$ 万一没有函数表达式呢?又万一有函数表达式但积分太难算了呢?

别急$,$ 数学家也想到了应用性的问题$,$ 理论总是要落地的嘛! 所以有了计算机之后$,$ 许多数值计算的方法也被发明了出来$,$ 其中就包括了数值积分的方法——计算机可以通过数值方法对许多不连续函数甚至是无法写出表达式的函数进行积分! 这便很好地解决了积分在实际应用上的问题.这将在我们的2.0课程中与大家见面👋

实战项目——人口增长问题

💼 任务来袭

搜集 $1950\sim 2020$ 年间美国人口数据$,$ 猜测其满足的函数关系$,$ 并综合数据预测美国 $2030$ 年的人口数.

公司希望实习生就以下的开放性问题给出自己的想法，公司要求实习生提交预测的思路$,$ 模型$,$ 算法以及结果.

面试官给出了如下提示 : 预测值与真实值存在误差$,$ 该问题如何转化为可用上述所学的知识解决的问题呢?

从你的数据来看$,$你猜想美国人口与时间的关系应该满足怎样的函数关系$($线性$,$二次$,$或是其他$)$? 这样的函数关系你如何确定? 又怎么衡量这个函数具有良好的预测效果?[请在下方Markdown方格中写下你的答案]

import pandas  as pd
popu = pd.read_excel('popu.xlsx')
popu.head()

	year	population
0	1950	151684000
1	1951	154287000
2	1952	156954000
3	1953	159565000
4	1954	162391000

popu.describe()

	year	population
count	71.000000	7.100000e+01
mean	1985.000000	2.422277e+08
std	20.639767	5.297308e+07
min	1950.000000	1.516840e+08
25%	1967.500000	1.997090e+08
50%	1985.000000	2.379240e+08
75%	2002.500000	2.888665e+08
max	2020.000000	3.294840e+08

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
years = popu['year']
popu_num = popu['population']
plt.scatter(years, popu_num)
plt.xlabel('Year')
plt.ylabel('Population#')
plt.show()

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-RY8k5Kxi-1655307548402)(output_36_0.png)]

from scipy import optimize
def test_func(x,a,b):
    return a*x+b
p,p_=optimize.curve_fit(test_func,years,popu_num)
print(p)

[ 2.56360259e+06 -4.84652345e+09]

# 最后，将拟合的模型对应到真实的数据中
def y__(x):
    return res.x[0] * x + res.x[1]
y_pre = y__(years)
plt.plot(years, popu_num, '-', years, y_pre, '--')
plt.xlabel('Year')
plt.ylabel('People Number')
plt.legend(['real data', 'predict data'])
plt.show()

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ROs1TTlH-1655307548405)(output_38_0.png)]

# 预测2020年和2021年的美国人口
year_2020 = y.evalf(subs={x:2020})
year_2021 = y.evalf(subs={x:2021})
print(year_2020)
print(331761000)
print("美国2020年的人口总数预测结果为为：{:.0f}，与实际值的误差为：{:.0f}".format(year_2020, np.abs(year_2020-331761000.0)))
print("美国2021年的人口总数预测结果为为：{:.0f}，与实际值的误差为：{:.0f}".format(year_2021, np.abs(year_2021-332213000)))
# 预测美国2030年人口总数
year_2030 = y.evalf(subs={x:2030})
print("美国2030年的人口总数预测结果为为：{:.0f}。".format(year_2030))

331953780.699434
331761000
美国2020年的人口总数预测结果为为：331953781，与实际值的误差为：192781
美国2021年的人口总数预测结果为为：334517383，与实际值的误差为：2304383
美国2030年的人口总数预测结果为为：357589807。

ref:https://blog.csdn.net/qq_44341554/article/details/125281730