Simple Monte-carlo approximation of $\pi$ and integration using Matlab.

Question

Hi I am new to programming and the following task that my tutor want me to handle was a bit to hard to begin with. I hope someone can help me and I need the code in Matlab.

1. You are throwing darts at a circular dartboard (radius one half of a unit length) that are fixed at a square background (with the side one unit length). The dart hits this constellation at point $(x,y)$ that are a random variable. Want to use this idea to approximate the value of $\pi$.

2. I want to use the same idea to approximate $\int \limits _0 ^1 \sqrt[3]{1-x^3} {\rm d} x$.

Thanks in advance.

Answer 1

hint :

N = 2^16;
X = rand(2,N)*2-1;
norm2_X = sum(X.*X,1);
nb_inside_the_circle = sum(norm2_X <= 1);
pi_estimated = 3.14;
fprintf('pi_estimated = %f\n',pi_estimated);

add one more line to get an estimation for $\pi$.

Answer 2

Some comments, hints, and code that should help you finish this assignment.

1. Here is some R code for the first part. Your question is difficult for me to understand. I hope my comments are relevant.

The idea is to put a million points at random into a square of area 4 centered at the origin. Then we find the percentage of points in that data cloud (translated 'constellation"?) that also lie in the inscribed circle, Finally we convert that percentage to provide an estimate of $\pi.$

 m = 10^6;  x = runif(m, -1, 1);  y = runif(m, -1, 1)
 in.circle = (x^2 + y^2 <= 1)  # logical variable (values T and F)
 mean(in.circle)   # mean of logical vector is proportion of T's
 ## 0.785467       # proportion of pts in square that are also in circle
 4*mean(in.circle) # approx. value of pi
 ## 3.141868

The plot below shows 30,000 points similarly generated in the square (light grey) and the ones of them also in the circle (dark blue). [A million points would have made a dense "blob" of points---not a readable plot, but more points give a better numerical estimate of $\pi.$]

I am not a good enough Matlab programmer to write exemplary code, and, in any case, do not want to do the assignment for you. Together with the Answer by @user1952009 that gives some Matlab code, maybe you can write suitable code, and annotate it with explanations.

2. For the second part, you should first sketch the area represented by the integral. Then you might ponder what fraction of points in the unit square with vertices $(0,0)$ and $(1,1)$ lies beneath the curve being integrated. This part illustrates one of several methods of 'Monte Carlo' integration. It is sometimes called an 'acceptance-rejection' method because points below the curve are accepted and those above are rejected.

By contrast, the R code below approximates the area under the curve by 3000 rectangles, all of width 1/3000. The height each rectangle is the value of the function at the midpoint of its base. This is numerical integration (with no randomness) based on the definition of the Riemann integral. Your simulation in the second part should result in an answer close to the one below.

 m = 3000;  a = 0;  b = 1
 w = (b-a)/m                      # width of rectangles
 g = seq(a+w/2, b-w/2, length=m)  # grid of centers
 h = (1 - g^3)^(1/3)              # heights at centers of rectangles
 sum(w * h)
 ## 0.8833213                     # deterministic approx of integral

Note: For a given number of iterations $m,$ an integral over an interval on the real line is frequently more accurately evaluated by Riemann approximation than by Monte Carlo integration. However, for multiple-dimensional integration, Monte Carlo methods are frequently superior.