Sharedtest.cOpen in CoCalc
Author: Harald Schilly
Views : 14
1
#include <stdlib.h>
2
#include <gsl/gsl_math.h>
3
#include <gsl/gsl_monte.h>
4
#include <gsl/gsl_monte_miser.h>
5
#include <gsl/gsl_monte_plain.h>
6
#include <gsl/gsl_monte_vegas.h>
7
8
/* Computation of the integral,
9
10
I = int (dx dy dz)/(2pi)^3 1/(1-cos(x)cos(y)cos(z))
11
12
over (-pi,-pi,-pi) to (+pi, +pi, +pi). The exact answer
13
is Gamma(1/4)^4/(4 pi^3). This example is taken from
14
C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
15
Volume 1", Section 1.1, p21, which cites the original
16
paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
17
1800 (1977) */
18
19
/* For simplicity we compute the integral over the region
20
(0,0,0) -> (pi,pi,pi) and multiply by 8 */
21
22
double exact = 1.3932039296856768591842462603255;
23
24
double g(double *k, size_t dim, void *params) {
25
(void)(dim); /* avoid unused parameter warnings */
26
(void)(params);
27
double A = 1.0 / (M_PI * M_PI * M_PI);
28
return A / (1.0 - cos(k[0]) * cos(k[1]) * cos(k[2]));
29
}
30
31
void display_results(char *title, double result, double error) {
32
printf("%s ==================\n", title);
33
printf("result = % .6f\n", result);
34
printf("sigma = % .6f\n", error);
35
printf("exact = % .6f\n", exact);
36
printf("error = % .6f = %.2g sigma\n", result - exact,
37
fabs(result - exact) / error);
38
}
39
40
int main(void) {
41
double res, err;
42
43
double xl[3] = {0, 0, 0};
44
double xu[3] = {M_PI, M_PI, M_PI};
45
46
const gsl_rng_type *T;
47
gsl_rng *r;
48
49
gsl_monte_function G = {&g, 3, 0};
50
51
size_t calls = 500000;
52
53
gsl_rng_env_setup();
54
55
T = gsl_rng_default;
56
r = gsl_rng_alloc(T);
57
58
{
59
gsl_monte_plain_state *s = gsl_monte_plain_alloc(3);
60
gsl_monte_plain_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
61
gsl_monte_plain_free(s);
62
63
display_results("plain", res, err);
64
}
65
66
{
67
gsl_monte_miser_state *s = gsl_monte_miser_alloc(3);
68
gsl_monte_miser_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
69
gsl_monte_miser_free(s);
70
71
display_results("miser", res, err);
72
}
73
74
{
75
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(3);
76
gsl_monte_vegas_integrate(&G, xl, xu, 3, 10000, r, s, &res, &err);
77
display_results("vegas warm-up", res, err);
78
79
printf("converging...\n");
80
81
do {
82
gsl_monte_vegas_integrate(&G, xl, xu, 3, calls / 5, r, s, &res, &err);
83
printf(
84
"result = % .6f sigma = % .6f "
85
"chisq/dof = %.1f\n",
86
res, err, gsl_monte_vegas_chisq(s));
87
} while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.5);
88
89
display_results("vegas final", res, err);
90
91
gsl_monte_vegas_free(s);
92
}
93
94
gsl_rng_free(r);
95
96
return 0;
97
}
98