CoCalc Shared Filesgsl-example.ipynbOpen in CoCalc with one click!
Author: Harald Schilly
Views : 62
Description: gsl in xeus c++ 17

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") #pragma cling load("gslcblas")
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=(pi,pi,pi)(+pi,+pi,+pi)dxdydz(2π)311cos(x)cos(y)cos(z)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)}

The exact answer is

Γ(1/4)44pi3\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 [ ]: