How to obtain a function from gsl library integration?

Click For Summary
SUMMARY

The discussion focuses on integrating functions using the GNU Scientific Library (GSL) with Monte Carlo methods, specifically addressing the challenge of incorporating dependencies on multiple variables. The user initially provided a simple integration example using GSL's Monte Carlo Vegas algorithm and sought guidance on how to evaluate a function of the form f(x,y) = x*y after integrating over x. A solution was presented, demonstrating how to define a function with parameters and integrate it while accounting for dependencies on additional variables.

PREREQUISITES
  • Familiarity with C++ programming language
  • Understanding of Monte Carlo integration techniques
  • Knowledge of the GNU Scientific Library (GSL) version 2.6 or later
  • Basic concepts of function parameters and struct usage in C++
NEXT STEPS
  • Explore GSL's Monte Carlo integration functions, specifically gsl_monte_vegas_integrate
  • Learn about struct usage in C++ for passing parameters to functions
  • Investigate advanced Monte Carlo techniques for multi-dimensional integration
  • Review GSL documentation for additional examples and best practices
USEFUL FOR

Researchers, physicists, and software developers working on numerical simulations or mathematical modeling that require integration of functions with multiple variables using the GSL library.

Fabio Kopp
Messages
16
Reaction score
0
I am having some problems to use the gsl library. I do the integration via monte carlo and I have obtained the right result, but a need to know how to use this result when I have a dependence on another variable. For example, f(x,y)=x*y, but the integration is only on x. Because I want to evaluate this function on y after the integration.
For the simple case, only one integration without the dependence in another variables, I obtained using this code:
Mod note: Edited to add [ code ] and [ /code ] tags.
C:
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <cmath>

using namespace std;

double g (double *k, size_t dim, void *params)
{
  double A;
   A=1*k[0];
    return A;}int main(void)
{
  double res, err;  double xl[10] = { 0.0, 0.0,0.0};
  double xu[10] = { 1.0, 1.0,1.0};

  const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G = { &g, 3, 0 };
  gsl_rng_env_setup ();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

{

gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);
gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
&res, &err);

do
{
size_t calls = 50000;

gsl_monte_vegas_integrate (&G, xl, xu, 3, calls, r, s,
&res, &err);
}
while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.1);

count<<"The result is "<< res <<"."<<endl;
gsl_monte_vegas_free (s);
}
  gsl_rng_free (r);

  return 0;
}

To compile the above code, I used:

c++ -Wall -I/usr/local/include -c teste00.cc
c++ -L/usr/local/lib teste00.o -lgsl -lgslcblas -lmI have a function that has the following form
A=x*k[0].
As I said before, I do not know how to implement this dependence in x to be evaluated after the integration on GSL.
Does anybody have any idea how to do this?
 
Last edited by a moderator:
Technology news on Phys.org
Fabio Kopp said:
I am having some problems to use the gsl library. I do the integration via monte carlo and I have obtained the right result, but a need to know how to use this result when I have a dependence on another variable. For example, f(x,y)=x*y, but the integration is only on x. Because I want to evaluate this function on y after the integration.
For the simple case, only one integration without the dependence in another variables, I obtained using this code:
Mod note: Edited to add [ code ] and [ /code ] tags.
C:
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>
#include <cmath>

using namespace std;

double g (double *k, size_t dim, void *params)
{
  double A;
   A=1*k[0];
    return A;}int main(void)
{
  double res, err;  double xl[10] = { 0.0, 0.0,0.0};
  double xu[10] = { 1.0, 1.0,1.0};

  const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G = { &g, 3, 0 };
  gsl_rng_env_setup ();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

{

gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);
gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
&res, &err);

do
{
size_t calls = 50000;

gsl_monte_vegas_integrate (&G, xl, xu, 3, calls, r, s,
&res, &err);
}
while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.1);

count<<"The result is "<< res <<"."<<endl;
gsl_monte_vegas_free (s);
}
  gsl_rng_free (r);

  return 0;
}

To compile the above code, I used:

c++ -Wall -I/usr/local/include -c teste00.cc
c++ -L/usr/local/lib teste00.o -lgsl -lgslcblas -lmI have a function that has the following form
A=x*k[0].
As I said before, I do not know how to implement this dependence in x to be evaluated after the integration on GSL.
Does anybody have any idea how to do this?
I have found the solution for this problem and the code is the following.
#include<iostream>
#include<gsl/gsl_math.h>
#include <gsl/gsl_monte_vegas.h>
#include <gsl/gsl_monte.h>

using namespace std;

struct my_f_params { double a; double b; double c; };

double my_f (double x[], size_t dim, void * p) {
struct my_f_params * fp = (struct my_f_params *)p;
if (dim != 2)
{ fprintf (stderr, "error: dim != 2");
abort ();
}
return fp->a * x[0] * x[0]
+ fp->b * x[0] * x[1]
+ fp->c * x[1] * x[1];
}int main(){

gsl_monte_function F;
for (double i=1.0; i<10;++i){
struct my_f_params params = { i, 2.0, 1.0 };
F.f = &my_f;
F.dim = 2;
F.params = &params;double res, err;
double xl[10] = { 0.0, 0.0,0.0};
double xu[10] = { 1.0, 1.0,1.0};

const gsl_rng_type *T;
gsl_rng *r;
gsl_rng_env_setup ();
T = gsl_rng_default;
r = gsl_rng_alloc (T);

{

gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);
gsl_monte_vegas_integrate (&F, xl, xu, 3, 10000, r, s,
&res, &err);

do
{
size_t calls = 50000;gsl_monte_vegas_integrate (&F, xl, xu, 3, calls, r, s,
&res, &err);
}
while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.1);

count<<" For i = "<< i << " the result is " << res << "."<<endl;gsl_monte_vegas_free (s);

}
gsl_rng_free (r);
}
return 0;
}
 

Similar threads

  • · Replies 11 ·
Replies
11
Views
2K
  • · Replies 2 ·
Replies
2
Views
2K
Replies
4
Views
11K
Replies
5
Views
2K
  • · Replies 4 ·
Replies
4
Views
2K
Replies
3
Views
2K