All geek questions in one place

Numerical integration in C ++

I need to integrate a function (from two variables). I know that I can do this using Fubini's theorem to integrate a single variable of a function, and then using numerical methods such as the Rectangle method or the Trapezoid rule.

But are there any pre-created functions for this in C ++? I need to integrate a triangle over the unit R2 ((0,0), (1,0), (0,1)) .

+4
source share
3 answers

You can use the GNU Science Library , which supports many of the functions of Numerical Analysis, including integration.

A very simple example of integration from 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 = &alpha; 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; }
+6
source

As far as I know, there are no functions of the type that you are looking for in the standard library. Here is one of the implementations:

Expanded trapezoidal rule.

For a fixed function f(x) , which must be integrated between the fixed limits a and b , we can double the number of intervals in the extended trapezoidal rule without losing the advantages of the previous work. The roughest implementation of the trapezoidal rule is to average the function at its ends a and b . The first step of refinement is to add to this average the value of the function in the middle of the path. The second refinement step is to add values ​​to the 1/4 and 3/4 tags.

enter image description here

Elementary quadrature algorithms include the addition of sequential refinement steps. It is convenient to encapsulate this function in the Quadrature structure:

 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 Trapzd structure Trapzd obtained from this as follows:

 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 fx / 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 * (ba) * (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; } } };

Using:

The Trapzd can be used in several ways. The simplest and toughest is the integration of the function according to the extended trapezoidal rule, where you know in advance the number of steps that you want. if you want 2^M + 1 , you can accomplish this with a snippet:

 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 Ftor is a functor containing an integrable function.

Bonus:

It is much better, of course, to clarify the trapezoid rule until a degree of accuracy is determined. For this function:

 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"); }

Type definitions.

 typedef double Doub; // 64 - bit floating point typedef int Int; // 32 - bit signed integer
+3
source

Source: https://habr.com/ru/post/1480386/

More articles:


All Articles