#! /usr/bin/env python # def triangle01_monte_carlo ( xy_num, triangle01_integrand, seed ): #*****************************************************************************80 # ## TRIANGLE01_MONTE_CARLO applies the Monte Carlo rule to integrate a function. # # Discussion: # # The function f(x,y) is to be integrated over the unit triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 04 November 2016 # # Author: # # John Burkardt # # Parameters: # # Input, integer XY_NUM, the number of sample points. # # Input, external TRIANGLE01_INTEGRAND, the integrand routine. # # Input/output, integer SEED, a seed for the random # number generator. # # Output, real RESULT(F_NUM), the approximate integrals. # import numpy as np from triangle01_sample import triangle01_sample area = 0.5 xy, seed = triangle01_sample ( xy_num, seed ) fxy = triangle01_integrand ( xy ) result = area * np.sum ( fxy ) / float ( xy_num ) return result, seed def triangle01_monte_carlo_test ( ): #*****************************************************************************80 # ## TRIANGLE01_MONTE_CARLO_TEST estimates integrals over the unit triangle. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 04 November 2016 # # Author: # # John Burkardt # import numpy as np from triangle01_monomial_integral import triangle01_monomial_integral global e m = 2 e_test = np.array ( [ \ [ 0, 0 ], \ [ 1, 0 ], \ [ 0, 1 ], \ [ 2, 0 ], \ [ 1, 1 ], \ [ 0, 2 ], \ [ 3, 0 ] ] ) print ( '' ) print ( 'TRIANGLE01_MONTE_CARLO_TEST01' ) print ( ' TRIANGLE01_MONTE_CARLO estimates an integral over' ) print ( ' the unit triangle using the Monte Carlo method.' ) seed = 123456789 print ( '' ) print ( ' N 1 X Y ' ), print ( ' X^2 XY Y^2 X^3' ) print ( '' ) n = 1 e = np.zeros ( m, dtype = np.int32 ) while ( n <= 65536 ): print ( ' %8d' % ( n ) ), for j in range ( 0, 7 ): e[0:m] = e_test[j,0:m] result, seed = triangle01_monte_carlo ( n, triangle01_integrand, seed ) print ( ' %14.6g' % ( result ) ), print ( '' ) n = 2 * n print ( '' ) print ( ' Exact' ), for j in range ( 0, 7 ): e[0:m] = e_test[j,0:m] result = triangle01_monomial_integral ( e[0], e[1] ) print ( ' %14.6g' % ( result ) ), print ( '' ) # # Terminate. # print ( '' ) print ( 'TRIANGLE01_MONTE_CARLO_TEST:' ) print ( ' Normal end of execution.' ) return def triangle01_integrand ( xy ): #*****************************************************************************80 # ## TRIANGLE01_INTEGRAND is a sample integrand function. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 04 November 2016 # # Author: # # John Burkardt # # Parameters: # # Input, real XY(XY_NUM), the number of evaluation points. # # Output, real FXY(XY_NUM), the function values. # from monomial_value import monomial_value m = 2 n = xy.shape[1] fxy = monomial_value ( m, n, e, xy ) return fxy if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) triangle01_monte_carlo_test ( ) timestamp ( )