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

定义1 给定连续函数$f(x)$，$x\in \Omega \subseteq \text{R}$。设 x , x 1 ∈ Ω x,x_1\\in\\Omega x,x1​∈Ω，
Δ x = x − x 1 \\Delta x=x-x_1 Δx=x−x1​
称为变量$x$的差分。此时， x = x 1 + Δ x x=x_1+\\Delta x x=x1​+Δx。称
f ( x 1 + Δ x ) − f ( x 1 ) Δ x \\frac{f(x_1+\\Delta x)-f(x_1)}{\\Delta x} Δxf(x1​+Δx)−f(x1​)​
为$f(x)$在 x 1 x_1 x1​处的前向差分。称
f ( x 1 ) − f ( x 1 − Δ x ) Δ x \\frac{f(x_1)-f(x_1-\\Delta x)}{\\Delta x} Δxf(x1​)−f(x1​−Δx)​
为$f(x)$在 x 1 x_1 x1​处的后向差分。称
f ( x 1 + Δ x / 2 ) − f ( x 1 − Δ x / 2 ) Δ x \\frac{f(x_1+\\Delta x/2)-f(x_1-\\Delta x/2)}{\\Delta x} Δxf(x1​+Δx/2)−f(x1​−Δx/2)​
为$f(x)$在 x 1 x_1 x1​处的中心差分。
若函数$f(x)$在 x 1 x_1 x1​处可微，则$f(x)$在 x 1 x_1 x1​处的导数为
f ′ ( x 1 ) = lim ⁡ Δ x → 0 f ( x ) − f ( x 1 ) Δ x = lim ⁡ Δ x → 0 f ( x 1 + Δ x ) − f ( x 1 ) Δ x f'(x_1)=\\lim_{\\Delta x\\rightarrow0}\\frac{f(x)-f(x_1)}{\\Delta x}=\\lim_{\\Delta x\\rightarrow0}\\frac{f(x_1+\\Delta x)-f(x_1)}{\\Delta x} f′(x1​)=Δx→0lim​Δxf(x)−f(x1​)​=Δx→0lim​Δxf(x1​+Δx)−f(x1​)​
即$f(x)$在 x 1 x_1 x1​处的导数 f ′ ( x 1 ) f'(x_1) f′(x1​)是$f(x)$在 x 1 x_1 x1​处的前向差分 f ( x 1 + Δ x ) − f ( x 1 ) Δ x \\frac{f(x_1+\\Delta x)-f(x_1)}{\\Delta x} Δxf(x1​+Δx)−f(x1​)​当$x$的差分$\Delta x$趋于0时的极限。因此，当$|\Delta x|$充分小时，
f ′ ( x 1 ) ≈ f ( x 1 ) − f ( x 1 − Δ x ) Δ x . f'(x_1)\\approx \\frac{f(x_1)-f(x_1-\\Delta x)}{\\Delta x}. f′(x1​)≈Δxf(x1​)−f(x1​−Δx)​.
在$f(x)$在 x 1 x_1 x1​处的后向差分 f ( x 1 ) − f ( x 1 − Δ x ) Δ x \\frac{f(x_1)-f(x_1-\\Delta x)}{\\Delta x} Δxf(x1​)−f(x1​−Δx)​中对变量差分$\Delta x$作变换 Δ x ′ = − Δ x \\Delta x'=-\\Delta x Δx′=−Δx，得
f ( x 1 ) − f ( x 1 − Δ x ) Δ x = f ( x 1 ) − f ( x 1 + Δ ′ x ) − Δ x ′ = f ( x 1 + Δ x ′ ) − f ( x 1 ) Δ x ′ \\frac{f(x_1)-f(x_1-\\Delta x)}{\\Delta x}=\\frac{f(x_1)-f(x_1+\\Delta'x)}{-\\Delta x'}=\\frac{f(x_1+\\Delta x')-f(x_1)}{\\Delta x'} Δxf(x1​)−f(x1​−Δx)​=−Δx′f(x1​)−f(x1​+Δ′x)​=Δx′f(x1​+Δx′)−f(x1​)​
由于 Δ x ′ → 0 \\Delta x'\\rightarrow0 Δx′→0当且进当$\Delta x\to 0$，故当$|\Delta x|$充分小时，有
f ′ ( x 1 ) ≈ f ( x 1 ) − f ( x 1 − Δ x ) Δ x . f'(x_1)\\approx\\frac{f(x_1)-f(x_1-\\Delta x)}{\\Delta x}. f′(x1​)≈Δxf(x1​)−f(x1​−Δx)​.
对于$f(x)$在 x 1 x_1 x1​处的中心差分
f ( x 1 + Δ x / 2 ) − f ( x 1 − Δ x / 2 ) Δ x = f ( x 1 + Δ x / 2 ) − f ( x 1 ) + f ( x 1 ) − f ( x 1 − Δ x / 2 ) Δ x = 1 2 f ( x 1 + Δ x / 2 ) − f ( x 1 ) + f ( x 1 ) − f ( x 1 − Δ x / 2 ) Δ x / 2 = 1 2 [ f ( x 1 + Δ x / 2 ) − f ( x 1 ) Δ x / 2 + f ( x 1 ) − f ( x 1 − Δ x / 2 ) Δ x / 2 ] \\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] Δxf(x1​+Δx/2)−f(x1​−Δx/2)​=Δxf(x1​+Δx/2)−f(x1​)+f(x1​)−f(x1​−Δx/2)​=21​Δx/2f(x1​+Δx/2)−f(x1​)+f(x1​)−f(x1​−Δx/2)​=21​[Δx/2f(x1​+Δx/2)−f(x1​)​+Δx/2f(x1​)−f(x1​−Δx/2)​]
所以
f ′ ( x 1 ) ≈ f ( x 1 + Δ x / 2 ) − f ( x 1 − Δ x / 2 ) Δ x . f'(x_1)\\approx\\frac{f(x_1+\\Delta x/2)-f(x_1-\\Delta x/2)}{\\Delta x}. f′(x1​)≈Δxf(x1​+Δx/2)−f(x1​−Δx/2)​.
令
g ( x ) = f ( x + Δ x / 2 ) − f ( x − Δ x / 2 ) Δ x g(x)=\\frac{f(x+\\Delta x/2)-f(x-\\Delta x/2)}{\\Delta x} g(x)=Δxf(x+Δx/2)−f(x−Δx/2)​
称其为$f(x)$的广义导数。对$f(x)$的广义导数$g(x)$计算其在 x 1 x_1 x1​处的中心差分，可得$f(x)$在 x 1 x_1 x1​处的二阶导数近似计算公式
f ′ ′ ( x 1 ) ≈ f ( x 1 + Δ x ) + f ( x 1 − Δ x ) − 2 f ( x 1 ) Δ x 2 . f''(x_1)\\approx\\frac{f(x_1+\\Delta x)+f(x_1-\\Delta x)-2f(x_1)}{\\Delta x^2}. f′′(x1​)≈Δx2f(x1​+Δx)+f(x1​−Δx)−2f(x1​)​.
利用 f ′ ( x 1 ) f'(x_1) f′(x1​)和 f " ( x 1 ) f"(x_1) f"(x1​)的近似公式，实现一元函数二阶可微的一、二阶导数计算的Python函数定义如下

def der_scalar(f,x1):						#f表示函数，x1表示自变量值dx=1e-7									#设置Δx=0.0000001f1=(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 x_1 x1​。第2行设置自变量增量 Δ x = 1 0 − 7 \\Delta x=10^{-7} Δx=10−7为dx。第3行按 f ′ ( x 1 ) f'(x_1) f′(x1​)的中心差分近似式计算 x 1 x_1 x1​处的一阶导数 f ′ ( x 1 ) f'(x_1) f′(x1​)为f1。第4行按广义导数的中心差分计算 x 1 x_1 x1​处的二阶导数 f ′ ′ ( x 1 ) f''(x_1) f′′(x1​)为f2。
例1 考虑函数 f ( x ) = sin ⁡ x 2 f(x)=\\sin{x^2} f(x)=sinx2，其一阶导数 f ′ ( x ) = 2 x cos ⁡ x 2 f'(x)=2x\\cos{x^2} f′(x)=2xcosx2，二阶导数 f ′ ′ ( x ) = 2 ( cos ⁡ x 2 − 2 x 2 sin ⁡ x 2 ) f''(x)=2(\\cos{x^2}-2x^2\\sin{x^2}) f′′(x)=2(cosx2−2x2sinx2)。故对 x 1 = π 3 x_1=\\frac{\\pi}{3} x1​=3π​可算得 f ′ ( x 1 ) = 0.9563 f'(x_1)=0.9563 f′(x1​)=0.9563， f ′ ′ ( x 1 ) = − 2.9893 f''(x_1)=-2.9893 f′′(x1​)=−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 ′ ( x ) f'(x) f′(x)和 f ′ ′ ( x ) f''(x) f′′(x)的解析式。第5行设置自变量的值 x 1 = π 3 x_1=\\frac{\\pi}{3} x1​=3π​。第6行输出解析计算结果，第7行调用der_scalar函数数值地计算 f ′ ( x 1 ) f'(x_1) f′(x1​)和 f ′ ′ ( x 1 ) f''(x_1) f′′(x1​)。运行程序，输出

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

输出的第1行为解析计算的结果，第2行为数值计算结果。比较两者可见der_scalar函数计算的一阶导数结果精度是很高的，二阶导数的计算精度稍差。这是因为，这是对“广义导数”进行差分计算的结果。换句话说，在计算二阶导数时，累积了两次计算误差，产生了误差的“放大”效应。尽管如此，der_scalar函数在大多数场合都能满足应用的要求。