CoCalc Public Filesgsl-example.ipynb
Author: Harald Schilly
Views : 147
Description: gsl in xeus c++ 17
Compute Environment: Ubuntu 18.04 (Deprecated)

# XEUS C++17 running GSL

Tipp: Use the doube-arrow ("fast forward") to restart and run all cells after any changes!

In [1]:
#pragma cling add_include_path("/usr/include/")

In [2]:
#pragma cling load("gsl")

In [3]:
#include <stdlib.h>
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>


Computation of the integral,

$I = \int_{(-pi,-pi,-pi)}^{(+pi, +pi, +pi)} \frac{\mathrm{d}x \mathrm{d}y \mathrm{d}z}{(2 \pi)^3} \frac{1}{1-cos(x) \cos(y) \cos(z)}$

$\frac{\Gamma(1/4)^4}{4 pi^3}$

This example is taken from C.Itzykson, J.M.Drouffe, "Statistical Field Theory - Volume 1", Section 1.1, p21, which cites the original paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74 1800 (1977)

For simplicity we compute the integral over the region (0,0,0) -> (pi,pi,pi) and multiply by 8

In [4]:
double exact = 1.3932039296856768591842462603255;

In [5]:
double g(double *k, size_t dim, void *params) {
(void)(dim); /* avoid unused parameter warnings */
(void)(params);
double A = 1.0 / (M_PI * M_PI * M_PI);
return A / (1.0 - cos(k[0]) * cos(k[1]) * cos(k[2]));
}

In [6]:
void display_results(char *title, double result, double error) {
std::cout << "================== " << title << " ==================\n";
std::cout << "result = " << result << "\n";
std::cout << "sigma  = " << error << "\n";
std::cout << "exact  = " << exact << "\n";
std::cout << "error  = " << result - exact << " = " << fabs(result - exact) / error << " sigma\n";
}

In [7]:
double res, err;

double xl[3] = {0, 0, 0};
double xu[3] = {M_PI, M_PI, M_PI};

const gsl_rng_type *T;
gsl_rng *r;

gsl_monte_function G = {&g, 3, 0};

size_t calls = 500000;

gsl_rng_env_setup();

T = gsl_rng_default;
r = gsl_rng_alloc(T);

{
gsl_monte_plain_state *s = gsl_monte_plain_alloc(3);
gsl_monte_plain_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
gsl_monte_plain_free(s);

display_results((char*)"plain", res, err);
}

{
gsl_monte_miser_state *s = gsl_monte_miser_alloc(3);
gsl_monte_miser_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
gsl_monte_miser_free(s);

display_results((char*)"miser", res, err);
}

{
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(3);
gsl_monte_vegas_integrate(&G, xl, xu, 3, 10000, r, s, &res, &err);
display_results((char*)"vegas warm-up", res, err);

printf("converging...\n");

do {
gsl_monte_vegas_integrate(&G, xl, xu, 3, calls / 5, r, s, &res, &err);
printf(
"result = % .6f sigma = % .6f "
"chisq/dof = %.1f\n",
res, err, gsl_monte_vegas_chisq(s));
} while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.5);

display_results((char*)"vegas final", res, err);

gsl_monte_vegas_free(s);
}

gsl_rng_free(r);


================== plain ================== result = 1.41221 sigma = 0.0134359 exact = 1.3932 error = 0.0190048 = 1.41448 sigma ================== miser ================== result = 1.39132 sigma = 0.00346056 exact = 1.3932 error = -0.00188235 = 0.543944 sigma ================== vegas warm-up ================== result = 1.39267 sigma = 0.00341041 exact = 1.3932 error = -0.000531339 = 0.155799 sigma ================== vegas final ================== result = 1.39328 sigma = 0.00036248 exact = 1.3932 error = 7.74556e-05 = 0.213682 sigma
In [ ]:


In [ ]: