C++ 中的数值积分
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/16512817/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Numerical integration in C++
提问by user2370139
I need to integrate a function (of two variables). I know I can do it by using Fubini theoremto integrate one variable functions, then using numerical methods such as the Rectangle methodor the Trapezoidal rule.
我需要集成一个函数（两个变量）。我知道我可以通过使用Fubini 定理对一个变量函数进行积分，然后使用诸如矩形方法或梯形规则之类的数值方法来做到这一点。
But are there any pre-built functions to do that in C++? I need to integrate over the unit R2triangle ((0,0), (1,0), (0,1)).
但是在C++ 中是否有任何预先构建的函数可以做到这一点？我需要对单位R2三角形进行积分((0,0), (1,0), (0,1))。
回答by Wagner Patriota
You can use the GNU Scientific Library, which supports many "Numerical analysis" functions including integration.
您可以使用GNU Scientific Library，它支持许多“数值分析”功能，包括积分。
A very simple example of integrationfrom the manual is just a few lines of code:
手册中一个非常简单的集成示例只是几行代码：
#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;
}
回答by Ziezi
As far as I know there are no functions of the type you are searching for in the standard library. Here is one implementation of the:
据我所知，标准库中没有您要搜索的类型的函数。下面是一个实现：
The extended trapezoidal rule.
扩展梯形规则。
For a fixed function
f(x)to be integrated between fixed limitsaandb, one can double the number of intervals in the extended trapezoidal rule without losing the benefit of previous work. The coarsest implementation of the trapezoidal rule is to average the function at its endpointsaandb. The first stage of refinement is to add to this average the value of the function at the halfway point. The second stage of refinement is to add the values at the1/4and3/4points.
对于
f(x)要在固定极限a和之间积分的固定函数，可以将b扩展梯形规则中的间隔数加倍，而不会失去先前工作的好处。梯形规则的最粗略的实现是在函数的端点a和处求平均b。细化的第一阶段是将函数在中间点的值添加到该平均值中。细化的第二阶段是在1/4和3/4点添加值。
A number of elementary quadrature algorithms involve adding
successive stages of refinement. It is convenient to encapsulate this feature in a Quadraturestructure:
许多基本正交算法涉及添加连续的细化阶段。将此功能封装在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.
};
Then the Trapzdstructure is derived from this as follows:
那么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;
}
}
};
Usage:
用法：
The Trapzdstructure could be used in several ways. The
simplest and crudest is to integrate a function by the extended trapezoidal rule where you know in advance the number of steps you want. If you
want 2^M + 1, you can accomplish this by the fragment:
该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();
with the answer returned as val. Here Ftoris a functor containing the function to be integrated.
返回的答案为val。这Ftor是一个包含要集成的函数的函子。
Bonus:
奖金：
Much better, of course, is to refine the trapezoidal rule until some specified degree of accuracy has been achieved. A function for this is:
当然，更好的是细化梯形规则，直到达到某种指定的准确度。一个函数是：
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");
}
Typedefs.
类型定义。
typedef double Doub; // 64 - bit floating point
typedef int Int; // 32 - bit signed integer
回答by duffymo
I would recommend Gaussian quadrature and linear shape functions:
我会推荐高斯正交和线性形状函数：
http://www.mathworks.com/matlabcentral/fileexchange/9230-gaussian-quadrature-for-triangles
http://www.mathworks.com/matlabcentral/fileexchange/9230-gaussian-quadrature-for-triangles
http://www.wolframalpha.com/input/?i=gaussian+quadrature+triangle
http://www.wolframalpha.com/input/?i=gaussian+quadrature+triangle
回答by djhanson
You might wish to check out the Boost Quadrature and Differentiation Library. In particular, they provide a version of the Trapezoid Rule:
您可能希望查看升压正交和微分库。特别是，他们提供了梯形规则的一个版本：
https://www.boost.org/doc/libs/1_69_0/libs/math/doc/html/math_toolkit/trapezoidal.html
https://www.boost.org/doc/libs/1_69_0/libs/math/doc/html/math_toolkit/trapezoidal.html
The Quadrature/Differentiation Library is very well written and compatible with modern C++, in that one can just use a lambda expression or a function object for the integrand. I was up and running with it very quickly.
正交/微分库编写得非常好并且与现代 C++ 兼容，因为可以只使用 lambda 表达式或函数对象作为被积函数。我很快就开始使用它。
Here's an example of approximating pi, by approximating the integral of 4/(1 + x^2), from x = 0 to x = 1, putting the integrand as a lambda expression.
这是一个近似 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;
I provided two examples, one using the default settings for the discretization of the range of integration and convergence, and the second using custom settings.
我提供了两个示例，一个使用用于积分和收敛范围离散化的默认设置，第二个使用自定义设置。
My results (just taking a copy of the output from the screen) are:
我的结果（只是从屏幕上复制输出）是：
Trapezoid Rule results for computing pi by integration
a) with defaults, and b) with tol and max_refine set :
3.14159, 3.14159
