C 代码 应用蒙特卡罗方法估计函数的积分 在2D普通圆盘的内部.rar

# include <math.h> # include <stdio.h> # include <stdlib.h> # include <time.h> # include "disk_monte_carlo.h" /******************************************************************************/ double disk_area ( double center[2], double r ) /******************************************************************************/ /* Purpose: DISK_AREA returns the area of a general disk in 2D. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 July 2018 Author: John Burkardt Parameters: Input, double CENTER[2], coordinates of the center. This information is not actually needed to compute the area. Input, double R, the radius of the disk. Output, double DISK_AREA, the area of the disk. */ { double area; const double r8_pi = 3.141592653589793; area = r8_pi * r * r; return area; } /******************************************************************************/ double *disk_sample ( double center[2], double r, int n, int *seed ) /******************************************************************************/ /* Purpose: DISK_SAMPLE uniformly samples a general disk in 2D. Licensing: This code is distributed under the GNU LGPL license. Modified: 05 July 2018 Author: John Burkardt Parameters: Input, double CENTER[2], coordinates of the center. Input, double R, the radius of the disk. Input, int N, the number of points. Input/output, int *SEED, a seed for the random number generator. Output, double X[2*N], the points. */ { int i; int j; double norm; double r2; double *v; double *x; x = ( double * ) malloc ( 2 * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { v = r8vec_normal_01_new ( 2, seed ); /* Compute the length of the vector. */ norm = sqrt ( pow ( v[0], 2 ) + pow ( v[1], 2 ) ); /* Normalize the vector. */ for ( i = 0; i < 2; i++ ) { v[i] = v[i] / norm; } /* Now compute a value to map the point ON the circle INTO the circle. */ r2 = r8_uniform_01 ( seed ); for ( i = 0; i < 2; i++ ) { x[i+j*2] = center[i] + r * sqrt ( r2 ) * v[i]; } free ( v ); } return x; } /******************************************************************************/ double disk01_monomial_integral ( int e[2] ) /******************************************************************************/ /* Purpose: DISK01_MONOMIAL_INTEGRAL returns monomial integrals in the unit disk in 2D. Discussion: The integration region is X^2 + Y^2 <= 1. The monomial is F(X,Y) = X^E(1) * Y^E(2). Licensing: This code is distributed under the GNU LGPL license. Modified: 03 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[2], the exponents of X and Y in the monomial. Each exponent must be nonnegative. Output, double DISK01_MONOMIAL_INTEGRAL, the integral. */ { double arg; int i; double integral; const double r = 1.0; double s; if ( e[0] < 0 || e[1] < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "DISK01_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] ); exit ( 1 ); } if ( ( e[0] % 2 ) == 1 || ( e[1] % 2 ) == 1 ) { integral = 0.0; } else { integral = 2.0; for ( i = 0; i < 2; i++ ) { arg = 0.5 * ( double ) ( e[i] + 1 ); integral = integral * tgamma ( arg ); } arg = 0.5 * ( double ) ( e[0] + e[1] + 2 ); integral = integral / tgamma ( arg ); } /* Adjust the surface integral to get the volume integral. */ s = e[0] + e[1] + 2; integral = integral * pow ( r, s ) / ( double ) ( s ); return integral; } /******************************************************************************/ 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 *r8vec_normal_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_NORMAL_01_NEW returns a unit pseudonormal R8VEC. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. Licensing: This code is distributed under the GNU LGPL license. Modified: 06 August 2013 Author: John Burkardt Parameters: Input, int N, the number of values desired. Input/output, int *SEED, a seed for the random number generator. Output, dou