gsl problem

[Sergey K]

Hello!

I need to calculate the integral:

http://s017.radikal.ru/i430/1601/17/cd5ecc24a40d.jpg

I use gsl library. Unfortunately, the library generates an error. what it might be? maybe I'm not properly convey koofitsenty in integrand?

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#include <gsl_integration.h>
#include <gsl_math.h>
#include <gsl/gsl_sf_pow_int.h>



double q = 120;
double h = 0.1;
int N_0 = 300;
double T = 300;

double plot = 7310;
double e_sub = 2.52;
double m_at = 114 * 1830;
int gam = 3;
int agam = 3;

const double h_pl = 1.05457e-34;
const double k = 1.3806503e-23;
const double m_e = 9.10938188e-31;

double E_sr;



double my_integrand (double x, void *params_ptr) {

// The next line recovers alpha from the passed params pointer
double alpha = *(double *) params_ptr;
double lower_limit = *(double *) params_ptr;

return ((pow(x,alpha-1)*pow(2.71828,-x)));
}

double NepGamma(double alpha, double lower_limit ) {
gsl_integration_workspace *work_ptr =
gsl_integration_workspace_alloc (1000);

/* start integral from 1 (to infinity) */
double abs_error = 1.0e-8; /* to avoid round-off problems */
double rel_error = 1.0e-8; /* the result will usually be much better */
double result; /* the result from the integration */
double error; /* the estimated error from the integration */

double expected = -0.16442310483055015762; /* known answer */

gsl_function My_function;
void *params_ptr = &alpha;

My_function.function = &my_integrand;
My_function.params = params_ptr;

gsl_integration_qagiu (&My_function, lower_limit,
abs_error, rel_error, 1000, work_ptr, &result,
&error);

return result;
}

double functionS(double E_int, void *param_ptr) {

double N1 = *(double *) param_ptr;
double n = *(double *) param_ptr;

double znam1;
double chisl;

double znam2;
double chis2;

znam1 = 2*(N1 * pow(q,2)/(2*m_at)-E_int);
chisl = (44.5466 * pow((E_sr/3),1.5));

znam2 = -3*(N1 * pow(q,2)/(2*m_at)-E_int);
chis2 = 2*E_sr;

return ( (sqrt(znam1)/chisl * exp(znam2/chis2))*1.1283791670955126 * NepGamma(1.5, 1.5*( (e_sub*pow(n,0.6666666666)*0.2)/(E_sr) )) );

}

double integr2(int N1, int n){

double down_limit = 0;
double up_limit = ((N1+n)*E_sr-h*e_sub*pow(N1+n,0.666));

int max_intervals = 10000;

gsl_integration_workspace * w = gsl_integration_workspace_alloc(max_intervals);


double abserr = 0.0;
double rellerr = 1e-8;
/* the result will usually be much better */
double result; /* the result from the integration */
double error; /* the estimated error from the integration */

double expected = -0.16442310483055015762; /* known answer */

gsl_function My_function;

void *param_ptr = &N1;

gsl_function F = {&functionS, 0};
F.params = param_ptr;


gsl_integration_qags(&F,down_limit,up_limit,abserr, rellerr , max_intervals, w, &result, &error);
gsl_integration_workspace_free(w);

return result;
}

int main(int argc, char *argv[]) {

E_sr = ((q*q)/(2*m_at))*27.21;

printf("\n Value Energy med per atom = %lf", E_sr);
printf("\n");
printf("\n");

printf("\n");
printf("\n Test: %e", integr2(1, 1) );
printf("\n");


getch();

return 0;
}

[Moschops]

Unfortunately, the library generates an error. what it might be?

No idea. We could be here trying to guess the error all day. It would be faster if you just told us the error.

[Sergey K]

You're right. Error here.

http://s019.radikal.ru/i628/1601/be/244b75f01c3c.png