二分法求根和极小值的C++实现

#ifndef BSEARCH_H
#define BSEARCH_H

typedef double (*func) (double);

//二分法求根
double bsearch(double a, double b, double eps, func f);

//二分法求单峰极小值
double bsearch_min(double a, double b, double eps, func f, double &x);

#endif // BSEARCH_H

#include "bsearch.h"
#include "assert.h"
#include "math.h"
#include <QDebug>

//QT中的实现
//static inline bool qFuzzyIsNull(double d)
//{
//    return qAbs(d) <= 0.000000000001;
//}

inline int mysign(double x)
{
    return (x > 0.0);
}

double bsearch(double a, double b, double eps, func f)
{
    double fa = f(a);
    double fb = f(b);
    eps = fabs(eps);

    if (mysign(fa) == mysign(fb) && !qFuzzyIsNull(fa))
    {
        qDebug() << fa << fb;
        assert(false);
        return a;
    }

    if (fabs(fa) <= eps)
    {
        return a;
    }
    else if (fabs(fb) <= eps)
    {
        return b;
    }

    const static int max_tiems = 3000;
    int times = 0;

    double m;
    while (true)
    {
        m = (a+b)/2.0;
        //        qDebug() << "TEST: "
        //                 << QString::number(a, 'f', 15)
        //                 << QString::number(b, 'f', 15)
        //                 << QString::number(m, 'f', 15);

        if (fabs(a-b) <= qMin(fabs(a), fabs(b)) * eps || ++times >= max_tiems)
        {
            //qDebug() << "times1:" << times;
            break;
        }

        double fm = f(m);
        if (qFuzzyIsNull(fm))
        {
            //qDebug() << "times2:" << times;
            break;
        }

        if (mysign(fa) == mysign(fm))
        {
            a = m;
        }
        else
        {
            b = m;
        }
    }
    return m;
}


static func f1_f = NULL;
static double f1_eps = 1e-14;

//类似求导函数，写法比较原始，可以有更优雅的方式
double f1(double x)
{
    double eps2 = fabs(x) * f1_eps / 3; //步长要小于容差
    double y1 = f1_f(x);
    double y2 = f1_f(x + eps2);
    return (y2 - y1) / eps2;
}

//class func_base{
//public:
//  virtual double operator() (double) = 0;
//};

double bsearch_min(double a, double b, double eps, func f, double &x)
{
    f1_f = f;
    f1_eps = eps;

    x = bsearch(a, b, eps, f1);
    return f(x);
}

#include <QCoreApplication>
#include <QDebug>
#include <QtCore>
#include "bsearch.h"


int T = 0;

//测试求根
double myfunc(double x)
{
    T++;
    return pow(x, 3) + pow(x, 1)  - 100;
}

//测试求极小值
double myfunc2(double x)
{
    T++;
    return pow(x, 2)  -200 * x + 10000.0;
}

void test_root(double a, double b, double eps, func f)
{
    T = 0;
    QElapsedTimer timer;
    timer.start();

    double x,x2;

    int times = 1000;
    T = 0;
    timer.restart();
    for (int i = 0; i < times; i++)
    {
        x2 = bsearch(a, b, eps, f);
    }
    qDebug () << "bsearch "
              << "x="  << QString::number(x2, 'E', 14)
              << "fx=" << QString::number(f(x2) , 'E', 14)
              << timer.nsecsElapsed() / 1000000.0
              << "ms, times:" << T/times;
}


void test_min(double a, double b, double eps, func f)
{
    QElapsedTimer timer;
    timer.start();

    double x = 0, fx = 0;
    double x2 = 0, fx2 = 0;
    double x3 = 0, fx3 = 0;

    int times = 1000;


    T = 0;
    timer.restart();
    for (int i = 0; i < times; i++)
    {
        fx2 = bsearch_min(a, b, eps, f, x2);
    }
    qDebug () << "bsearch_min      "
              << "x="  << QString::number(x2, 'E', 14)
              << "fx=" << QString::number(fx2 , 'E', 14)
              << timer.nsecsElapsed() / 1000000.0
              << "ms, times:" << T/times;
}


int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);

    test_root(-1e100, 1e100, 1.0e-14, myfunc);
    test_min(-1e100, 1e100, 1.0e-14, myfunc2);

    return a.exec();
}

输出：

bsearch x= "4.56978016293266E+0" fx= "5.12617726045050E-13" 26.2534 ms, times: 368

bsearch_min x= "1.00000234750956E+2" fx= "5.51080114874480E-8" 42.0539 ms, times: 689