numerical_integration_6
37 Pages

numerical_integration_6

Course Number: PH 4390, Fall 2008

College/University: Mich Tech

Word Count: 2730

Rating:

Document Preview

Numerical Integration Monte Carlo Integration We want to Monte Carlo integrate the function 2 0 sin (x)dx = 1.0 In addition to the previous example we also compute the standard deviation, , and include it in the output. Then we want to vary (increase) the total number of MC evaluation points, N , and see how the results will change. A First Course in Computational Physics, PL Vries, Chp 4, ISBN...

Unformatted Document Excerpt
Coursehero >> Michigan >> Mich Tech >> PH 4390

Course Hero has millions of student submitted documents similar to the one
below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

Course Hero has millions of student submitted documents similar to the one below including study guides, practice problems, reference materials, practice exams, textbook help and tutor support.

Integration Monte Numerical Carlo Integration We want to Monte Carlo integrate the function 2 0 sin (x)dx = 1.0 In addition to the previous example we also compute the standard deviation, , and include it in the output. Then we want to vary (increase) the total number of MC evaluation points, N , and see how the results will change. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.1/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.2/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.2/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.2/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.2/30 Numerical Integration Monte Carlo Integration MC converges very slowly ! We need to go to very large values of N to get an overlap of the condence interval with the expected result. The standard deviation becomes smaller with increasing N. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.2/30 Numerical Integration Monte Carlo Integration The error in a Monte Carlo integration is fundamentally different to that discussed in previous numerical integration methods. When we compare it to the trapezoidal rule, we recall that there the error was found to be proportional to h2 . Thus when we decrease h in the trapezoidal scheme h like h we have actually halved the error. 2 h Making h smaller as however means making N 2 larger like N 2N . A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.3/30 Numerical Integration Monte Carlo Integration In contrast to the trapezoidal rule, the error in the Mone Carlo scheme is a probabilistic error. Including the calculation of allows us to say that our estimate is 68.3 % of the time within one standard deviation of the correct answer. Adding more points means that the average will improve following this 68.3 % criterion. Nevertheless in the case of a multidimensional integral, the error will remain to decrease according to the probabilistic nature of the averaging process and not according to the dimensionality of the integral ! A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.4/30 Numerical Integration Monte Carlo Integration Exactly this independence of dimensionality of the integration makes Mone Carlo integration an interesting technique for many problems. Considering again the trapezoidal analogue we would have to increase N by a factor of 2 for a two-dimensional integral in order to halve the error ( 2N in x-dimension and 2N in y -dimension ). In general we had to increase N in the trapezoidal technique as 2N for each additional dimension, thus some d-dimensional integral would require an increase d of N by a factor of 2 2 to halve the error. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.5/30 Numerical Integration Monte Carlo Integration We have already seen that the convergence of Mone Carlo integration is slow and goes as N . So to halve the error in MC schemes we had to increase N by a factor of 4. But since this happens regardless of the dimensionality of the integration, Monte Carlo becomes superior to e.g. the trapezoidal rule for dimensions greater than d = 4. So here is the actual strength of MC: faster convergence than any other scheme for high-dimensional integrals. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.6/30 Numerical Integration Monte Carlo Integration Multidimensional integrals are frequently encountered in many areas of Physics, especially in Statistical Mechanics. For example the expectation value of some macroscopic property u is given as phase space integral, i.e. u= R R ue kT dx3N dv 3N e kT dx3N dv 3N E E where dx3N , dv 3N means integration over the 3 components of position coordinates and velocity coordinates of each particle N . Already for a smaller value of N such a process is virtually impossible to all but the MC method. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.7/30 Numerical Integration Monte Carlo Integration As an example we want to carry out a multidemensional integration using the Monte Carlo method. Let us try to solve the 9 dimensional integral 1 1 1 .... 0 dax day daz dbx dby dbz dcx dcy dcz 0 0 (a+b)c We again want to monitor the evolution of the standard deviation with increasing numbers of evaluation points N . A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.8/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.9/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.9/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.9/30 Numerical Integration Monte Carlo Integration A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.9/30 Numerical Integration Monte Carlo Simulations Monte Carlo methods are also a common way to describe random processes. One of the earliest applications was to study scattering neturons in nuclear reactors. In particular, how much shielding is necessary to stop these neutrons. Many problems in physics are similar in spirit to random movement. Random processes are said to be stochastic . A rather well known example is the drunken sailor problem: After an extended bar tour and lots of partying the sailor is supposed to get back to the ship. But his choices of going into a certain direction happen arbitrarily and random. The question is how far will he have traveled after having taken N steps ? A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.10/30 Numerical Integration Monte Carlo Simulations Simulation of the random walk can be done on a square grid. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.11/30 Numerical Integration Monte Carlo Simulations A path similar to the one shown on the previous slide can be generated with the help of random numbers. Say the probability to move in any of the four directions is equal. A random number shall be created and if it is between 0 and 0.25 the movement shall be towards north. If the random number is between 0.25 and 0.5 the direction to move forward shall be east and so on... Of course, a single path wont tell us much and we need a number of pathes (repeat the simulation a number of times) to build up a distribution. After that we can state, how far the sailor is likely to travel. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.12/30 Numerical Integration Monte Carlo Simulations The histogram for a MC simulation of the Drunken Sailor Scenario of N = 100 steps. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.13/30 Numerical Integration Monte Carlo Simulations The drunken sailor example has some relevance to real physics problems. For example the mobility of an atom within the lattice points of a crystal could be described this way. Or also, the problem of diffusion could be addressed with MC simulations. In the case of diffusion, we know that at room temperature molecules in the air have velocities the on order of hundreds of meters per second. But if a bottle of perfume is opened at one end of the classroom, it takes several minutes to perceive the aroma at the other end of the classroom. Why ? A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.14/30 Numerical Integration Monte Carlo Simulations The explanation of the limited progress of the aroma of the perfume is that the molecules move around very much like the drunken sailor. There will be a large number of collisions with other molecules and each of these collisions will change the direction of the moving particle so that there will be little net gain in actually traveled distance. If we want to simulate that process we need to allow for random movement in all different directions. So we need a method to obtain uniformly distributed directional movements on the unit sphere. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.15/30 Numerical Integration Monte Carlo Simulations We can get the requested uniform distribution of solid angles by rst selecting from a uniform distribution of random variables on the interval [0, 2] followed by choosing a second random variable g from the interval [1, 1] and derive as = arccos (g) In order to simulate the diffusion process we would let the molecule move some given distance before some collision would occur. The collision would result in a random change of direction of continued movement. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.16/30 Numerical Integration Monte Carlo Simulations Of course the distance of unperturbed traveling is another critical component and should not be taken as a constant. But as a rst approximation we keep this as a xed constant and say it is the mean free path. So the scenario looks as follows: a molecule travels a mean free path distance, collides with another molecule, which leads to scattering in a random direction, travels another mean free path distance, and so on... We can simulate an example of real diffusion with a modied version of our drunken sailor MC program. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.17/30 Numerical Integration Monte Carlo Simulations For example, let us consider the diffusion of an aromatic molecule in air. The velocity of the molecule shall be 500 m and the mean free path shall be = 1m. We s want to calculate the distance < d > a molecule moves in one second. Our average shall include 100 trial runs of a single aromatic moecule. We nd a net distance traveled much smaller than the 500m that theoretically were possible. Our results show some clear accumulation around 19m. The importance of MC methods is not so much related to exact results but more to the fact that they can yield qualitatively valid results. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.18/30 Numerical Integration Monte Carlo Simulations The histogram for the diffusion of some particle as obtained from MC simulation. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.19/30 Numerical Integration Monte Carlo Simulations Although the results seem to be statistically solid, they still are physically invalid. It is true that the net displacement is considerably less than the free displacement. However, the estimated 19m per second still seem to be too much. From experience we would estimate the diffusion time for the molecules of a perfume traveling for about 20m to be on the order of minutes and not just one second. We could rerun the MC simulation with a different value for the mean free path. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.20/30 Numerical Integration Monte Carlo Simulations We might have actually found something more general. The net displacement in units of the mean free path length after 500 collisions. It would be interesting to see if there is a relationship between the magnitude of the net displacement and the number of collisions. One simple way of investigating whether there is some relation between the number of collisions and the net displacement is to modify our MC program and monitor the average net displacement as a function of the number of collisions. A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.21/30 Numerical Integration Monte Carlo Simulations If we do the mentioned modied MC simulation we nd the following behaviour ( 1000 different paths ). A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.22/30 Numerical Integration Monte Carlo Simulations The obtained plot is very smooth ! We could argue that there should be a relationship of the kind <d> Nq and the interesting question is what value do we get for the exponent q . A simple plot on logarithmic axes would yield a value of q 1. 2 We therefore have discovered that the displacement is proportional to the square root of N . A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.23/30 Numerical Integration Uniform Charge Distribution Let us assume we have a square region in the xy -plane, such that 1 x 1 and 1 y 1. In this square region there shall be a charge distribution . The charge shall be distributed uniformly. We are interested in the electrostatic potential at the point (xp , yp ). y (xp,yp) x A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.24/30 Numerical Integration Uniform Charge Distribution Since the charge is distributed uniformly every innitesimal volume element will contribute the same fraction of source charge and we can get the net effect via integration. y (xp,yp) x A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.25/30 Numerical Integration Uniform Charge Distribution Thus the electrostatic potential is (xp , yp ) = 40 1 1 1 1 dxdy (xp x)2 +(yp y)2 40 and we may take the prefactor y to be 1. (xp,yp) x A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.26/30 Numerical Integration Uniform Charge Distribution On the other hand we could make use of the set of orthogonal functions to expand the electrostatical potential, i.e. (r, ) = i=0 i (r)Pi (cos ) where Pi denotes Legendre polynomials in the unnormalized form. The coefcients i (r) cover the radial dependence, while the Legendre functions will do the angular dependence . A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.27/30 Numerical Integration Uniform Charge Distribution The initial unnormalized Legendre polynomials are P0 (x) P1 (x) P2 (x) P3 (x) P4 (x) P5 (x) = 1 = x 3x2 1 = 2 3 3x = 5x 2 35x4 30x2 +3 = 8 3 5 = 63x 70x +15x 8 A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.28/30 Numerical Integration Uniform Charge Distribution The orthogonality condition for these unnormalized functions is 0 Pm (cos )Pn (cos ) sin d = 2 2m+1 m,n With this condition we can derive the rule to obtain the radial expansion coefcients, i (r), for the electrostatic potential, (r, ) in the standard way, i.e. multiplication of the basic expansion with Pj (cos ) followed by integration over which yields, j (r) = 2j+1 0 2 (r, )Pj (cos ) sin d A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.29/30 Numerical Integration Uniform Charge Distribution With the just presented expression we can get a more general method to calculate the electrostatic potential. We just need to determine the expansion coefcients i (r) and we have xp = rp cos (p ) and yp = rp sin (p ). y (rp, p) x A First Course in Computational Physics, PL Vries, Chp 4, ISBN 0-471548689-3 PH4390: Computational Methods in Physics, MTU Fall 2006 p.30/30
MOST POPULAR MATERIALS FROM PH
MOST POPULAR MATERIALS FROM Mich Tech