Actual source code: ex26.c

      1: static char help[] = "Calculates moments for Gaussian functions.\n\n";

      3: #include <petscviewer.h>
      4: #include <petscdt.h>
      5: #include <petscvec.h>

      7: #include <gsl/gsl_sf_hermite.h>
      8: #include <gsl/gsl_randist.h>

     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;

     16:   PetscFunctionBeginUser;
     17:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

     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));

     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:   }

     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 */

     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

     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: }

     62: /*TEST

     64:   build:
     65:     requires: gsl double

     67:   test:

     69: TEST*/