python实现黄金分割搜索算法+动态展示

前言

数值算法是跟数学关系比较密切的一门课程，主要是用计算机程序实现一些数学公式 和数学算法。我觉得计算机处理数据的方式与数学最大的不同就是它是离散化，虽然计算机处理数据具有一定的误差性，但是由于机器的精度远远大于实际工程中应用的误差，所以用计算机计算出的数值大多数是可以利用的数值。
离散化是计算机设计和算法中一种重要的思想。

要求

用黄金分割搜索算法求cos(x),x∈[-π/2,π/2]的最大值，设计出具体的程序，使之能够动态演示搜索过程。

黄金分割搜索算法原理

设一条线段AB的长度为a，C点在靠近B点的黄金分割点上，且AC为b，则a比b就是黄金数。

以上是计算的过程（来自百度）。
二分之根号五减一就是0.618。给我们一个区间乘以0.618就是他的黄金分割点。
那么咋么利用黄金分割点求单峰最大值呢？
求最大值的题目一般会给一个区间，我们这里是[-π/2,π/2]。想求得其中的最大值，需要在这个区间中选择两个点，计算它们的函数值，使其中的一个点替换边界的一个点，根据内部点选定规则再计算新的内部点，然后不断的计算直至两个内部点计算出的函数值的差小于规定的误差，取这两个函数值的平均就是 我们要求的最大值。
黄金分割 就是在此基础之上定义了两个内部点的选定规则。
在搜索范围[a,b]中选择两点，r1=a+0.382(b-a), r2=a+0.618(b-a)。计算f(a),f(r1),f(r2),f(b)的函数值，为了使下一次计算值接近极值，我们选择用r1或者r2替换a或者b，根据所求最大还是最小情况有所不同，如下图中，f(r2)>f(r1),最小值更靠近f(r1)，舍弃b的值，使b=r2，在计算新的r2，如此反复计算，直到误差小于要求值。

（图源：网络）

伪代码

1: Input（a,b,e）
2: x1=a+(b-a)0.382
3: x2=a+(b-a)0.618
4: f1=f(x1)
5: f2=f(x2)
6: while(b-a)>e{
7: if(f2>f1)
8: a=x1;
9: x1=x2;
10: f1=f2;
11: x2=a+0.618(b-a);
12: f2=f(x2);
13：else
14： b=x2;
15： x2=x1;
16: f2=f1;
17: x1=a+0.382(b-a);
18: f1=f(x1);
19: x max=(a+b)/2;
20: output(xmax)

代码编写

要求是要动态展示搜索过程，我觉得python的图形库还是很好用的就用了python的库matplotlib。用matlab也可以实现。但由于对matlab不是很熟就用了python。
要动态实现就要用animation，这也是python提供好的一个动态展示的方法。但是有一些局限性。
numpy库一般用来进行数据处理。也是很强大的一个库。
想要保存动图还要下载ImageMagick Display

动态结果

python代码

import numpy as np 
import matplotlib.pyplot as plt
from matplotlib import animation

def f(x):
    y=np.cos(x)
    return y
#创建一个画板和画布
fig, ax = plt.subplots()
x = np.linspace(-np.pi/2, np.pi/2, 200)
y = f(x)
#创建坐标轴
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data',0))
l = ax.plot(x, y)
dot, = ax.plot([], [], 'r.')
xdata, ydata =  [], []

#初始化函数	
def init():
    ax.set_xlim(-np.pi/2, np.pi/2)
    ax.set_ylim(-1, 2)
    return l
#生成节点函数
def gen_dot():
    tdata,ldata= [], []
    xdata=[-np.pi/2, np.pi/2]
    ydata=f(xdata)
    ydata=ydata.tolist()
    a=xdata[0]
    b=xdata[1]
    x1=a+0.382*(b-a)
    x2=a+0.618*(b-a)
    e=0.0005
    f1=f(x1)
    f2=f(x2)
    xdata.append(x1)
    xdata.append(x2)
    ydata.append(f1)
    ydata.append(f2)
    while(b-a>e):
        if(f2>f1):
            a=x1
            x1=x2
            f1=f2
            x2=a+0.382*(b-a)
            f2=f(x2)
            xdata.append(x2)
            ydata.append(f2)
            
        else:
            b=x2
            x2=x1
            f2=f1
            x1=a+0.618*(b-a)
            f1=f(x1)
            xdata.append(x1)
            ydata.append(f1)
    l=(a+b)/2
    xdata.append(l)
    ydata.append(f(l))
    print(l)
    print(f(l))
    for i in range(len(xdata)):
        tdata.append(xdata[i])
        ldata.append(ydata[i])
        newdot=[tdata,ldata]
        yield newdot
    xdata, ydata ,tdata,ldata= [], [], [], []
#更新节点函数
def update_dot(newd):
    dot.set_data(newd[0], newd[1])
    return dot,
#生成动态图像的节点
ani = animation.FuncAnimation(fig, update_dot, frames = gen_dot, interval = 1000, init_func=init)
ani.save('sin_dot.gif', writer='imagemagick', fps=30)

plt.show()