Integration issues

ROOT
1

Hi there,

I’d like to integrate a rectangular function (see below).
I’ve tested several configurations (integration type, function range (xmin,xmax), function position).
The results are given in results.png (non zero values are highlighted in red). The code is also attached.
Only the GaussLegendreIntegrator works whatever the range and position used (GaussLegendreIntegrator allows to specify the number of steps used in the integration process).
Unfortunately, results are approximate.
Correct results (integral=1) can only be obtained with the GSLIntegrator but this integrator seems very sensitive to the function range and position (see results.png).
Is there a way to cope with that (e.g. adjust the number of steps in GSLIntegrator)?

Furthermore, I didn’t manage to use the plug-in manager with usual “fit-devoted” functions: Double_t F_myFct ( Double_t *x, Double_t *par )
In the examples given in
root.cern.ch/drupal/content/function-integration
only “basic” functions are used.
Switching to “fit-devoted” functions makes ACLIC complain.
How to deal with that too?

TIA,
Z


results.png1999×475 153 KB


zzztest.c (8.43 KB)

2

Hello,
Thank you for your nice studies of all the different integration algorithm.

In your case the function is having a discontinuity in the derivatives and that can cause problems in the numerical integration.
The best is to use the GSL algorithm with the method QAGP where the singular points are passed to the integrand
See project-mathlibs.web.cern.ch/pro … ae6cab31bd

In this case you will always get the correct results, otherwise depending on how the intervals are constructed, you will include or not your rectangular region. You cannot change in this algorithm the number of points, just changing the rule.

Concerning your second question you would like to use directly a function defined as f(double *x, double *p) ?
You can use a TF1 or

You can do also:
double par[] = { …}; // parameter vector
ROOT::Math::Functor func(std::bind2nd(std::ptr_fun(f),par),1 );
ROOT::Math::OneDimFunctionAdapter<> wf(func);

and use wf in the integrator

Best Regards

Lorenzo

3

Thanks Lorenzo!

For the first point, everything’s now fine :smiley:
project-mathlibs.web.cern.ch/pro … e4c2dbbc05

For the second, I don’t even manage to compile the following piece of code :cry:

[code]#include <Math/Functor.h>
#include <Math/OneDimFunctionAdapter.h>
#include
using namespace std ;

Double_t F_rect ( Double_t *x, Double_t *p ) { return (p[0]*x[0]) ; }

void zzzztest ()
{
Double_t par[] = {2.} ;
ROOT::Math::Functor func ( std::bind2nd ( std::ptr_fun(F_rect) , par ), 1 ) ;
//ROOT::Math::OneDimFunctionAdapter<> wf_ (func) ;
//ROOT::Math::Integrator ig (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR, 1e-9, 1e-6, 1000) ;
//printf ("\n%e",ig.Integral(xmin,xmax)) ;
return ;
}
[/code]
Cheers,
Z

4

Sorry, I typed too fast, this is a working code

#include <Math/Functor.h> 
#include <Math/OneDimFunctionAdapter.h> 
#include <functional> 
#include "Math/Integrator.h"
using namespace std ; 

Double_t F_rect ( const Double_t *x, const Double_t *p ) { return (p[0]*x[0]) ; } 

void testBind (double xmin=0, double xmax = 1) 
{ 
  Double_t par[] = {2.} ; 
  ROOT::Math::Functor func ( std::bind2nd ( std::ptr_fun(F_rect) , par ), 1 ) ; 
  ROOT::Math::OneDimMultiFunctionAdapter<> wf_ (func) ; 
  ROOT::Math::Integrator ig (wf_, ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR, 1e-9, 1e-6, 1000) ; 
  printf ("\n%e\n",ig.Integral(xmin,xmax)) ; 
  return ; 
}

Your function must be of type f(const double *, const double *)
Cheers
Lorenzo

5

Thanks Lorenzo :smiley:
Attached’s a new program that needs to be compiled and that takes into account every case (with and without singularities specified).
The results are given in the .pdf file.
Z

PS: Many function numbers of GaussIntegrator and GaussLegendreIntegrator are said to be “not implemented”.
Should one read not YET implemented?
zzztest.c (9.76 KB)
zzztest.pdf (23.7 KB)

6

Hello,

Thank you for the results and the test program. If you don;t mind I’will probably create a full test program which I will add in our test suite.

Just few comments:
[ul]
ROOT::Math::Integrator is an interface for all types of integrators, including GSLIntegrator, so the results are the same if the same type of algorithm is used (ADAPTIVE, NONADAPTIVE, etc…)

I noticed just that the difference is that the default of GSLIntegrator is ADAPTIVESINGULAR while for ROOT::Math::Integrator is just ADAPTIVE. I will need to change it to make it consistent. This explains some of the difference you observe between using the two classes

Changing the rule makes sense only when you use ADAPTIVE

You can integrate with singularity only when using ADAPTIVESINGULAR as algorithm. This explains why you obtained a zero results when you run with ROOT::Math::Integrator
[/ul]

Some of the functionality is not provided in the GaussIntegrator or LegendreIntegrator (like integrating open intervals, singularity, etc…) so these functions are and will not be implemented, although they are available in the interface class.
I should verify that when these functions are called for those integrators a status code different than zero is set.

Cheers

Lorenzo

7

[quote] If you don’t mind I will probably create a full test program which I will add in our test suite.[/quote]No prob.

Ok for the comments and thanks for your help.
Cheers,
Z