1 static char help[] = "Calculates moments for Gaussian functions.\n\n";
2
3 #include <petscviewer.h>
4 #include <petscdt.h>
5 #include <petscvec.h>
6
7 #include <gsl/gsl_sf_hermite.h>
8 #include <gsl/gsl_randist.h>
9
main(int argc,char ** argv)10 int main(int argc, char **argv)
11 {
12 int s, n = 15;
13 PetscInt tick, moment = 0, momentummax = 7;
14 PetscReal *zeros, *weights, scale, h, sigma = 1 / sqrt(2), g = 0, mu = 0;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
18
19 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
20 PetscCall(PetscOptionsGetInt(NULL, NULL, "-moment_max", &momentummax, NULL));
21 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma", &sigma, NULL));
22 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &mu, NULL));
23
24 /* calculate zeros and roots of Hermite Gauss quadrature */
25 PetscCall(PetscMalloc1(n, &zeros));
26 zeros[0] = 0;
27 tick = n % 2;
28 for (s = 0; s < n / 2; s++) {
29 zeros[2 * s + tick] = -gsl_sf_hermite_zero(n, s + 1);
30 zeros[2 * s + 1 + tick] = gsl_sf_hermite_zero(n, s + 1);
31 }
32
33 PetscCall(PetscDTFactorial(n, &scale));
34 scale = exp2(n - 1) * scale * PetscSqrtReal(PETSC_PI) / (n * n);
35 PetscCall(PetscMalloc1(n + 1, &weights));
36 for (s = 0; s < n; s++) {
37 h = gsl_sf_hermite(n - 1, (double)zeros[s]);
38 weights[s] = scale / (h * h);
39 }
40 /* zeros and weights verified up to n = 5 with http://mathworld.wolfram.com/Hermite-GaussQuadrature.html */
41
42 for (moment = 0; moment < momentummax; moment++) {
43 /* http://www.wouterdenhaan.com/numerical/integrationslides.pdf */
44 /* https://en.wikipedia.org/wiki/Gauss-Hermite_quadrature */
45 /*
46 int_{-infinity}^{infinity} \frac{1}{sigma sqrt(2pi)} exp(- \frac{(x - mu)^2}{2 sigma^2) h(x) dx
47
48 then approx equals 1/pi sum_i w_i h( sqrt(2) sigma x_i + mu)
49 */
50 g = 0;
51 for (s = 0; s < n; s++) g += weights[s] * PetscPowRealInt(sqrt(2) * sigma * zeros[s] + mu, moment);
52 g /= sqrt(PETSC_PI);
53 /* results confirmed with https://en.wikipedia.org/wiki/Normal_distribution#Moments sigma^p * (p-1)!!*/
54 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Moment %" PetscInt_FMT " %g \n", moment, (double)g));
55 }
56 PetscCall(PetscFree(zeros));
57 PetscCall(PetscFree(weights));
58 PetscCall(PetscFinalize());
59 return 0;
60 }
61
62 /*TEST
63
64 build:
65 requires: gsl double
66
67 test:
68
69 TEST*/
70