【Python讲优化】S06E07 一元函数的极值:运用迭代法

在上一讲的末尾我们谈到,在实际的工程当中我们常常借助计算机程序,利用迭代法进行极值的求取,这里我们首先从一元函数入手,看看如何通过这种方法找到一元函数的极值点。

1.起步:用牛顿法解方程

1.1.原理分析

在介绍求取函数$f(x)$的极值方法前,我们首先谈一下解方程的问题。

在解一元函数的高阶方程,形如$ax^n+bx^{n-1}+cx^{n-2}+...+1=0$时,大家肯定会想到用因式分解或者是用求根公式来获取方程的解析解,但是事实总不会是那么一帆风顺的,有时候复杂的方程我们没办法通过上述方法求取他的解,那么我们就得用迭代的方法来获取他的近似解,牛顿法就是一种常用的方法。

他的理论基础还是基于函数的泰勒近似:

对于方程式$f(x)=0$,我们先选择一个点$x_0$作为首轮迭代的初始点,在$x_0$处对$f(x)$进行一阶泰勒展开:

$f(x)\approx f(x_0)+(x-x_0)f'(x_0)$

大家知道,$f(x)$在$x=x_0$处的一阶泰勒展开$f(x_0)+(x-x_0)f'(x_0)$是函数$f(x)$在$x=x_0$处的切线,并以这条切线作为函数$f(x)$在$x=x_0$邻域的近似,那我们就通过这个直线方程$f(x_0)+(x-x_0)f'(x_0)=0$的根,来近似的作为$f(x)=0$的根。

具体获得的解如下:

$f(x_0)+(x-x_0)f'(x_0)=0\Rightarrow x=x_0-\frac{f(x_0)}{f'(x_0)}$

作为一条直线,获得这个切线方程的解是一件非常容易的事情,但是就这么求出来的解能够直接用来作为原方程$f(x)=0$的解吗?看看这误差,想想也不太敢直接用吧,下面这幅图显示的就一清二楚:
图1.用切线获得的近似解

图1.用切线获得的近似解

显然,利用一阶泰勒展开式求得的近似解$x_1=x_0-\frac{f(x_0)}{f'(x_0)}$,距离原方程$f(x)=0$真正的解距离还很远,但是他的的确确是比初始点$x_0$要更加接近于方程的真实解$x^*$,这也是一个好的趋势。

那么好,我们按照$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$这个方法继续迭代下去,就会产生一个牛顿法的近似解序列$x_1,x_{2},x_{3},...$,总的来说,他会越来越试图靠近方程的实际解析解,直到最后得到的两个相邻近似解之差小于预设的阈值$\epsilon$时(即此时满足$|x_{k+1}-x_k|<\epsilon$),这个迭代过程才停止,整个迭代的过程如下图所示:

图2.牛顿法解方程的迭代过程

图2.牛顿法解方程的迭代过程

从图中我们可以看到,牛顿法迭代出来的序列$x_0,x_1,x_2,x_3$正在不断的逼近方程实际的解:$x$。

1.2.实际举例

我们来实际看一个例子,我们用牛顿法来求方程$x^3-12x^2+8x+40=0$的一个根,我么选择$x_0=14$进行初始迭代:

代码片段:

import matplotlib.pyplot as plt
import numpy as np
import seaborn
seaborn.set()

def f(x):
    return x**3-12*x**2+8*x+40

def df(x):
    return 3*x**2-24*x+8

def newton(x):
    return x - f(x)/df(x)

x0 = 14
x_list = []
x_list.append(x0)

while(True):
    x = newton(x_list[-1])
    if abs(x-x_list[-1])<=1e-5:
        break
    x_list.append(x)

fig, ax = plt.subplots(2, 1)

x = np.linspace(0, 15, 1000)
plt.xlim(0, 15)
ax[0].plot(x, f(x))

x = np.linspace(10.5, 14.5, 1000)
plt.xlim(10.5, 14.5)
ax[1].plot(x, f(x))
ax[1].plot(x_list, [0]*len(x_list), 'ko')
print('x={},f(x)={}'.format(x_list[-1], f(x_list[-1])))
print(x_list)
plt.show()

运行结果:

x=10.933723301726161,f(x)=0.0003467730414286052
[14, 11.907692307692308, 11.079933151128879, 10.937805495069655, 10.933723301726161]

图3.用牛顿法解方程

图3.用牛顿法解方程

我们看到,在运行结果的两幅图中,上面的第一幅图是原函数$f(x)=x^3-12x^2+8x+40$的图像,下面的第二幅图是在$[10.5,14.5]$的小区间范围内的$f(x)$图像,$x$轴上的黑色点显示出:通过四轮迭代的解序列为
$[14, 11.9076, 11.0799, 10.9378, 10.9337]$。他们逐步逼近方程的一个根,最终收敛到一个方程的近似根,即:$x=10.9337$。经过验算,此时的函数取值$f(x)$已经相当接近于$0$了。这个过程,恰是计算机程序所擅长的。

不过需要注意的是,在利用牛顿法进行迭代求根的时候,初始值$x_0$的选取很重要,如果选取的不恰当,使得比值$\frac{f(x_0)}{f'(x_0)}$较大,会使得迭代的结果在真实根的附近不断摆动,可能会使得求解的过程失效。

2.利用牛顿法求一元函数极值

2.1.原理分析

到了这里,大家应该已经明白了如何利用牛顿法去迭代的求取一个方程的根,但是大家可能会有疑惑,我们不是在讨论如何求取函数的极值吗?这和解方程有什么关系?

这里面的关系非常紧密,我们接下来仔细分析一下。

试想,如果我们想要求取一个函数$f(x)$的极小值点,很显然可以利用之前介绍过的极值存在的条件,去寻找满足一阶导数$f'(x)=0$的$x$点,作为极小值的最终备选项,那么思路就很清晰了:

寻找使得$f(x)$取得极小值点$x$的问题,就被转化成了寻找满足一阶导数$f'(x)=0$的解方程的问题,这就又回到了利用牛顿法来解方程的问题了。

我们令$g(x)=f'(x)$,我们要求的就是方程$g(x)=0$的解,方法就完全一样了:也是先寻找一个合适的迭代初始点$x_0$,然后按照$x_{k+1}=x_k-\frac{g(x_k)}{g'(x_k)}$这个方法继续迭代下去,直到迭代出来的相邻两个解$x_k$和$x_{k+1}$之差小于阈值。

这里有一点需要转化一下:由于$g(x)=f'(x)$,因此实际操作中的迭代公式就变成了:$x_{k+1}=x_k-\frac{f'(x_k)}{f''(x_k)}$。

2.2.实际举例

我们来实际举一个例子,利用上述方法求取函数$f(x)=x^2-2sinx$的极值:

代码片段:
```
import matplotlib.pyplot as plt
import numpy as np
import seaborn
seaborn.set()

def f(x):
return x*2-2np.sin(x)

def df(x):
return 2x-2np.cos(x)

def d2f(x):
return 2+2*np.sin(x)

def newton(x):
return x - df(x)/d2f(x)

x0 = 3
x_list = []
x_list.append(x0)

top Created with Sketch.