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