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