1.Intro
求函数的极值是科学计算中的基本问题。在介绍求函数极值的方法之前,我们先介绍一种求方程解的方法。
对于方程f(x)=0,满足如下条件的a,b称之为root,
a<b and f(a)>0,f(b)<0
我们可以用如下算法求方程的解:
1.令x=(a+b)/2,计算f(x),如果
a). f(x)>0,则令a=x
b). else,令b=x
2.重复1,直到b-a小于指定的精度,最后的x即为函数的一个根
类似的,求函数极值时我们可以找出满足如下条件的a,b,c三个值,称之为一个bracket
a < b < c and f (a) > f (b) and f (c) > f (b)
我们按照如下算法来确定函数的解。
1.给定初始的(a,b,c)
2.没有聚合(converge)之前:
a). 在a,b或b,c之前找一点x,计算f(x)
b). 根据f(x)和f(b)的关系,在a,b,c,x四点中找三点以满足上面的特性。
2.Note
- 不能只在a,b或b,c之间求x,否则可能陷入死循环,无法收敛。我们可以在比较大的区间上求x,以尽早收敛。
- 求x的方法有多种,比较直观的是取中点,黄金分割会稍微快一点。
- 这个算法只能算出一个极小值,不能算出最小值。
3.Code
下面是Python实现的黄金分割求极小值的代码,用sinx做测试。
#-*- coding:utf-8 -*-
import math
def minimise(func,a,b,c,tol):
GOLD=(3-math.sqrt(5))/2
if(a<b and b<c and func(a)>func(b) and func(b)<func(c)): #check wether it's valid bracket or not
while c-a>tol:
if b-a>c-b:
x=(1-GOLD)*a+GOLD*b # evalute x
if(func(x)>func(b)):
a=x
else:
c=b
b=x
else:
x=GOLD*b+(1-GOLD)*c # evalute x
if(func(x)>func(b)):
c=x
else:
a=b
b=x
return b,func(b) #return the value of x that minimises func(x) and minimum of func(x)
else:
raise Exception,'Invalid bracket.'
def func(x):
return math.sin(x)
if __name__=="__main__":
print minimise(func,0,math.pi+1,math.pi*2,0.1**10)
函数的输出为:
lonfee:~/tmp$ python minimise.py
(4.712388981257224, -1.0)