| #include <stdlib.h> |
| #include <gsl/gsl_math.h> |
| #include <gsl/gsl_monte.h> |
| #include <gsl/gsl_monte_plain.h> |
| #include <gsl/gsl_monte_miser.h> |
| #include <gsl/gsl_monte_vegas.h> |
| |
| /* Computation of the integral, |
| |
| I = int (dx dy dz)/(2pi)^3 1/(1-cos(x)cos(y)cos(z)) |
| |
| over (-pi,-pi,-pi) to (+pi, +pi, +pi). The exact answer |
| is 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 */ |
| |
| double exact = 1.3932039296856768591842462603255; |
| |
| double |
| g (double *k, size_t dim, 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])); |
| } |
| |
| void |
| display_results (char *title, double result, double error) |
| { |
| printf ("%s ==================\n", title); |
| printf ("result = % .6f\n", result); |
| printf ("sigma = % .6f\n", error); |
| printf ("exact = % .6f\n", exact); |
| printf ("error = % .6f = %.1g sigma\n", result - exact, |
| fabs (result - exact) / error); |
| } |
| |
| int |
| main (void) |
| { |
| 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 ("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 ("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 ("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, s->chisq); |
| } |
| while (fabs (s->chisq - 1.0) > 0.5); |
| |
| display_results ("vegas final", res, err); |
| |
| gsl_monte_vegas_free (s); |
| } |
| |
| gsl_rng_free (r); |
| |
| return 0; |
| } |