我需要整合一个函数(两个变量)。我知道我可以通过使用Fubini 定理来整合一个变量函数,然后使用数值方法,例如Rectangle 方法或Trapezoidal rule。
但是在C++中是否有任何预构建的函数可以做到这一点?我需要对单位R2
三角形进行积分((0,0), (1,0), (0,1))
。
我需要整合一个函数(两个变量)。我知道我可以通过使用Fubini 定理来整合一个变量函数,然后使用数值方法,例如Rectangle 方法或Trapezoidal rule。
但是在C++中是否有任何预构建的函数可以做到这一点?我需要对单位R2
三角形进行积分((0,0), (1,0), (0,1))
。
您可以使用GNU Scientific Library,它支持许多“数值分析”功能,包括集成。
手册中一个非常简单的集成示例只是几行代码:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
double f (double x, void * params) {
double alpha = *(double *) params;
return log(alpha*x) / sqrt(x);
}
int
main (void)
{
double result, error;
double expected = -4.0;
double alpha = 1.0;
gsl_integration_workspace * w
= gsl_integration_workspace_alloc (1000);
gsl_function F;
F.function = &f;
F.params = α
gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000,
w, &result, &error);
printf ("result = % .18f\n", result);
printf ("exact result = % .18f\n", expected);
printf ("estimated error = % .18f\n", error);
printf ("actual error = % .18f\n", result - expected);
printf ("intervals = %d\n", w->size);
gsl_integration_workspace_free (w);
return 0;
}
据我所知,标准库中没有您要搜索的类型的函数。这是一种实现:
f(x)
对于要在固定极限a
和之间积分的固定函数,可以b
将扩展梯形规则中的区间数加倍,而不会失去先前工作的好处。梯形规则最粗略的实现是在函数的端点a
和处平均b
。细化的第一阶段是将函数在中间点的值添加到该平均值中。细化的第二阶段是添加1/4
和3/4
点的值。
许多基本正交算法涉及添加连续的细化阶段。将这个特性封装在一个Quadrature
结构中很方便:
struct Quadrature
{
//Abstract base class for elementary quadrature algorithms.
Int n; // Current level of refinement.
virtual Doub next() = 0;
//Returns the value of the integral at the nth stage of refinement.
//The function next() must be defined in the derived class.
};
然后Trapzd
结构由此衍生如下:
template<class T>
struct Trapzd: Quadrature
{
Doub a, b, s; // Limits of integration and current value of integral.
T &func;
Trapzd() { };
// func is function or functor to be integrated between limits: a and b
Trapzd(T &funcc, const Doub aa, const Doub bb)
: func(funcc), a(aa), b(bb)
{
n = 0;
}
// Returns the nth stage of refinement of the extended trapezoidal rule.
// On the first call (n = 1), the routine returns the crudest estimate
// of integral of f x / dx in [a,b]. Subsequent calls set n=2,3,... and
// improve the accuracy by adding 2n - 2 additional interior points.
Doub next()
{
Doub x, tnm, sum, del;
Int it, j;
n++;
if (n == 1)
{
return (s = 0.5 * (b-a) * (func(a) + func(b)));
}
else
{
for (it = 1, j = 1; j < n - 1; j++)
{
it <<= 1;
}
tnm = it;
// This is the spacing of the points to be added.
del = (b - a) / tnm;
x = a + 0.5 * del;
for (sum = 0.0,j = 0; j < it; j++, x += del)
{
sum += func(x);
}
// This replaces s by its refined value.
s = 0.5 * (s + (b - a) * sum / tnm);
return s;
}
}
};
该Trapzd
结构可以以多种方式使用。最简单粗暴的方法是通过扩展梯形规则集成一个函数,您可以在其中预先知道您想要的步数。如果你愿意2^M + 1
,你可以通过片段来完成:
Ftor func; // Functor func here has no parameters.
Trapzd<Ftor> s(func, a, b);
for(j = 1 ;j <= m + 1; j++) val = s.next();
答案返回为val
. 这Ftor
是一个包含要集成的函数的函子。
当然,更好的方法是细化梯形规则,直到达到特定的精确度。一个功能是:
template<class T>
Doub qtrap(T &func, const Doub a, const Doub b, const Doub eps = 1.0e-10)
{
// Returns the integral of the function or functor func from a to b.
// The constants EPS can be set to the desired fractional accuracy and
// JMAX so that 2 to the power JMAX-1 is the maximum allowed number of
// steps. Integration is performed by the trapezoidal rule.
const Int JMAX = 20;
Doub s, olds = 0.0; // Initial value of olds is arbitrary.
Trapzd<T> t(func, a, b);
for (Int j = 0; j < JMAX; j++)
{
s = t.next();
if (j > 5) // Avoid spurious early convergence.
{
if (abs(s - olds) < eps * abs(olds) || (s == 0.0 && olds == 0.0))
{
return s;
}
}
olds = s;
}
throw("Too many steps in routine qtrap");
}
typedef double Doub; // 64 - bit floating point
typedef int Int; // 32 - bit signed integer
您可能希望查看 Boost 正交和微分库。特别是,它们提供了梯形规则的一个版本:
https://www.boost.org/doc/libs/1_69_0/libs/math/doc/html/math_toolkit/trapezoidal.html
Quadrature/Differentiation Library 编写得非常好,并且与现代 C++ 兼容,因为可以只使用 lambda 表达式或函数对象作为被积函数。我很快就开始使用它。
这是一个逼近 pi 的示例,通过逼近 4/(1 + x^2) 的积分,从 x = 0 到 x = 1,将被积函数作为 lambda 表达式。
#include <boost/math/quadrature/trapezoidal.hpp>
#include <iostream>
using boost::math::quadrature::trapezoidal;
using std::cout;
using std::endl;
// Put inside a test function:
auto f = [](double x)
{
return 4.0 / (1.0 + x * x);
};
double appPi = trapezoidal(f, 0.0, 1.0);
double tol = 1e-6;
int max_refinements = 20;
double appPi2 = trapezoidal(f, 0.0, 1.0, tol, max_refinements);
cout << "Trapezoid Rule results for computing pi by integration" << endl;
cout << "a) with defaults, and b) with tol and max_refine set : " << endl;
cout << appPi << ", " << appPi2 << endl << endl;
我提供了两个示例,一个使用默认设置来离散化积分和收敛范围,第二个使用自定义设置。
我的结果(只是从屏幕上获取输出的副本)是:
Trapezoid Rule results for computing pi by integration
a) with defaults, and b) with tol and max_refine set :
3.14159, 3.14159