黄金分割搜索法求单峰极小值C++实现

#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