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 109371c9d4SSatish Balay int main(int argc, char **argv) { 11e0f85f4dSBarry Smith int s, n = 15; 12e0f85f4dSBarry Smith PetscInt tick, moment = 0, momentummax = 7; 13e0f85f4dSBarry Smith PetscReal *zeros, *weights, scale, h, sigma = 1 / sqrt(2), g = 0, mu = 0; 14e0f85f4dSBarry Smith 15327415f7SBarry Smith PetscFunctionBeginUser; 169566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 17e0f85f4dSBarry Smith 189566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 199566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-moment_max", &momentummax, NULL)); 209566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-sigma", &sigma, NULL)); 219566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &mu, NULL)); 22e0f85f4dSBarry Smith 23e0f85f4dSBarry Smith /* calulate zeros and roots of Hermite Gauss quadrature */ 249566063dSJacob Faibussowitsch PetscCall(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 329566063dSJacob Faibussowitsch PetscCall(PetscDTFactorial(n, &scale)); 33e0f85f4dSBarry Smith scale = exp2(n - 1) * scale * PetscSqrtReal(PETSC_PI) / (n * n); 349566063dSJacob Faibussowitsch PetscCall(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 /* http://www.wouterdenhaan.com/numerical/integrationslides.pdf */ 43e0f85f4dSBarry Smith /* https://en.wikipedia.org/wiki/Gauss-Hermite_quadrature */ 44e0f85f4dSBarry Smith /* 45e0f85f4dSBarry Smith int_{-infinity}^{infinity} \frac{1}{sigma sqrt(2pi)} exp(- \frac{(x - mu)^2}{2 sigma^2) h(x) dx 46e0f85f4dSBarry Smith 47e0f85f4dSBarry Smith then approx equals 1/pi sum_i w_i h( sqrt(2) sigma x_i + mu) 48e0f85f4dSBarry Smith */ 49e0f85f4dSBarry Smith g = 0; 50*ad540459SPierre Jolivet for (s = 0; s < n; s++) g += weights[s] * PetscPowRealInt(sqrt(2) * sigma * zeros[s] + mu, moment); 51e0f85f4dSBarry Smith g /= sqrt(PETSC_PI); 52e0f85f4dSBarry Smith /* results confirmed with https://en.wikipedia.org/wiki/Normal_distribution#Moments sigma^p * (p-1)!!*/ 5363a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Moment %" PetscInt_FMT " %g \n", moment, (double)g)); 54e0f85f4dSBarry Smith } 559566063dSJacob Faibussowitsch PetscCall(PetscFree(zeros)); 569566063dSJacob Faibussowitsch PetscCall(PetscFree(weights)); 579566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 58b122ec5aSJacob Faibussowitsch return 0; 59e0f85f4dSBarry Smith } 60e0f85f4dSBarry Smith 61e0f85f4dSBarry Smith /*TEST 62e0f85f4dSBarry Smith 63e0f85f4dSBarry Smith build: 64e0f85f4dSBarry Smith requires: gsl double 65e0f85f4dSBarry Smith 66e0f85f4dSBarry Smith test: 67e0f85f4dSBarry Smith 68e0f85f4dSBarry Smith TEST*/ 69