我正在尝试使用自适应梯形规则近似积分。
我有一个粗略的积分近似值:
//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函数时,我的所有结果都太小了。
例如,trap(square,0,2.0,1.0e-2)给出输出0.0424107,其中square函数定义为x ^ 2。但是,输出应约为2.667。这远比在整个时间间隔内运行一次fine_app(值3)要糟糕得多。
从概念上讲,我相信我已经正确实现了它,但是有关C++递归的某些事情并没有达到我的期望。
首次使用C++进行编程,因此欢迎所有改进。
最佳答案
我假设您在其他地方定义了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的错误条款检查与此相关的确切错误。关于c++ - 积分近似函数的递归,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/7705034/