Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How to obtain a function from gsl library integration?

  1. Nov 6, 2015 #1
    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.
    Code (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);

    cout<<"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 -lm


    I 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: Nov 6, 2015
  2. jcsd
  3. Nov 8, 2015 #2

    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);

    cout<<" For i = "<< i << " the result is " << res << "."<<endl;


    gsl_monte_vegas_free (s);

    }
    gsl_rng_free (r);
    }
    return 0;
    }
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: How to obtain a function from gsl library integration?
Loading...