Integral with an infinite upper limit

Hi

Somebody can help me?

I want to define a function in this way

Sqrt[3/(2rh) (1+(R/rh)^2)^2 NIntegrate[(M[y,Ms,rs]y^(2\[Beta]-2) FF[y,R,\[Beta]])/(1+(y/rh)^2)^(5/2),{y,R,Infinity}]]

where rh = 241 and Ms, rs and Beta are parameters to fit.

But I have problems with the integration

I define inside of my root file.c

Double_t f1(float x, float R,  Double_t *par){

    Double_t value3 = ((par[0]/(1.-(TMath::Pi())/4))*((x/par[1]) - TMath::ATan(x/par[1])))*TMath::Power(x,2.0*par[2]-2.0)*(0.5*TMath::Power(R,1.0 -2.0*par[2])*(TMath::Sqrt(TMath::Pi())*TMath::Gamma(-0.5+par[2])/TMath::Gamma(par[2]) - TMath::BetaIncomplete(TMath::Power(R/x,2), -0.5 + par[2], 0.5) + par[2]*TMath::BetaIncomplete(TMath::Power(R/x,2), 0.5 + par[2], 0.5)  -TMath::Sqrt(TMath::Pi())*TMath::Gamma(0.5+par[2])/TMath::Gamma(par[2])))/TMath::Power(1.0 + (TMath::Power(x/196.,2)), 2.5);
    return value3;
}

I really need help to define now $\int_R^\infty{f1(x,R,par)dx}$ inside the code.

Please help me

Thanks

Hi,

What are the problems?

  • Do you ask us to verify that you did the integration properly?
  • Do you have a compilation error?
  • Do you get wrong results?

Please be more specific, we cannot probe for all possible errors for all of our users :slight_smile:

Axel.

Ok. Sorry … I’m not good in this.

In my case these par (parameters p[0],p[1],p[2]) will be determinated by gMinuit but gsl_integration_qagiu need only a function f(x) to obtain a result. So, when I do this

Double_t integrand(float R, Double_t *par){
  
    gsl_integration_workspace *work_ptr =
    gsl_integration_workspace_alloc (1000);
    double result;		/* the result from the integration */
    double error;			/* the estimated error from the integration */
    
    gsl_function F;
    
    F.function = &f1;
    

    Double_t value6= gsl_integration_qagiu(&F,R,1.0e-8,1.0e-8,1000,work_ptr, &result,&error);
    return result;
}

has the problem of f(x,R,par). Root gives me this message

assigning to ‘double (*)(double, void )’ from incompatible type 'Double_t ()(float, Double_t *)’:
type mismatch at 1st parameter (‘double’ vs ‘float’)
F.function = &f1;

I need to define this function before doing the fit.

I hope this help to show what my problem is.

Thanks a lof for answer.

Change Double_t f1(float x, float R, Double_t *par){ to

Double_t f1(float x, float R,  void *vpar){
   Double_t* par = (Double_t*)vpar;

Does that help?

Cheers, Axel.

Hi, I try it, but it doesn’t work:

Double_t f1(float x, float R,  void *vpar){
    Double_t*par = (Double_t*)vpar;
    Double_t value3 = ((par[0]/(1.-(TMath::Pi())/4))*((x/par[1]) - TMath::ATan(x/par[1])))*TMath::Power(x,2.0*par[2]-2.0)*(0.5*TMath::Power(R,1.0 -2.0*par[2])*(TMath::Sqrt(TMath::Pi())*TMath::Gamma(-0.5+par[2])/TMath::Gamma(par[2]) - TMath::BetaIncomplete(TMath::Power(R/x,2), -0.5 + par[2], 0.5) + par[2]*TMath::BetaIncomplete(TMath::Power(R/x,2), 0.5 + par[2], 0.5)  -TMath::Sqrt(TMath::Pi())*TMath::Gamma(0.5+par[2])/TMath::Gamma(par[2])))/TMath::Power(1.0 + (TMath::Power(x/196.,2)), 2.5);
    return value3;
}

I get this

error: assigning to 'double (*)(double, void *)' from incompatible type 'Double_t (*)(float, float, void *)':
      different number of parameters (2 vs 3)
    F.function = &f1;

I think float R is a problem too, because

gsl_integration_qagiu(&F,R,1.0e-8,1.0e-8,1000,work_ptr, &result,&error); 

needs only a function f(x).

Cheers

Let’s see wether @moneta knows how that gsl_function is supposed to work…

Thank you so much.

I will do a change of variable y/R to avoid the problem of R inside of the ingration.

But I still want to define a function with gsl_integration_qagiu, that will be used after by gMinuit to find the parameters and minimize chisquare.

I really appreciate your help.

Hi,

If you want to use directly gel you need to use the right type. Look at the User guide of GSL as gal_function is defined, and you need to define your function exactly in the same way.

Otherwise if you use the Integrator class in ROOT, which is based on GSL you don;'t need all this and it is much simpler to use.

See https://root.cern.ch/doc/v608/classROOT_1_1Math_1_1IntegratorOneDim.html
or https://root.cern.ch/doc/v608/classROOT_1_1Math_1_1GSLIntegrator.html.

Examples are provided in https://root.cern.ch/doc/master/mathmoreIntegration_8C.html

Lorenzo

Thanks Lorenzo @moneta for your advice, but I feel lost

I need to use ROOT::Math::GSLIntegrator::IntegralUp

so I try to do this

Double_t Integrand(Double_t *par)
{

return ROOT::Math::GSLIntegrator::IntegralUp(&f1,1) ;
}

where

Double_t f1(float x, Double_t *par){

Double_t value3 =(x/par[1] - TMath::ATan(x/par[1]))*(TMath::Sqrt(TMath::Pi())*TMath::Gamma(-0.5+par[2])/TMath::Gamma(par[2]) - TMath::BetaIncomplete(TMath::Power(1./x,2), -0.5 + par[2], 0.5) + par[2]*TMath::BetaIncomplete(TMath::Power(1/x,2), 0.5 + par[2], 0.5)  -TMath::Sqrt(TMath::Pi())*TMath::Gamma(0.5+par[2])/TMath::Gamma(par[2]))/TMath::Power(1.0 + (TMath::Power(x/par[3],2)), 2.5);

return value3;
}

I have no idea how to use it.

I need to define this function Integrand because the Final function needs it to fit the parameters.

Lizbeth

@moneta Do you need the macro ??

Cheers

Yes, if you have a macro, it would be easier for me.

To use the IntegratorOneDim, you need to define the integrand function with this signature;

double f( double x) {

// here you implement your integrand function
}

Lorenzo

@moneta

Hi, thanks for the help

this is the macrofit.C (7.0 KB)
car.txt (2.1 KB)

@moneta

I will appreciate any help.

Thanks.

@moneta

Please help me with the integrand function.

Hi,

Thank you for the code. But I really do not understand completely what you would like to to.
Your function Sigma calls Integrand which seems to try to compute an integral between 1 and + infinity of the function f1(x, parameters). But you do not pass the parameters from Sigma to Integrand.

Can you maybe explains in English or Math formula (even better) what you would like to do ?

One more comment, I would not mix floating precision with double. Just use double everywhere instead of float.

Lorenzo

@moneta

Hi Lorenzo, thanks for your answer

I want to implement these equations
\begin{equation}
\sigma_{los}^2 ® = \frac{3}{2r_{half}}\left( 1 + \left(\frac{R}{r_{half}}\right)^2\right)^2\int_R^\infty \frac{M(y)y^{2\beta-2}FF(\beta,R,y)}{\left(1+ \left(\frac{y}{r_{half}}\right)^2\right)^{5/2}} dy
\end{equation}
where $ \beta, r_{half}, M_s, r_s$ are my fit parameters and

\begin{equation}
\text{FF}(R,y,\beta )=\frac{1}{2} R^{1-2 \beta } \left(\frac{\sqrt{\pi } \Gamma \left(\beta -\frac{1}{2}\right)}{\Gamma (\beta )}-\frac{\sqrt{\pi } \Gamma \left(\beta +\frac{1}{2}\right)}{\Gamma (\beta )}+\beta B\left[\frac{R^2}{y^2},\beta +\frac{1}{2},\frac{1}{2}\right]-B\left[\frac{R^2}{y^2},\beta -\frac{1}{2},\frac{1}{2}\right]\right),
\end{equation}
\begin{equation}
M®=\frac{M_s \left(\frac{r}{r_s}-\tan ^{-1}\left(\frac{r}{r_s}\right)\right)}{1-\frac{\pi }{4}},
\end{equation}

where $\sqrt{\sigma_{los}}$ is the function to fit the observational data.

I will appreciate your help and time.

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.