Page - 84 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

84 3 Computing Integrals 3.7.3 MonteCarloIntegrationforComplex-ShapedDomains Repeateduseofone-dimensional integration rules tohandledoubleand triple inte- grals constitute aworking strategy only if the integration domain is a rectangle or box. For any other shape of domain, completely differentmethodsmust be used. A common approach for two- and three-dimensional domains is to divide the do- main intomanysmall trianglesor tetrahedraandusenumerical integrationmethods for each triangle or tetrahedron. The overall algorithmand implementation is too complicated to be addressed in this book. Instead,we shall employ an alternative, very simple and general method, calledMonte Carlo integration. It can be im- plemented in half a page of code, but requires orders ofmagnitudemore function evaluations indouble integralscompared to themidpoint rule. However,MonteCarlo integration ismuchmore computationally efficient than themidpoint rulewhencomputinghigher-dimensional integrals inmore than three variables over hypercube domains. Our ideas for double and triple integrals can easilybegeneralized tohandlean integral inmvariables.Amidpoint formula then involvesm sums. Withn cells in each coordinate direction, the formula requires nm functionevaluations. That is, the computationalwork explodesas an exponen- tial function of the number of space dimensions. MonteCarlo integration, on the other hand, does not suffer from this explosion of computationalwork and is the preferredmethod for computing higher-dimensional integrals. So, itmakes sense in a chapter on numerical integration to addressMonte Carlo methods, both for handlingcomplexdomainsand forhandling integralswithmanyvariables. TheMonteCarlo integrationalgorithm The ideaofMonteCarlo integrationofRb a f.x/dx is to use themean-value theorem fromcalculus, which states that the integral Rb a f.x/dx equalsthelengthoftheintegrationdomain,hereb a, timesthe averagevalueoff , Nf , in Œa;b . The averagevalue canbe computedby sampling f at a set of random points inside the domain and take themean of the function values. In higher dimensions, an integral is estimated as the area/volume of the domain times theaveragevalue,andagainonecanevaluate the integrandatasetof randompoints in thedomainandcompute themeanvalueof thoseevaluations. Letus introducesomequantities tohelpusmake thespecificationof the integra- tionalgorithmmoreprecise. Supposewehavesome two-dimensional integral Z ˝ f.x;y/dxdx; where˝ is a two-dimensionaldomaindefinedviaahelp functiong.x;y/: ˝ Df.x;y/jg.x;y/ 0g That is, points .x;y/ for which g.x;y/ 0 lie inside˝, and points for which g.x;y/ <˝ are outside˝. The boundaryof the domain@˝ is givenby the im- plicit curveg.x;y/D0. Such formulationsofgeometrieshavebeenverycommon during the last coupleofdecades, andonerefers tog asa level-set functionand the