#ifndef GOLDEN_SEARCH_H
#define GOLDEN_SEARCH_H
/**
* @brief golden_min 黄金分割搜索,求最小值;要求f(x)在[a,b]区间内有单峰最小值
* 每次问题规模减少(1-GR)=0.381966011250
* @param a 函数f的自变量x的上边界
* @param b 函数f的自变量x的上边界
* @param eps 自变量容差
* @param f f(x)函数
* @param x 求得最小值时的x
* @return 最小值
*/
double golden_min(double a, double b, double eps, double f ( double x ), double &x);
#endif // GOLDEN_SEARCH_H
#include "golden_search.h"
#include "math.h"
#include "assert.h"
#include <QtCore>
double golden_min(double a, double b, double eps, double f(double), double &x)
{
//a----c--d----b
//ac为小头,bc为大头, ac=db, bc/ab为黄金分割比例
double result = NAN;
const double GR = (sqrt(5.0) - 1.0)/2.0; //黄金分割比例, golden ratio
double c = a + (b-a) * (1 - GR);
double fc = f(c);
bool bTure = (fc <= f(a) && fc <= f(b));
if (!bTure)
{
//初始参数不满足f(c) <= f(a)和f(b)
assert(false);
return result;
}
//保护措施,防止函数不满足条件等导致的死循环
const int max_times = 3000;
int times = 0;
bool bSwapCD = false;
double d, fd;
while (true)
{
//安全保护
times++;
if (times > max_times)
{
result = fc;
break;
}
c = a + (b-a) * (1 - GR);
d = a + (b-a) * GR;
if (bSwapCD)
{
fc = f(c);
}
else
{
fd = f(d);
}
if (c > d)
{
assert(false);
}
if (fabs(a - b) < eps * qMin(fabs(a), fabs(b)))
{
x = (a + b) / 2;
result = f(x);
//std::cout << "times:" << times;
break;
}
if (fc <= fd)
{
b = d;
fd = fc;
bSwapCD = true;
}
else
{
a = c;
fc = fd;
bSwapCD = false;
}
}
return result;
}
#include <QCoreApplication>
#include <QtCore>
#include "golden_search.h"
static int T = 0;
typedef double (*func) (double);
static double myfunc2(double x)
{
T++;
return pow(x, 2) -200 * x + 10000.0;
}
static void test_min(double a, double b, double eps, func f)
{
QElapsedTimer timer;
timer.start();
double x3 = 0, fx3 = 0;
int times = 1000;
T = 0;
timer.restart();
for (int i = 0; i < times; i++)
{
fx3 = golden_min(a, b, eps, f, x3);
}
qDebug () << "golden_min "
<< "x=" << QString::number(x3, 'E', 14)
<< "fx=" << QString::number(fx3 , 'E', 14)
<< timer.nsecsElapsed() / 1000000.0
<< "ms, times:" << T/times;
}
int main(int argc, char *argv[])
{
QCoreApplication a(argc, argv);
test_min(-1e100, 1e100, 1.0e-14, myfunc2);
return a.exec();
}
运行输出:
golden_min x= "9.99999999789266E+1" fx= "0.00000000000000E+0" 22.8964 ms, times: 543