函数作为函数的参数，使用二分法C++

Function as parameter of a function, using bisection method C++

本文关键字：函数 二分 C++ 参数      更新时间：2023-10-16
 // Function as parameter of a function
#include <iostream>
#include <cmath>
#include <cassert>
using namespace std;
const double PI = 4 * atan(1.0); // tan^(-1)(1) == pi/4 then 4*(pi/4)== pi
typedef double(*FD2D)(double);
double root(FD2D, double, double); //abscissae of 2-points, 
//For the algorithm to work, the function must assume values of opposite sign in these two points, check
// at point 8
double polyn(double x) { return 3 - x*(1 + x*(27 - x * 9)); }
int main() {
double r;
cout.precision(15);
r = root(sin, 3, 4);
cout << "sin:       " << r << endl
    << "exactly:  " << PI << endl << endl;
r = root(cos, -2, -1.5);
cout << "cos:       " << r << endl
    << "exactly:  " << -PI/2 << endl << endl;
r = root(polyn, 0, 1);
cout << "polyn:       " << r << endl
    << "exactly:  " << 1./3 << endl << endl;
/* 
we look for the root of the function equivalent to function polyn
but this time defined as a lambda function 
*/
r = root([](double x) -> double { 
            return 3 - x*(1 + x*(27 - x * 9)); 
        }, 0, 1);
cout << "lambda:       " << r << endl
    << "exactly:  " << 1. / 3 << endl << endl;
return 0;
}
// Finding root of function using bisection.
// fun(a) and fun(b)  must be of opposite sign
double root(FD2D fun, double a, double b) {
static const double EPS = 1e-15; // 1×10^(-15)
double f, s, h = b - a, f1 = fun(a), f2 = fun(b);
if (f1 == 0) return a;
if (f2 == 0) return b;
assert(f1*f2<0); // 8.- macro assert from header cassert.
do {
    if ((f = fun((s = (a + b) / 2))) == 0) break;
    if (f1*f2 < 0) { 
        f2 = f; 
        b = s;
    }
    else {
        f1 = f; 
        a = s; 
    }
} while ((h /= 2) > EPS);
    return (a + b) / 2;
}

有人可以解释一下双根函数中的循环是如何工作的吗？我似乎不是 100%，我在网上检查了这个二分法并在纸上尝试，但我无法从这个例子中弄清楚。提前感谢！

如果将"聪明"行拆分为多行，则更容易理解。以下是一些修改，以及评论：

double root(FD2D fun, double a, double b) {
    static const double EPS = 1e-15; // 1×10^(-15)
    double fa = fun(a), fb = fun(b);
    // if either f(a) or f(b) are the root, return that
    // nothing else to do
    if (fa == 0) return a;
    if (fb == 0) return b;
    // this method only works if the signs of f(a) and f(b)
    // are different. so just assert that
    assert(fa * fb < 0); // 8.- macro assert from header cassert.

    do {
        // calculate fun at the midpoint of a,b
        // if that's the root, we're done
        // this line is awful, never write code like this...
        //if ((f = fun((s = (a + b) / 2))) == 0) break;
        // prefer:
        double midpt = (a + b) / 2;
        double fmid = fun(midpt);
        if (fmid == 0) return midpt;
        // adjust our bounds to either [a,midpt] or [midpt,b]
        // based on where fmid ends up being. I'm pretty
        // sure the code in the question is wrong, so I fixed it
        if (fa * fmid < 0) { // fmid, not f1
            fb = fmid; 
            b = midpt;
        }
        else {
            fa = fmid; 
            a = midpt; 
        }
    } while (b-a > EPS); // only loop while
                         // a and b are sufficiently far
                         // apart
    return (a + b) / 2;  // approximation
}
相关文章：