[Help-gsl] GSL Integration on Data Set

[Amol Holkundkar]

Greetings,


I have used gsl_integration_qag to integrate the data provided by 
vector<double>.

In practical scenario this actually happens, so just attempted to translate the

vector<double> as gsl_fucntion and then do the integration. Hope it would be

useful for some users.


/*******************************************************************************************************/


//Author: Amol Holkundkar, Department of Physics, BITS Pilani, Pilani, India.
//This code demonstrate the GSL integration routine on some array. In applied
//scenario we need to integrate some set of data points. Here,
//the function DoIntegration(from,to,dx,fx,&result,&error) will directly
//return the result for a function which is passed  as vector<double>.
//$ g++ -std=c++11 -lgslcblas -lgsl this_file.cpp

#include <iostream>
#include <fstream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <vector>
#include <gsl/gsl_integration.h>
using namespace std;
typedef vector<double> VD;

double dx = 1e-5;
double a = 0.0;            /* lower limit a */
double b = M_PI;    /* upper limit b */

double Function(double x)
{
 return x*x*sin(x); //Analytical Integration Can be done, result = pi*pi - 4
}
double exact_result = M_PI*M_PI - 4.0 ;

void DoIntegration(double from,double to,double dx0,VD fx,double &result,double 
&error)
{

  double abs_err = dx0*100; // to avoid round off errors
  double rel_err = dx0*100; // to avoid roung off errors

  VD alpha = fx;
  int N = fx.size();

  alpha.push_back(dx0);

  gsl_integration_workspace *w = gsl_integration_workspace_alloc (fx.size());
  gsl_function F;
  void *param = &alpha;
  auto integrand = [] (double x, void *p){VD al = *(VD *) p;return 
al[int(x/al[al.size()-1])];};

  F.function = integrand;
  F.params = param;

  gsl_integration_qag(&F,from,to,abs_err,rel_err,N,GSL_INTEG_GAUSS41,w, 
&result,&error);

  gsl_integration_workspace_free (w);
}


int main ()
{
  VD fx;

  double res,err;

  for(double x = a;x <= b;x += dx)
  {
   fx.push_back(  Function(x) ) ;
  }

  DoIntegration(a,b,dx,fx,res,err);

  cout << setprecision(15) << scientific << uppercase;
  cout << endl << "Result           = " <<  res;
  cout << endl << "Exact Result     = " << exact_result;
  cout << endl << "Estimated Error  = " <<  err;
  cout << endl << "Actual Error     = " << abs(res-exact_result);
  cout << endl << endl;

  return 0;
}


/*******************************************************************************************************/

/*******************************************************************************************************/



Cheers

Amol Holkundkar

gsl_integration_vector.cpp
Description: gsl_integration_vector.cpp