最优化方法Python计算：一元函数导数的数值计算

定义1 给定连续函数$f(x)$，$x\in \Omega \subseteq \text{R}$。设$x,x_1\in \Omega$，
$$\Delta x=x-x_1$$
称为变量$x$的差分。此时，$x=x_1+\Delta x$。称
$$\frac{f(x_1+\Delta x)-f(x_1)}{\Delta x}$$
为$f(x)$在$x_1$处的前向差分。称
$$\frac{f(x_1)-f(x_1-\Delta x)}{\Delta x}$$
为$f(x)$在$x_1$处的后向差分。称
$$\frac{f(x_1+\Delta x/2)-f(x_1-\Delta x/2)}{\Delta x}$$
为$f(x)$在$x_1$处的中心差分。
若函数$f(x)$在$x_1$处可微，则$f(x)$在$x_1$处的导数为
$$f^{\prime }(x_1)=\underset{\Delta x\to 0}{\lim }\frac{f(x)-f(x_1)}{\Delta x}=\underset{\Delta x\to 0}{\lim }\frac{f(x_1+\Delta x)-f(x_1)}{\Delta x}$$
即$f(x)$在$x_1$处的导数$f^{\prime }(x_1)$是$f(x)$在$x_1$处的前向差分$\frac{f(x_1+\Delta x)-f(x_1)}{\Delta x}$当$x$的差分$\Delta x$趋于0时的极限。因此，当$|\Delta x|$充分小时，
$$f^{\prime }(x_1)\approx \frac{f(x_1)-f(x_1-\Delta x)}{\Delta x}.$$
在$f(x)$在$x_1$处的后向差分$\frac{f(x_1)-f(x_1-\Delta x)}{\Delta x}$中对变量差分$\Delta x$作变换$\Delta x^{\prime }=-\Delta x$，得
$$\frac{f(x_1)-f(x_1-\Delta x)}{\Delta x}=\frac{f(x_1)-f(x_1+{\Delta }^{\prime }x)}{-\Delta x^{\prime }}=\frac{f(x_1+\Delta x^{\prime })-f(x_1)}{\Delta x^{\prime }}$$
由于$\Delta x^{\prime }\to 0$当且进当$\Delta x\to 0$，故当$|\Delta x|$充分小时，有
$$f^{\prime }(x_1)\approx \frac{f(x_1)-f(x_1-\Delta x)}{\Delta x}.$$
对于$f(x)$在$x_1$处的中心差分
$$\begin{gathered} \frac{f(x_1+\Delta x/2)-f(x_1-\Delta x/2)}{\Delta x}=\frac{f(x_1+\Delta x/2)-f(x_1)+f(x_1)-f(x_1-\Delta x/2)}{\Delta x} \\ =\frac{1}{2}\frac{f(x_1+\Delta x/2)-f(x_1)+f(x_1)-f(x_1-\Delta x/2)}{\Delta x/2} \\ =\frac{1}{2}\left[\frac{f(x_1+\Delta x/2)-f(x_1)}{\Delta x/2}+\frac{f(x_1)-f(x_1-\Delta x/2)}{\Delta x/2}\right] \end{gathered}$$
所以
$$f^{\prime }(x_1)\approx \frac{f(x_1+\Delta x/2)-f(x_1-\Delta x/2)}{\Delta x}.$$
令
$$g(x)=\frac{f(x+\Delta x/2)-f(x-\Delta x/2)}{\Delta x}$$
称其为$f(x)$的广义导数。对$f(x)$的广义导数$g(x)$计算其在$x_1$处的中心差分，可得$f(x)$在$x_1$处的二阶导数近似计算公式
$$f^{\prime \prime }(x_1)\approx \frac{f(x_1+\Delta x)+f(x_1-\Delta x)-2f(x_1)}{\Delta x^2}.$$
利用$f^{\prime }(x_1)$和$f"(x_1)$的近似公式，实现一元函数二阶可微的一、二阶导数计算的Python函数定义如下

def der_scalar(f,x1):						#f表示函数，x1表示自变量值
   dx=1e-7									#设置Δx=0.0000001
   f1=(f(x1+dx/2)-f(x1-dx/2))/dx			#计算一阶导数
   f2=(f(x1+dx)+f(x1-dx)-2*f(x1))/dx**2		#计算二阶导数
   return f1,f2

程序中第1~5定义函数der_scalar。两个参数之一f表示二阶可微函数$f(x)$，x1表示自变量值$x_1$。第2行设置自变量增量$\Delta x=10^{-7}$为dx。第3行按$f^{\prime }(x_1)$的中心差分近似式计算$x_1$处的一阶导数$f^{\prime }(x_1)$为f1。第4行按广义导数的中心差分计算$x_1$处的二阶导数$f^{\prime \prime }(x_1)$为f2。
例1 考虑函数$f(x)=\sin x^2$，其一阶导数$f^{\prime }(x)=2x\cos x^2$，二阶导数$f^{\prime \prime }(x)=2(\cos x^2-2x^2\sin x^2)$。故对$x_1=\frac{\pi }{3}$可算得$f^{\prime }(x_1)=0.9563$，$f^{\prime \prime }(x_1)=-2.9893$。下列代码比较直接由导函数解析式计算和调用上列程序定义的der_scalar函数计算的结果。

import numpy as np									#导入numpy
f=lambda x:np.sin(x**2)								#设置函数f(x)
f1=lambda x:2*x*np.cos(x**2)						#设置一阶导数f'(x)
f2=lambda x:2*(np.cos(x**2)-2*np.sin(x**2)*x**2)	#设置二阶导数f"(x)
x1=np.pi/3											#设置自变量x1
print((f1(x1),f2(x1)))								#解析计算
print(der_scalar(f,x1))								#数值计算

程序的2、3、4行分别设置$f(x)$、$f^{\prime }(x)$和$f^{\prime \prime }(x)$的解析式。第5行设置自变量的值$x_1=\frac{\pi }{3}$。第6行输出解析计算结果，第7行调用der_scalar函数数值地计算$f^{\prime }(x_1)$和$f^{\prime \prime }(x_1)$。运行程序，输出

(0.9563079109495481, -2.9893240816612874)
(0.9563079098976839, -3.0000693287648676)

输出的第1行为解析计算的结果，第2行为数值计算结果。比较两者可见der_scalar函数计算的一阶导数结果精度是很高的，二阶导数的计算精度稍差。这是因为，这是对“广义导数”进行差分计算的结果。换句话说，在计算二阶导数时，累积了两次计算误差，产生了误差的“放大”效应。尽管如此，der_scalar函数在大多数场合都能满足应用的要求。
例2: 下图给出了移动无线通信系统的简化模型。有两个同为$1$个高度单位相矩$2$个长度单位的相邻基站，手机用户所处位置为$x$，手机收到的基站信号的功率与用户与基站之间的距离平方成反比。试确定用户最为合适的位置，使得信号干扰比最大，即用户接收到来自基站$1$的信号功率与来自基站$2$的信号功率之间的比值最大。

解：按题意，手机收到的基站信号的功率与用户与基站之间的距离平方成反比。即收到基站1和基站2的功率分别为
$$\begin{array}{l}{w}_1=\frac{1}{1+x^2}\\ w_2=\frac{1}{1+(2-x{)}^2}\end{array}$$
为使得信号干扰比最大，即用户接收到来自基站1的信号功率与来自基站2的信号功率之间的比值最大，求函数$f(x)=\frac{w_1}{w_2}=\frac{1+(2-x{)}^2}{1+x^2}$的最大值点$x_0$。
由于$f(x)$二阶连续可微，为计算$f(x)$的最大值点，先求得$f(x)$驻点。
$$\begin{gathered} f^{\prime }(x)=\frac{-2(2-x)(1+x^2)-2x(1+(2-x{)}^2)}{(1+x^2{)}^2} \\ =\frac{4(x^2-2x-1)}{(1+x^2{)}^2}=\frac{4((x-1{)}^2-2)}{(1+x^2{)}^2} \\ =\frac{4((x-1-\sqrt{2})(x-1+\sqrt{2}))}{(1+x^2{)}^2} \end{gathered}$$
令$f^{\prime }(x)=0$，即得驻点$x_1=1-\sqrt{2}$和$x_2=1+\sqrt{2}$。按二阶可微函数极值点的必要条件，需计算$f(x)$的二阶导数$f^{\prime \prime }(x)$在驻点处的值，通过二阶导数值的正负以判断是极大值点或是极小值点。下列代码完成计算过程。

import numpy as np					#导入numpy
f=lambda x:(1+(2-x)**2)/(1+x**2)	#设置函数f(x)
x1=1-np.sqrt(2)						#驻点x1
x2=1+np.sqrt(2)						#驻点x2
print('%.1f,%.4f'%der_scalar(f,x1))	#计算x1处一、二阶导数
print('%.1f,%.4f'%der_scalar(f,x2))	#计算x2处一、二阶导数

利用代码中的注释信息不难理解程序。仅提请注意第3、4行分别设置的是驻点$x_1=1-\sqrt{2}$和$x_2=1+\sqrt{2}$。运行程序，输出：

0.0,-8.2107
0.0,0.2566

第1行输出的是$f(x)$在$x_1=1-\sqrt{2}$处的一、二阶导数值，说明$x_1$确实为驻点，且由于在$x_1$处的二阶导数值近似为-8.2107，根据极值点的充分条件知$x_1$为$f(x)$的唯一极大值点，即为本题所求。第二行的输出表明$x_2=1+\sqrt{2}$是$f(x)$的极小值点，也是最小值点。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！