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 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 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 17 18 PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL)); 19 PetscCall(PetscOptionsGetInt(NULL,NULL,"-moment_max",&momentummax,NULL)); 20 PetscCall(PetscOptionsGetReal(NULL,NULL,"-sigma",&sigma,NULL)); 21 PetscCall(PetscOptionsGetReal(NULL,NULL,"-mu",&mu,NULL)); 22 23 /* calulate zeros and roots of Hermite Gauss quadrature */ 24 PetscCall(PetscMalloc1(n,&zeros)); 25 zeros[0] = 0; 26 tick = n % 2; 27 for (s=0; s<n/2; s++) { 28 zeros[2*s+tick] = -gsl_sf_hermite_zero(n,s+1); 29 zeros[2*s+1+tick] = gsl_sf_hermite_zero(n,s+1); 30 } 31 32 PetscCall(PetscDTFactorial(n, &scale)); 33 scale = exp2(n-1)*scale*PetscSqrtReal(PETSC_PI)/(n*n); 34 PetscCall(PetscMalloc1(n+1,&weights)); 35 for (s=0; s<n; s++) { 36 h = gsl_sf_hermite(n-1, (double) zeros[s]); 37 weights[s] = scale/(h*h); 38 } 39 /* zeros and weights verfied up to n = 5 with http://mathworld.wolfram.com/Hermite-GaussQuadrature.html */ 40 41 for (moment=0; moment<momentummax; moment++) { 42 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++) { 52 g += weights[s]*PetscPowRealInt(sqrt(2)*sigma*zeros[s] + mu,moment); 53 } 54 g /= sqrt(PETSC_PI); 55 /* results confirmed with https://en.wikipedia.org/wiki/Normal_distribution#Moments sigma^p * (p-1)!!*/ 56 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Moment %" PetscInt_FMT " %g \n",moment,(double)g)); 57 58 } 59 PetscCall(PetscFree(zeros)); 60 PetscCall(PetscFree(weights)); 61 PetscCall(PetscFinalize()); 62 return 0; 63 } 64 65 /*TEST 66 67 build: 68 requires: gsl double 69 70 test: 71 72 TEST*/ 73