在C++中寻找根的二分法

Question

我写了一个简短的C/C++代码，通过二分法查找根。(这是一种简单的迭代数值方法，允许找到方程的根，即x，使f(x) = 0)。二分法

标题简单地由警卫和下列行组成：

#include <functional>
double bisection(double x1, 
                 double x2, 
                 double e, 
                 std::function<double(double)>& f);

为了得到简单的C代码，必须将函数的签名中的std::function<double(double)>转换为double(*f)(double)；

我有两个问题。

第一个问题是，这段代码是否可以改进(轻微的改进)。第二个问题是递归是否是一个好的选择。

namespace { 
    double bisection(double x1, 
                     double x2, 
                     double e, 
                     std::function<double(double)>& f, 
                     double fx1, 
                     double fx2){
        double mid = x1*0.5 + x2*0.5;
        if (x2 - x1 < e)
            return mid;
        double fmid = f(mid);
        if ((fmid>0.0) == (fx1 > 0.0))
            return bisection(mid, x2, e, f, fmid, fx2);
        return bisection(x1, mid, e, f, fx1, fmid);
    }
}

double bisection(double x1, 
                 double x2, 
                 double e, 
                 std::function<double(double)>& f){
    double fx1 = f(x1);
    double fx2 = f(x2);
    return bisection(x1, x2, e, f, fx1, fx2);
}

main使用我的代码的一个例子是：

#include "bisection.hpp"
int main(){
    bisection(0.0, 1.0, 1.0E-6, [](double x){return x*x-0.5;})

    return 0;
}       

请注意，预计x1和x2将被订购: x1 < x2。此外，用户应该意识到存在一个log_2((x2-x1)/tolerance)调用堆栈。

Answer 1

• 在标准C++中，临时不能绑定到非const引用(由于编译器扩展，MSVC++将编译它)：使f成为const引用
• Const正确性:尽可能使用const
• 可以避免一些临时变量。例如，fmid
• 你做x1 * 0.5 + x2 * 0.5而不是(x1 + x2) * 0.5有什么原因吗？如果出于一些数字上的原因，那就好了。否则，我发现后者更直观。
• 你不需要fx2
• 递归确实更昂贵(如果不是编译器优化的话)和更难阅读。

考虑到以上所有考虑，下面是一个简单的非递归版本：

double bisection(
    double x1,
    double x2,
    const double e,
    const std::function<double(double)>& f)
{
    while(x2 - x1 >= e)
    {
        const double mid = (x1 + x2) * 0.5;
        if((f(mid) > 0.0) == (f(x1) > 0.0))
            x1 = mid;
        else
            x2 = mid;
    }
    return (x1 + x2) * 0.5;
}

下面是缓存f(x1)的版本：

double bisection(
    double x1,
    double x2,
    const double e,
    const std::function<double(double)>& f)
{
    double fx1 = f(x1);
    while(x2 - x1 >= e)
    {
        const double mid = (x1 + x2) * 0.5;
        const double fmid = f(mid);
        if((fmid > 0.0) == (fx1 > 0.0))
        {
            x1 = mid;
            fx1 = fmid;
        }
        else
        {
            x2 = mid;
        }
    }
    return (x1 + x2) * 0.5;
}

Answer 2

除了来自其他答案的建议之外，惯用的C++将为函数对象类型使用模板参数，而不是std::function。无论出于什么原因，当您需要在运行时后期绑定函数时，std::function是很棒的，但它也是一个巨大的优化障碍。我想这会快得多：

template <typename Func>
double bisection(double x1, double x2, double e, Func f)
{
    double fx1 = f(x1);
    double mid = x1 * 0.5 + x2 * 0.5;
    while (x2 - x1 >= e) {
        double fmid = f(mid);
        if ((fx1 > 0.0) == (fmid > 0.0)) {
            x1 = mid;
            fx1 = fmid;
        } else {
            x2 = mid;
        }
        mid = x1 * 0.5 + x2 * 0.5;
    }
    return mid;
}

(大部分实现都归功于AMA的回答。)

您甚至可以通过参数化值类型，而不是假设double，使其更具一般性。