#! /usr/bin/env python3 # def two_flies_square ( n ): #*****************************************************************************80 # ## two_flies_square(): average distance between two random flies on unit square. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 09 February 2023 # # Author: # # John Burkardt # # Input: # # integer n: the number of observations. # # Output: # # real d(n): the observed distances. # import numpy as np d = np.zeros ( n ) for i in range ( 0, n ): [ x1, y1 ] = np.random.random ( 2 ) [ x2, y2 ] = np.random.random ( 2 ) d[i] = np.linalg.norm ( [ x1 - x2, y1 - y2 ] ) return d def two_flies_square_test ( ): #*****************************************************************************80 # ## two_flies_square_test() tests two_flies_square(). # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 08 July 2022 # # Author: # # John Burkardt # import matplotlib.pyplot as plt import numpy as np min_exact = 0.0 mean_exact = ( np.sqrt ( 2.0 ) + 2.0 + 5.0 * np.log ( 1.0 + np.sqrt ( 2.0 ) ) ) / 15.0 max_exact = np.sqrt ( 2.0 ) var_exact = ( 69.0 - 4.0 * np.sqrt ( 2.0 ) \ - 10.0 * ( 2.0 + np.sqrt ( 2.0 ) ) * np.arcsinh ( 1.0 ) \ - 25.0 * ( np.arcsinh ( 1.0 ) )**2 ) / 225.0 std_exact = np.sqrt ( var_exact ) print ( '' ) print ( 'two_flies_square_test()' ) print ( ' Two flies land at random on a square plate of radius 1.' ) print ( ' What is the average distance between them?' ) print ( '' ) print ( ' n dmin dmean dmax dvar dstd ||dmean-exact||' ) print ( '' ) for n in [ 10, 100, 1000, 10000 ]: d = two_flies_square ( n ) d_min = np.min ( d ) d_mean = np.mean ( d ) d_max = np.max ( d ) d_var = np.var ( d ) d_std = np.std ( d ) print ( ' %5d %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f' \ % ( n, d_min, d_mean, d_max, d_var, d_std, np.abs ( d_mean - mean_exact ) ) ) print ( '' ) print ( ' Exact %6.4f %6.4f %6.4f %6.4f %6.4f' \ % ( min_exact, mean_exact, max_exact, var_exact, std_exact ) ) # # Make a histogram. # n = 10000 d = two_flies_square ( n ) plt.clf ( ) plt.hist ( d, bins = 50 ) plt.grid ( True ) plt.xlabel ( '<-- Distance -->' ) plt.ylabel ( '<-- Frequency -->' ) filename = 'two_flies_square_hist.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.close ( ) # # Compare observed data (converted to a density) # against theoretical PDF. # pdf = two_flies_square_pdf ( d ) plt.clf ( ) plt.hist ( d, bins = 20, rwidth = 0.95, density = True ) plt.plot ( d, pdf, 'r-', linewidth = 3 ) plt.grid ( True ) plt.xlabel ( '<-- Distance -->' ) plt.ylabel ( '<-- Probability -->' ) filename = 'two_flies_square_pdf.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.close ( ) return def two_flies_square_pdf ( d ): #*****************************************************************************80 # ## two_flies_square_pdf() evaluates the PDF for the two flies distance. # import numpy as np n = len ( d ) d.sort ( ) pdf = np.zeros ( n ) i1 = np.where ( d <= 1.0 ) pdf[i1] = 2 * d[i1] * ( d[i1]**2 - 4 * d[i1] + np.pi ) i2 = np.where ( 1.0 < d ) pdf[i2] = 2 * d[i2] * (4 * np.sqrt ( d[i2]**2 - 1 ) - ( d[i2]**2 + 2 - np.pi ) \ - 4.0 * np.arctan ( np.sqrt ( d[i2]**2 - 1 ) ) ) return pdf def timestamp ( ): #*****************************************************************************80 # ## timestamp() prints the date as a timestamp. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 21 August 2019 # # Author: # # John Burkardt # import time t = time.time ( ) print ( time.ctime ( t ) ) return if ( __name__ == '__main__' ): timestamp ( ) two_flies_square_test ( ) timestamp ( )