Monte Carlo Methods - Massachusetts Institute of …

1 Chapter 17 Monte Carlo Methods59 A taste of Monte Carlo methodMonte Carlo Methods is a class of numerical Methods that relies on random sampling. For example, the followingMonte Carlo method calculates the value of :1. Uniformly scatter some points over a unit square[0,1] [0,1], as in Figure ??.2. For each point, determine whether it lies inside the unit circle, the red region in Figure ??.3. The percentage of points inside the unit circle is an estimate of the ratio of the red area and the area of the square,which is /4. Multiply the percentage by 4 to estimate .The following Matlab script performs the Monte Carlo calculation:1n = 1000000; # number of Monte Carlo samples2x = rand(n, 1); # sample the input random variable x3y = rand(n, 1); # sample the input random variable y4isInside = (x.)

2 2 + y. 2 < 1); # is the point inside a unit circle?5percentage = sum(isInside) / n; # compute statistics: the inside percentage6piEstimate = percentage*4 Exercise each value ofn=100,n=10000 andn=1000000, run the script 3 times. How accurate is theestimatedpi?This example represent a general procedure of Monte Carlo Methods : First, the input random variables (xandy)are sampled. Second, for each sample, a calculation is performed to obtain the outputs (whether the point is inside ornot). Due to the randomness in the inputs, the outputs are also random variables. Finally, the statistics of the outputrandom variables (the precentage of points inside the circle) are computed, which estimates the a Matlab script to use Monte Carlo method to estimate thevolume of a 3-dimensional a Matlab script to use Monte Carlo method to estimate thevolume of a 10-dimensional Monte Carlo method in Engineering: Colloid thrusterIn many engineering problems, the inputs are inheriently random.

3 As an example of Monte Carlo method for theseengineering applications , we study a space propulsion device, the colloid thruster. It use electrostatic acceleration ofcharged particles for propulsion. The charged particles are produced by an electrospray process, and have randominitial velocity trajectory of the particle inside the electrical field 0 x Lis governed by the set of ODEsdxdt=vx,drdt=vr,dvxdt=Ax(x),dvrdt=0, (141)whereAx=emd (x)dxis the acceleration induced by the electrical field. Here we use the polynomial approximationAx=(Ax0(1 3x2+2x3)x L0x>L(142)The initial position of the charged particle is atx(0) =r(0) =0. The initial velocity isvx(0) =V0cos 0,vr(0) =V0sin 0(143)In this equation,V0is the initial speed of the charged particle. It is assumed tobe a fixed value in this section.)

4 0is theinitial angle, and israndom. Here, we assume that 0is auniform random variableU(0, max). We are interestedin thestatistical distributionof velocity(vx,vr)after the charged particle leave the electrical field tox> simple simulation of charge particle acceleration can be performed using the following Matlab function. Forgiven initial velocityV0, initial angle 0, Coulumb accelerationAx, length of accelerationLand time step size t, thefunction returns the final [vx, vr] = thruster(V0, a0, Ax0, L, dt)2% initial condition3x = 0; r = 0;4vx = V0*cos(a0); vr = V0*sin(a0);5% time integration using Forward Euler6while(x L)7x = x + dt*vx;8r = r + dt*vr;9vx = vx + dt*Ax0*(1 - 3*x. 2 + 2*x. 3);10endIn our problem, there are large amounts of charged particles, and the initial angle of each charged particle is role of probabilistic Methods is to quantify the impact of this type of randomness on properties of interest ( theterminal velocity).

5 The results of the probabilistic analysis can take many forms depending on the specific the example of the thurster where the terminal velocity iscritical to the performance and efficiency of the thruster,the following information might be desired from a probabilistic analysis: The distribution of the terminal velocityvx,vrthat would be observed in the population of charged particles. The probability thatvxorvris above some critical value ( , indicating the particlewill hit part of the thruster). Instead of determining the entire distribution of terminal velocity, sometimes knowing the mean values vx, vrissufficient, for calculating the thrust. To have some indication of the variability ofvxandvywithout requiring accurate estimation of the entire distribu-tion, the standard deviation, vxand vy, can be Monte Carlo method is based on the idea of taking a small, randomly-drawn sample from a population andestimating the desired outputs from this sample.

6 For the outputs described above, this would involve: Replacing the distribution ofvx,vrthat would be observed over the entire population of particles with the distribu-tion ( histogram) of those observed in the random Replacing the probability thatvxorvris above a critical value for the entire population of charged particles withthe fraction of particles in the random sample that havevxorvrgreater than the critical value. Replacing the mean value ofvxandvrfor the entire population with the mean value of the random sample. Replacing the standard deviation ofvxandvrfor the entire population with the standard deviation of this exactly what is done in the field of statistics, theanalysis of the Monte Carlo method is a direct applicationof summary, the Monte Carlo method involves essentially three steps:1.

7 Generate a random sample of the input parameters according to the (assumed) distributions of the Analyze (deterministically) each set of inputs in the Estimate the desired probabilistic outputs, and the uncertainty in these outputs, using the random following Matlab code performs the Monte Carlo simulation for our thruster1% Deterministic (non-random) parameters2V0 = ; Ax = ; L = ; dt = ;3a0 Max = 60.*pi/ 180.;4% 1. Generate a random sample of the input parameters5a0MC =rand(10000, 1)*a0 Max;6% Array to store Monte Carlo outputs7vxMC = [];8vrMC = [];9% 2. Analyze (deterministically) each set of inputs in the sample10fori = 1:1000011[vx, vr] = thruster(V0, a0MC(i), Ax, L, dt);12vxMC = [vxMC; vx];13vrMC = [vrMC; vr];14end15% 3. Estimate the desired probabilistic outputs16% histogram17figure;hist(vxMC);18figure;hi st(vrMC);19% mean20muVx =mean(vxMC)21muVr =mean(vrMC)22% standard deviation23sigmaVx =std(vxMC)24sigmaVr =std(vrMC)Exercise Monte Carlo method to estimate the probabilityP(vr>C).

8 61 Uniform and Non-Uniform Random VariablesIn the previous examples, the random input parameters have uniform distribution. A uniform distribution is definedby the two parameters,aand Ib, which are the minimum and maximum values the random variable can possibly the interval(a,b), all values are equally probable. The distribution is oftenabbreviatedU(a,b): Its probabilitydensity function (pdf) isf(x) =(1b aa x b0x<aorx>b(144)It s cumulative distribution function (cdf), the integralof its pdf, is100F(x) = x ab aa x b0x<a1x>b(145)Uniform random variable is special in Monte Carlo Methods and in computation most psuedo random numbergenerators are designed to generate uniform random numbers. In MATLAB, for example, the following commandgenerates anmbymarray ofU(0,1)uniform random (m,n);To generate anU(a,b)uniformrandom numbers, one can simply scale theU(0,1)random numbers byx=rand(m,n)*(b-a)+a;Almost all otherlanguages used for scientific computation have similar random number is themean, varianceandstandard deviationof aU(a,b)random variable?)

9 Non-uniform distributions are those whose probability density functions are not constant. Several simple but im-portant non-uniform distributions are Triangular distribution. It is characterized by three parametersa,b,c. The probability density function isf(x) = 2(x a)(b a)(c a)a x c2(b x)(b a)(c a)c x b0x<aorx>b(146) Exponential distribution. It is characterized by a singleparameter . The probability density function isf(x) =( e xx 00x<0(147) Normal distribution, also known as Gaussian Sampling non-uniform random Rejection methodRejection for triangular Transformation methodUniform sampling for square root of 0. Triangular distributionInverse cummulative density to exponential random variable10162 Risk Computing probability of Charged particle breakup and errosion in a Colloid thruster63 Error Estimation for Monte Carlo Law of Large Number of Central Limit Error Estimator via Central Limit Application to Mean and Application to Risk Error Estimator via Bootstrap64 Importance Sampling65 Linear Sensitivity Method (Delta Method))