C 代码 应用蒙特卡罗方法估计函数的积分 在 3D 单元球的内部.rar

# include <math.h> # include <stdlib.h> # include <stdio.h> # include <time.h> # include "ball_monte_carlo.h" /******************************************************************************/ double ball01_monomial_integral ( int e[3] ) /******************************************************************************/ /* Purpose: BALL01_MONOMIAL_INTEGRAL returns monomial integrals in the unit ball. Discussion: The integration region is X^2 + Y^2 + Z^2 <= 1. The monomial is F(X,Y,Z) = X^E(1) * Y^E(2) * Z^E(3). Licensing: This code is distributed under the GNU LGPL license. Modified: 02 January 2014 Author: John Burkardt Reference: Philip Davis, Philip Rabinowitz, Methods of Numerical Integration, Second Edition, Academic Press, 1984, page 263. Parameters: Input, int E[3], the exponents of X, Y and Z in the monomial. Each exponent must be nonnegative. Output, double BALL01_MONOMIAL_INTEGRAL, the integral. */ { int i; double integral; const double r = 1.0; double s; if ( e[0] < 0 || e[1] < 0 || e[2] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "BALL01_MONOMIAL_INTEGRAL - Fatal error!\n" ); fprintf ( stderr, " All exponents must be nonnegative.\n" ); fprintf ( stderr, " E[0] = %d\n", e[0] ); fprintf ( stderr, " E[1] = %d\n", e[1] ); fprintf ( stderr, " E[2] = %d\n", e[2] ); exit ( 1 ); } if ( ( e[0] % 2 ) == 1 || ( e[1] % 2 ) == 1 || ( e[2] % 2 ) == 1 ) { integral = 0.0; } else { integral = 2.0; for ( i = 0; i < 3; i++ ) { integral = integral * tgamma ( 0.5 * ( double ) ( e[i] + 1 ) ); } integral = integral / tgamma ( 0.5 * ( double ) ( e[0] + e[1] + e[2] + 3 ) ); } /* The surface integral is now adjusted to give the volume integral. */ s = e[0] + e[1] + e[2] + 3; integral = integral * pow ( r, s ) / ( double ) ( s ); return integral; } /******************************************************************************/ double *ball01_sample ( int n, int *seed ) /******************************************************************************/ /* Purpose: BALL01_SAMPLE uniformly samples points from the unit ball. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 January 2014 Author: John Burkardt Reference: Russell Cheng, Random Variate Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, pages 168. Reuven Rubinstein, Monte Carlo Optimization, Simulation, and Sensitivity of Queueing Networks, Krieger, 1992, ISBN: 0894647644, LC: QA298.R79. Parameters: Input, int N, the number of points. Input/output, int *SEED, a seed for the random number generator. Output, double X[3*N], the points. */ { const double exponent = 1.0 / 3.0; int i; int j; double norm; double r; double *x; x = r8mat_normal_01_new ( 3, n, seed ); for ( j = 0; j < n; j++ ) { /* Compute the length of the vector. */ norm = sqrt ( pow ( x[0+j*3], 2 ) + pow ( x[1+j*3], 2 ) + pow ( x[2+j*3], 2 ) ); /* Normalize the vector. */ for ( i = 0; i < 3; i++ ) { x[i+j*3] = x[i+j*3] / norm; } /* Now compute a value to map the point ON the sphere INTO the sphere. */ r = r8_uniform_01 ( seed ); for ( i = 0; i < 3; i++ ) { x[i+j*3] = pow ( r, exponent ) * x[i+j*3]; } } return x; } /******************************************************************************/ double ball01_volume ( ) /******************************************************************************/ /* Purpose: BALL01_VOLUME returns the volume of the unit ball. Licensing: This code is distributed under the GNU LGPL license. Modified: 02 January 2014 Author: John Burkardt Parameters: Output, double BALL01_VOLUME, the volume of the unit ball. */ { const double r = 1.0; const double r8_pi = 3.141592653589793; double volume; volume = 4.0 * r8_pi * pow ( r, 3 ) / 3.0; return volume; } /******************************************************************************/ double *monomial_value ( int m, int n, int e[], double x[] ) /******************************************************************************/ /* Purpose: MONOMIAL_VALUE evaluates a monomial. Discussion: This routine evaluates a monomial of the form product ( 1 <= i <= m ) x(i)^e(i) where the exponents are nonnegative integers. Note that if the combination 0^0 is encountered, it should be treated as 1. Licensing: This code is distributed under the GNU LGPL license. Modified: 08 May 2014 Author: John Burkardt Parameters: Input, int M, the spatial dimension. Input, int N, the number of points at which the monomial is to be evaluated. Input, int E[M], the exponents. Input, double X[M*N], the point coordinates. Output, double MONOMIAL_VALUE[N], the value of the monomial. */ { int i; int j; double *v; v = ( double * ) malloc ( n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { v[j] = 1.0; } for ( i = 0; i < m; i++ ) { if ( 0 != e[i] ) { for ( j = 0; j < n; j++ ) { v[j] = v[j] * pow ( x[i+j*m], e[i] ); } } } return v; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1]. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the GNU LGPL license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { int i4_huge = 2147483647; int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ double *r8mat_normal_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_NORMAL_01_NEW returns a unit pseudonormal R8MAT. Licensing: This code is distributed under the GNU LGPL license. Modified: 03 October 2005 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Si