积分逼近函数的递归

Question

我正在尝试使用自适应梯形规则来近似积分。

我有一个粗略的积分近似：

//Approximates the integral of f across the interval [a,b]
double coarse_app(double(*f)(double x), double a, double b) {
    return (b - a) * (f(a) + f(b)) / 2.0;
}

我有一个很好的积分近似：

//Approximates the integral of f across the interval [a,b]
double fine_app(double(*f)(double x), double a, double b) {
    double m = (a + b) / 2.0;
    return (b - a) / 4.0 * (f(a) + 2.0 * f(m) + f(b));
}

这是通过对给定间隔的递减部分的近似求和来实现自适应的，直到递归级别太高或粗略和精细近似彼此非常接近：

//Adaptively approximates the integral of f across the interval [a,b] with
//    tolerance tol.
double trap(double(*f)(double x), double a, double b, double tol) {
    double q = fine_app(f, a, b);
    double r = coarse_app(f, a, b);
    if ((currentLevel >= minLevel) && (abs(q - r) <= 3.0 * tol)) {
        return q;
    } else if (currentLevel >= maxLevel) {
        return q;
    } else {
        ++currentLevel;
        return (trap(f, a, b / 2.0, tol / 2.0) + trap(f, a + (b / 2.0), b, tol / 2.0));
    }
}

如果我通过将积分分解为多个部分并在其上使用fine_app 来手动计算积分，我会得到一个非常好的近似值。但是，当我使用应该为我执行此操作的陷阱功能时，我的所有结果都太小了。

例如，trap(square, 0, 2.0, 1.0e-2) 给出输出 0.0424107，其中 square 函数定义为 x^2。但是，输出应该在 2.667 左右。这比在整个时间间隔上执行一次fine_app 差得多，它的值是 3。

从概念上讲，我相信我已经正确实现了它，但是关于 C++ 递归的某些东西并没有达到我的预期。

第一次使用 C++ 编程，因此欢迎所有改进。

Answer 1

我假设您在其他地方定义了 currentLevel。你不想那样做。您还错误地计算了中点。

取 a = 3，b = 5：

[a, b / 2.0] = [3, 2.5]
[a + b / 2.0, b] = 2.5, 3]

正确的点应该是[3, 4]和[4, 5]

代码应如下所示：

double trap(double(*f)(double x), double a, double b, double tol, int currentLevel) {
    double q = fine_app(f, a, b);
    double r = coarse_app(f, a, b);
    if ((currentLevel >= minLevel) && (abs(q - r) <= 3.0 * tol)) {
        return q;
    } else if (currentLevel >= maxLevel) {
        return q;
    } else {
        ++currentLevel;
        return (trap(f, a, (a + b) / 2.0, tol / 2, currentLevel) + trap(f, (a + b) / 2.0, b, tol / 2, currentLevel));
    }
}

您可以添加一个辅助函数，这样您就不必指定 currentLevel：

 double integrate(double (*f)(double x), double a, double b, double tol)
 {
     return trap(f, a, b, tol, 1);
 }

如果我将其称为integrate(square, 0, 2, 0.01)，我会得到 2.6875 的答案，这意味着您需要更低的容差才能收敛到8/3 = 2.6666...7 的正确结果。您可以使用Simpson's method 的错误术语来检查确切的错误范围。