Fibonacci法与黄金分割法

Fibonacci法与黄金分割法均适用于求解单峰函数的极小值问题。

单峰函数

定义：

设f(x)f(x)f(x)是定义在某实数集合D上的实值函数。

如果存在x0∈Dx_0\in Dx0∈D，对任意的x1,x2∈Dx_1,x_2 \in Dx1,x2∈D，

当x1≤x2≤x0x_1 \leq x_2 \leq x_0x1≤x2≤x0时，有f(x1)≤f(x2)f(x_1)\leq f(x_2)f(x1)≤f(x2)

当x0≤x1≤x2x_0\leq x_1 \leq x_2x0≤x1≤x2时，有f(x1)≥f(x2)f(x_1)\geq f(x_2)f(x1)≥f(x2)

则称f(x)f(x)f(x)为单峰函数（如下图所示，函数图像只有一个峰）

如果在下单峰函数的定义域[a,b][a,b][a,b]中，任取两点a1、b1，a1<b1a_1、b_1，a_1 <b_1a1、b1，a1<b1，对于该函数的极小值点位置可能出现两种情况：

f(a1)<f(b1)f(a_1) < f(b_1)f(a1)<f(b1)，则极小值点xxx位于[a,b1][a,b_1][a,b1]
f(a1)≥f(b1)f(a_1) \geq f(b_1)f(a1)≥f(b1)，则极小值点x位于[a1,b][a_1, b][a1,b]

通过这种方式，我们就可以将含有极小值点的区间不断缩短，当缩短到一定程度时，就可以将区间端点作为极小值点的近似值。

Fibonacci法和黄金分割法就是通过这种思想逐步缩短区间的长度，从而得到极小值点的近似值

Fibonacci法

Fibonacci法通过Fibonacci数列得到每一次迭代的区间缩短率。

算法描述

具体的算法步骤如下：

输入区间[a,b][a,b][a,b]，相对精度ϵ\epsilonϵ
根据Fn≥1ϵF_n\geq \frac{1}{\epsilon}Fn≥ϵ1，得到最少的试点个数nnn
令k=0k=0k=0
根据Fibonacci数列确定试点的位置，t1=a+Fn−kFn−k+1(a−b),t2=a+Fn−kFn−k+1(b−a)t_1=a+\frac{F_{n-k}}{F_{n-k+1}}(a-b),t_2=a+\frac{F_{n-k}}{F_{n-k+1}}(b-a)t1=a+Fn−k+1Fn−k(a−b),t2=a+Fn−k+1Fn−k(b−a)
比较两试点的函数值大小。若f(t1)<f(t2)f(t_1)<f(t_2)f(t1)<f(t2)，令a=a,b=t2,t2=t1a=a,b=t_2,t_2=t_1a=a,b=t2,t2=t1;若f(t1)≥f(t2)f(t_1)\geq f(t_2)f(t1)≥f(t2)，令a=t1,b=b,t1=t2a=t_1,b=b,t_1=t_2a=t1,b=b,t1=t2；k=k+1k=k+1k=k+1，转至第4步，直到k=n−1k=n-1k=n−1，转第6步
当k=n−1k=n-1k=n−1时，有t1=t2=12(a+b)t_1=t_2=\frac{1}{2}(a+b)t1=t2=21(a+b)，此时无法根据f(t1)、f(t2)f(t_1)、f(t_2)f(t1)、f(t2)的大小确定极小值点，为此，取
t1=12(a+b)t2=a+(12+δ)(b−a)，常取δ=0.01t_{1}=\frac{1}{2}(a+b) \\ t_2 = a + (\frac{1}{2}+\delta)(b-a)，常取\delta=0.01 \\ t1=21(a+b)t2=a+(21+δ)(b−a)，常取δ=0.01
比较f(t1)、f(t2)f(t_1)、f(t_2)f(t1)、f(t2)的值，将较小者作为函数极小值的近似值，对应的点作为近似的极小值点。

算法结束。
注意，虽然算法迭代了n次，但使用的试点数并不是2n个，而是n个。因为在每一迭代中，有一个试点是上一步迭代中得到的。

具体实例及代码实现

对于函数f(t)=et+e−t,t∈[−1,1]f(t)=e^t+e^{-t},t\in [-1,1]f(t)=et+e−t,t∈[−1,1]。其函数图像如下：

代码实现如下：

import math
"""
Fibonacci法求解凸函数极小值
"""
class Fibonacci():def __init__(self, a0, b0, delta, epsilon):""":param a0: 区间左端点:param b0: 区间右端点:param delta: 相对精度:param epsilon: 误差"""self.delta = deltaself.a0 = a0self.b0 = b0self.epsilon = epsilondef Fib(self, n):"""生成长度为n的Fibonacci数列:param n: Fibonacci数列长度:return: Fibonacci数列"""result = []result.append(0)result.append(1)result.append(1)for i in range(3, n + 1):result.append(result[i - 1] + result[i - 2])return resultdef get_n(self, delta):"""根据给定的相对精度返回所需Fibonacci数列的长度:param delta: 相对精度:return: Fibonacci数列的长度"""F = self.Fib(1000)for i in range(1, len(F)):if F[i] > 1 / delta or F[i] == 1 / delta:print("试点数为{0}".format(i))return iprint("the length of fib is low")def fibonacci(self):"""求解凸函数的极小值点以及极小值，并检验是否满足相对精度:return: 函数的极小值点和极小值，以及检验结果"""a = self.a0b = self.b0n = self.get_n(delta)F = self.Fib(n)for i in range(1, n - 1):t1 = b + (F[n - i] / F[n - i + 1]) * (a - b)t2 = a + (F[n - i] / F[n - i + 1]) * (b - a)if f(t1) < f(t2):b = t2t2 = t1else:a = t1t1 = t2t1 = 0.5 * (a + b)t2 = a + (0.5 + epsilon) * (b - a)if f(t1) < f(t2):b = t2print("最优点为{0},最优值为{1}".format(t1, f(t1)))print("检验结果为{0}".format((b - a) / (self.b0 - self.a0) < self.delta))else:a = t1print("最优点为{0},最优值为{1}".format(t2, f(t2)))print("检验结果为{0}".format((b - a) / (self.b0 - self.a0) < self.delta))# 待求解的下单峰函数
def f(t):""":param t: 自变量:return: 待求解函数"""return math.exp(t) + math.exp(-t)# return -1 * ((math.sin(t) ** 6) * (math.tan(1 - t) * math.exp(30 * t)))if __name__ == "__main__":# 实例演练delta = 0.00001epsilon = 0.01a0 = 0b0 = 1fib = Fibonacci(a0, b0, delta, epsilon)fib.fibonacci()

求解结果为：

黄金分割法

黄金分割法的思路与Fibonacci法类似，只不过Fibonacci法的区间缩短率是通过Fibonacci数列来得到。而黄金分割法的区间缩短率是固定的，为0.618.

算法描述

具体实例及代码实现

对于下单峰函数f(t)=sin(t)6tan(1−t)e30t,t∈[0,1]f(t)=sin(t)^6tan(1-t)e^{30t},t\in [0,1]f(t)=sin(t)6tan(1−t)e30t,t∈[0,1]；其函数图像为：

具体代码如下：

import math
"""
黄金分割法法求解凸函数极小值
"""
class GoldenSection():def __init__(self, a0, b0, delta, epsilon):"""初始化参数:param a0: 区间左端点:param b0: 区间右端点:param delta: 相对精度:param epsilon:"""self.delta = deltaself.a0 = a0self.b0 = b0self.epsilon = epsilondef solve(self):"""求解凸函数的 满足相对精度 极小值点以及极小值:return: 函数的极小值点和极小值"""a = self.a0b = self.b0# n = self.get_n(delta)while True:if (b - a) / (self.b0 - self.a0) < self.delta:breakt1 = a + 0.382 * (b - a)t2 = a + 0.618 * (b - a)if f(t1) < f(t2):b = t2else:a = t1t1 = 0.5 * (a + b)t2 = a + (0.5 + epsilon) * (b - a)if f(t1) < f(t2):print("最优点为{0},最优值为{1}".format(t1, f(t1)))else:print("最优点为{0},最优值为{1}".format(t2, f(t2)))# 待求解的上单峰函数
def f(t):""":param t: 自变量:return: 待求解函数"""#return math.exp(t) + math.exp(-t)return -1 * ((math.sin(t) ** 6) * (math.tan(1 - t) * math.exp(30 * t))) # 乘以-1，转换为下单峰函数if __name__ == "__main__":# 实例演练delta = 0.00001a0 = -1b0 = 1epsilon = 0.01fib = GoldenSection(a0, b0, delta, epsilon)fib.solve()

求解结果为：