Montecarlo Methods to Efficently Compute Volume of Solids Kamil Czapla1 , Mariusz PleszczyΕ„ski2 1 Academic High School of Silesian University of Technologies, Rudzka 13, Rybnik, 44–200, Poland 2 Faculty of Applied Mathematics, Silesian University of Technology, Kaszubska 23, Gliwice, 44–100, Poland Abstract In this work we present and compare chosen Monte Carlo methods, used to calculate solid’s volume. Such solids are defined with a certain surface 𝜚(πœ™, πœƒ) wich results in classical methods being hard or impossible to apply. In designating those volume, due to the desire to shorten calculation time, we will apply parallel calculations. Keywords Double integral, Monte Carlo methods, parallel computations 1. Introduction above and below respectively. If the domain D s normal (e.g) with respect to the x-axis, i.e. can be described by In classic task of calculating solid’s volume we use double dependencies: integral (or triple integral), often utilizing variable’s con- version. Most often they are polar, cylindrical or spherical coordinates. The use of this approach most often occurs, when surfaces limiting given solid are given in open form (2) 𝑧 = 𝑓 (π‘₯, 𝑦)en however, whether for curves in a plane or for surfaces in space, such limitation are given different hen the volume (1) can be described using (2), in the (non-classical) way. For curve this could be for example: form: parametric or polar form. In case of a surface from 𝑅3 surface, the situation may look similar, but then there is an additional variable. Wanting to calculate area of a given flat surface or volume of a surface in space, appro- priate formulas should be designated. Authors discussed their use in [1] for flat areas and in [2] for spatial areas. Especially in the 𝑅3 space, where surface limiting given (3) volume is dependent from 𝜚(πœ™, πœƒ) suitable integral de- scribing this volume may be hard or even impossible to To describe such domain also can be used polar coordi- determine analytically. In this work we propose the use nates, cylindrical or spherical. Remembering the Jacobian of certain approximate methods on calculating numer- determinant of given variables swap and describing right ically given volume, independent from the form of the Ddomain (for flat surfaces) or Space (in case of describing funcion 𝜚(πœ™, πœƒ) solid, e.g. using spherical coordinates) we get appropri- ate, known to us equations for calculating volume of the solid with given coordinates. We want to focus on the 2. Domains in 𝑅3 space case, in which area limiting given solid is defined by a In which the volume V is described by the integral: function 𝜚(πœ™, πœƒ) this case any point of this surface is de- fined by three variables. The first of them is the distance 𝜚 from this point to the beginning of the coordinate sys- tem. The second is an angle πœ™ between positive Ox is and a view, on a Oxy face of a segment joining the begin- (1) ning of a coordinate system and given point. The third where continuous in D domain functions 𝑓1 and 𝑓2 one is an angle πœƒ between positive Oz is and a segment, are limiting the space above and below D domain, from joining the beginning of a coordinate system and given point (view Fig. 2). If for certain angles πœ™ and for certain ICYRIME 2021 @ International Conference of Yearly Reports on angles πœƒ e define non-negative function 𝜚(πœ™, πœƒ) then all Informatics Mathematics and Engineering, online, July 9, 2021 points received this way will define a certain surface in " kamil.czapla@alorybnik.polsl.pl (K. Czapla); mariusz.pleszczynski@polsl.pl (M. PleszczyΕ„ski) 𝑅3 space. Β© 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). It turns out that then the volume of a solid defined this CEUR Workshop Proceedings http://ceur-ws.org ISSN 1613-0073 CEUR Workshop Proceedings (CEUR-WS.org) way can e described by integral [2]: 1 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 Figure 2: Monte Carlo method illustration Figure 1: Spherical coordinates (5) where c β©Ύ min g(x), d ≀ max f(x) The next step is to draw n points belonging to the rectangle R. If this choice is random (from uniform distribution) and independent, then these points should cover the rectangle R evenly. Therefore, the ratio of the number of m points lying (4) between the graphs of the functions f and g to the number In this connection that the form of the function 𝜚 can of all points n is the same as the ratio of the fields: S– lead to integrals that do not express themselves with ele- between the graphs of the functions f and g |R| – area of mentary functions, thereby the values of them cannot be the rectangle R We have: determined precisely, therefore we will use some numeric methods to calculate approximate volume. 3. Monte Carlo methods (6) Monte Carlo methods appeared (most frequently this is assumed, although the theory describing these methods In Fig. 2 we present these dependencies graphically, appeared the most frequently) in 1949, after the publica- with green color marking the points that went to the tion of the work [3], and allowed their creators to solve a area of interest, and red color those that did not go to computationally complex problem, that emerged during this area. the work on the construction of the atomic bomb. The Monte Carlo methods family provides an approximate 3.1.1. Example 1 solution to a wide class of computationally complex prob- We will show how it is possible to use this method to lems, including numerical integration. We will focus on find the approximate value of the number 𝑒. For this we this group of problems. The computational complexity consider the integral: and error estimation obtained by this method show [4] that as the dimension increases, Monte Carlo methods become more and more useful in relation to classical methods (e.g. to quadrature based methods). We will deal with the double integral, but the ideas of the Monte (7) Carlo methods, for the sake of clarity, will be shown on Knowing its exact value, we can estimate the value of the example of the integral. the number 𝑒 rising the formula (4). Let’s take a rectangle: 3.1. Monte Carlo methods no. 1 In classical terms, the domain between the graphs of the (8) functions f(x) and g(x)we assume that 𝑓 (π‘₯) β©Ύ 𝑔(π‘₯) for and let’s conduct a series of experiments involving 𝑏 β©Ύ π‘₯ β©Ύ π‘Ž can be approximated as follows: enter the the drawing and verification of 𝑛 points on the rectan- first graphs of the functions f and g into the rectangle R: gle R For each of the n values presented in Table 1, we 2 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 Figure 4: Monte Carlo method no.I – dependence of the Ξ” error on the number of ndraws Figure 3: Monte Carlo method no.I conducted 50 draws (we want to check the stability of is divided evenly, that is Ξ”π‘₯𝑖 = Ξ”π‘₯ = π‘βˆ’π‘Ž . Thanks to 𝑛 the method) and determined the following values: min – this, formula (5) will take the form: smallest area value, max – largest area value, S- average area value (as the arithmetic mean of the 50 results ob- tained), 𝜎 – standard deviation of the obtained results, Ξ” – calculation error in the formula Ξ” = |𝑆 βˆ’ 𝑆 Β― | the area (10) S is the exact value of the area, and thus the number |e| and 𝑆 Β― s the obtained value. where πœ‰π‘– , 𝑖 1, 2, . . . n are the drawn points of the The conducted experiments show that this method is interval [a, b. On this basis, the field described by formula stable and gives a good approximation. The Monte Carlo (4) will take the form: error theory says that the error obtained is proportional to the number βˆšπ‘π‘› , where 𝑐 < 0 is a constant and 𝑛 s the number of points drawn. Figure 3 shows that this relationship is preserved (in this figure, 3 has the value 0.5) (11) 3.2. Monte Carlo methods no. 2 3.2.1. Example In this method, we have a slightly different approach to determining the value of the integral previously ap- Now let’s do the same experiment as in Example 3.2.1. proximated by the formula (4). This time we refer to the Table 2 shows the results of these experiments, and Fig. definition of the Riemann definite integral: 4 we present the equivalent of Fig. 3, where the constant c as the value 0.1. Intuitively, it is difficult to estimate which of the for- mulas (4) or (6) is represented by the greater error. It (9) turns out that the formula (6) is more useful, and the where πœ‰ , 𝑖 = 1, 2, ...𝑛 are any points up to the subin- results summarized in Tables 1 and 2 confirm this theory 𝑖 terval [π‘₯𝑖 , π‘₯𝑖+1 ],𝑖 = 0, 1, . . . 𝑛 of the interval [π‘Ž, 𝑏] These sub-intervals are determined by a certain natural divi- 3.3. Other Monte Carlo methods sion, Ξ”π‘₯𝑖 is the length of such a sub-interval. If we run In addition to the many advantages of Monte Carlo meth- ndraws again (the assumptions for the drawing are the ods, such as simplicity of implementation, regardless of same as before), we can assume that the interval [π‘Ž, 𝑏] the dimension of the task, or the non-growing of errors 3 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 as the dimension increases, there are also disadvantages 3.3.2. Importance sampling of this method. The most frequently mentioned is the This method [5] is based on the assumption that more value of the obtained errors, which is proportional to √1𝑛 points should be drawn in more significant parts of the , where nis the number of points drawn. To increase the interval [π‘Ž, 𝑏]. To do this, one needs to find a function error by one row, e.g. from 1/10 to 1/100 , the number of 𝑝(π‘₯) which in the interval[π‘Ž, 𝑏] is the probability density, points drawn should be increased from 100 to 104 . The i.e. 𝑝(π‘₯) > 0 methods discussed by us in this part of the work want Of course, there are infinitely many such functions, to counteract this problem. There are many different and its choice has an impact on the resulting error. Most approaches to such a task, we will cover the three most often, choosing the best such function is impossible, be- commonly used. cause it is related to the value of the integral (which is unknown). In practice, the approach is to select the 3.3.1. Separating of the main part function 𝑝(π‘₯)in such a way that, similar to the previous This method consists in finding a certain function which example, it is close to the function 𝑓 (π‘₯) i.e. that 𝑓𝑝(π‘₯) (π‘₯) β‰ˆπ‘ is close to the integral function, but for which we can find where c is a constant. After finding such a function, we the exact value of the integral. To refer to our example can use the formula: of finding the approximate value of e where we integrate the function 𝑓 (π‘₯ = 𝑒π‘₯ + 1), we can use the well-known Maclaurin series expansion of the function 𝑒π‘₯ (17) which leads to the formula for the approximate value (12) of this integral: Also note that: (18) (13) where πœ‰ = 1, 2, . . . 𝑛 are random interval points [π‘Ž, 𝑏] In the discussed example, using the previously found Since the difference between the functions f(x and g(x) function 𝑔(π‘₯) = 2 + π‘₯ + 0.5π‘₯2 is determined, so: we can take: 𝑝(π‘₯) = 38 𝑔(π‘₯) (19) (14) So we will get the following formula which allows us to find the approximate value of the number 𝑒 and in the discussed example, this formula takes the form: (20) (15) where πœ‰ = 1, 2, . . . 𝑛 are random distribution 𝑝(π‘₯) which leads to the formula for the approximate value (You can read about the drawing of points from a given of this integral: distribution in e.g. [6, 7, 8]). 3.3.3. Stratified sampling This approach [9] is similar to the previous one. This time we divide the integration interval [a,b into k subin- tervals and randomize π‘›π‘˜ points in each of them. The (16) numbersπ‘›π‘˜ e chosen so that their sum is 𝑛 and that their size is proportional to the rate of change of the integral where πœ‰ = 1, 2, . . . 𝑛 are random interval points [π‘Ž, 𝑏] value of the function 𝑓 (π‘₯) Based on this approach, the from uniform distribution. sought integral can be approximated by the formula: 4 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 Figure 5: Comparison of other Monte Carlo methods (21) (𝑗) where πœ‰π‘– is the i-th number drawn from the j-th subinterval, and 𝑑𝑗 is the length of the jth sub-interval. 4. Monte Carlo methods for If the interval [0, 1] we will divide it into 5 equal parts double integrals and choose the next sub-ranges accordingly 1/25 , 3/25 , 5/25 , 7/25 and 9/25 the number 𝑛 of all drawn points, The advantages of the Monte Carlo method, highlighted then the formula (9) adjusted to approximate the value by us previously, are perfect for double integrals, both of the number e will take the form: the simplicity of implementation (for both classical meth- ods and their improvements) and the estimated error are transferred to double integrals (in general, to multiple integrals). As we mentioned in the introduction, we want to use Monte Carlo methods for integrals (3). This inte- gral is formed when the surface bounding a given spatial domain is given by the function𝜚(πœ™, πœƒ) [2] authors dis- cussed the formula (3) and applied it to several surfaces, incl. peanut, bluebell, shell and exotic fruit. In this paper, (22) we will use three approaches to the problem of determin- ing these integrals using the Monte Carlo method I, the 3.3.4. Result for other Monte Carlo methods Monte Carlo method II and parallel calculations used in the Monte Carlo method II, the aim of which is to shorten As for the previous methods, for the same values of n we the time needed to determine a given integral. Moreover, carry out the process of approximating the value of the as some integrals (3) are reduced to a single integral in number e using the formulas (7), (8) and (9). In tables the above-mentioned surfaces, we will consider domains 3, 4 and 5 we will collect the appropriate results and in both on the plane and in space. Parallel computations drawing 5 we illustrate the absolute error Ξ” similarly as will be performed using a multithreaded process, and the we did in the drawings 3 and 4. aaabbbccc process itself is very simple. As the points The conducted research confirms that, compared to drawn are independent and do not affect each other, as the best known Monte Carlo method I, all other methods well as the calculations of the function values in these are more effective. For the same number n of drawn points, it is enough to divide the points, they give a better (statistically) approximation of number of draws n into k equal parts and each thread the value of e will perform its calculations in parallel. After making the calculations by all threads, the results obtained are summed up to the final result. Note that the formula (6) can be implemented in parallel – the number n of all points is divided into k equal parts, where k is the number of threads. Each thread will do its own sum: 5 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 Figure 6: Comparison of other Monte Carlo methods Figure 7: Bell area (23) where β„Ž(π‘₯) = 𝑓 (π‘₯) βˆ’ 𝑔(π‘₯) are random points from This time, the equivalent of formula (4) should be a the interval [π‘Ž, 𝑏] And then we will create the final sum: formula adapted to the size, so the volume V of the sought solid would be approximate by the formula: (24) This process is shown in Fig. 6. (27) The calculations were carried out on a computer with a 12-thread Intel Core i7-3930k CPU 3.20 GHz processor, where, like before, m is the number of points that hit 16 Gb RAM and a 64 bit Windows 7 system. Programs the inside of the bell, and n is the number of all drawn that implement all Monte Carlo methods were written points of the cuboid 𝐢 However, using the spaces per- in Mathematica 12.2 (including those for parallel com- formed in the work [2] the volume of the bell is described puting) [10]. We will now present the application of the by the single integral: proposed Monte Carlo methods on the example of two surfaces introduced and discussed in [2] – bell and exotic flower (28) So we have a situation where (without taking into ac- 4.1. Bell Area count the unit) the volume of the bell equals (numerically) This area is represented by the relationship: 𝜚(πœ™, πœƒ) = the area under the curve: π‘Ž + 𝑠𝑖𝑛(2πœƒ) + π‘π‘œπ‘ (3πœƒ) and its graph, for a= 2 is shown in Fig. 7. As we showed in [2] the volume of the solid limited (29) by this area is: for x in [0, πœ‹], which means that it is enough to draw n points of the rectangle R. (25) (30) which for a= 2 gives the value 39.382572945... (the number 12 can be determined on the basis in Fig. 8 which shows the graph of the function 𝑓 (π‘₯) along with 4.1.1. Monte Carlo method no. I the drawn 𝑛 = 103 points) and use the formula (4). According to the idea of this method, the surface of the After conducting the experiment consisting in draw- bell will be entered in a cuboid 𝐢 ing 𝑛 = 9! points of the rectangle 𝑅 e obtained the following approximate result of the bell volume value: 𝑉 = 39.3611 what the error corresponds to Ξ” = 0.0214 (26) This result, which will be important when comparing the 6 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 Figure 8: Area under the curve f(x) corresponding to the Figure 9: Dependence of t time on the number of k threads volume of the bell that we will use parallel computations with the use of π‘˜ = 2, 3, ..., 10 threads. The results of these experiments are collected in Table 8, we will pay special attention to the average time ¯𝑑 – its change is to prove the usefulness of parallel computations. As indicated by the data in table 8, the introduction of calculations, on the one hand does not worsen the obtained volume results, and on the other hand reduces the time needed to obtain this result. methods, was achieved in 𝑑 = 3.19802 seconds. Pro- The comparison of the dependence of this time t to the ceeding similarly to chapter 3 we then conducted a series number of threads k is shown in Fig. 9. of 50 such experiments, and their results (designations analogous to the previous ones with an additional value of the calculation time ¯𝑑 are presented in Table 6. 4.2. Exotic fruit surface This area will be represented by the relationship: 4.1.2. Monte Carlo method no.II As before, we will go straight to the calculation of the area (31) under the curve (10) and use the formula (6). Therefore, we will select 𝑛 = 9! points from the interval [0, πœ‹]and 0 ≀ πœ™ ≀ 2πœ‹, 0 ≀ πœƒ ≀ πœ™ and π‘Ž β‰₯ 1, and its graph for a searched area (and thus the volume of the bell) we a= 1 is shown in Fig. 10. will determine approximately. This time we got an area As we showed in [2], the volume of a solid limited by (volume) equal to 39.3762, which corresponds to error this area is: Ξ” = 0.0064. We obtained this result in 𝑑 = 1.903 sec- onds. After conducting 50 such experiments, we obtained the results summarized in Tab. 7. (32) 4.1.3. Parallel Monte Carlo method no.II what for a=1 gives value of 7.120943348. Monte Carlo methods are perfect for parallelization. We can do it because the independent draws of n points from 4.2.1. Monte Carlo method no. I the interval [π‘Ž, 𝑏] is equivalent to k draws of n/k points from this interval (we sum the k set obtained in this way, This time we will refer to the formula (4) which says that obtaining one n-elementary set). Let us now carry out the surface of an exotic fruit is limited by an area whose experiments analogous to the above experiments, except 7 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 Figure 11: Exotic fruit surface. The volume limited by the area of 𝑓 (π‘₯, 𝑦) corresponding to the volume of an exotic fruit Figure 10: Exotic fruit surface volume corresponds to the volume of the area limited from above by the area 𝑓 (π‘₯, 𝑦) above the rectangle (33) where (34) which means that it is enough to draw n points of the cuboid 𝐢 (35) (the number 8/3 can be determined on the basis in Fig. 11, which shows the graph of the function 𝑓 (π‘₯, 𝑦) with the drawn 𝑛 = 104 points) and use the formula (6), which this time will have to change by one form dimension: (37) where πœŒπ‘– , 𝑖 = 1, 2, ...𝑛 are the points of the rectangle on which the integration takes place, we will get the (36) results (we have carried out 50 such experiments again) where m is the number of points that hit the inside of which are summarized in Tab. 10. the fruit (and thus are in the cuboid C below the surface 𝑓 (π‘₯, 𝑦)) 4.2.3. Parallel Monte Carlo method After conducting a series of 50 experiments described above, assuming n= 9!, we obtained the results presented In this method we will act in a similar way to what we in table Tab. 9. did in subsection 4.1.3. Again we will draw 𝑛 parallel points, using k threads, each of which will draw 𝑛/π‘˜ points. As before, we will use π‘˜ = 2, 3, ..., 10 threads, 4.2.2. Monte Carlo method no. II and the results of these experiments are presented in Tab. Proceeding as before, i.e. drawing n= 9! points and using 11. Again, it turned out that multi-threading does not the formula (11), after adjusting it to the dimension of deteriorate the quality of the solution, while reducing the task: the computation time. A graphical interpretation of the 8 Kamil Czapla et al. CEUR Workshop Proceedings 1–9 dependence of this time on the number of threads used [9] The stochastic collocation monte carlo sampler: is shown in Fig. 12. Highly efficient sampling from ’expensive’ distribu- tions,Quantitative Finance, 19(2), 339-356. [10] PoΕ‚ap D., WΓ³zniak M., Napoli C., Tramontana E., Is 5. Conclusions Swarm Intelligence Able to Create Mazes? (2015) International Journal of Electronics and Telecom- The authors in [2] introduced a certain formula for de- munications, 61 (4), pp. 305 - 310, doi: 10.1515/eletel- termining the volume of solids. They also calculated 2015-0039. the volume of the selected solids.Since the formula is [11] PoΕ‚ap D., WΓ³zniak M., Borowik G., Napoli C., A non-standard (in this work it is formula (4)) the authors First Attempt to Cloud-Based User Verification in decided to verify it. for this purpose, we used selected Distributed System (2015) Proceedings - 2015 Asia- Monte Carlo methods that are well suited for this purpose. Pacific Conference on Computer-Aided System En- We discussed 5 such methods, and to verify the formula gineering, APCASE 2015, art. no. 7287024, pp. 226 - (4) we used 2 of them and an additional sixth method, 231, doi: 10.1109/APCASE.2015.47. which is a parallel version of the second method. The [12] Steinbrecher, G., and Shaw, W. T. (2008). Quantile conducted research has shown that Monte Carlo methods mechanics. European journal of applied mathemat- are effective in this type of tasks, that their modifications ics, 19(2), 87-112. are more useful, that formula (4) is correct and that the [13] Botev, Z., and Ridder, A. (2017). Variance reduction. disadvantage of Monte Carlo methods can be effectively Wiley statsRef: Statistics reference online, 1-6. dealt with. The introduced parallel method does not have [14] Morrison, Michael, Cliff Hastings, and Kelvin a negative impact on the quality, and at the same time Mischo. "Hands-on start to wolfram mathematica." shortens the calculation time. (2015). References [1] K. Czapla, M. PleszczyΕ„ski, Some formulas for de- termining areas and volumes, part I: areas in space 𝑅2 , MINUT 3 (2021) 208–228. In Polish [2] K. Czapla, M. PleszczyΕ„ski, Some formulas for de- termining areas and volumes, part II: areas in space 𝑅3 , MINUT 3 (2021) 229–244. In Polish. [3] N. Metropolis, S. Ulam, The monte carlo method, Journal of the American Statistical Association 44 (1949) 335–341. [4] N. P. Buslenko, D. I. Golenko, I. M. Sobol, W. G. Sragowicz, J. A. Szrejder, The Monte Carlo Method, PWN, Warsaw, 1967. In Polish. [5] Capizzi G., Sciuto G.L., Napoli C., Shikler R., Woz- niak M., Optimizing the organic solar cell manu- facturing process by means of AFM measurements and neural networks (2018) Energies, 11 (5), doi: 10.3390/en11051221. [6] Kloek, T., and Van Dijk, H. K. (1978). Bayesian esti- mates of equation system parameters: an applica- tion of integration by Monte Carlo. Econometrica: Journal of the Econometric Society, 1-19. [7] Devroye, Luc. "Sample-based non-uniform random variate generation." Proceedings of the 18th confer- ence on Winter simulation. 1986. [8] WrΓ³bel M., Starczewski J.T., Napoli C., Handwrit- ing recognition with extraction of letter fragments (2017) Lecture Notes in Computer Science (includ- ing subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 10246 LNAI, pp. 183 - 192, doi: 10.1007/978-3-319-59060-8_18. 9