To do this, we define a square, whose area we know, that bounds the circle. So, let the side length of the square equal the diameter of the circle.
We conceptually throw darts randomly at the square.
The fraction of the darts that land in the circle multiplied by the area of the square is the area of the circle:
$$A_{circle} = A_{square}\frac{n_{in}}{n_{tot}},$$
where $n_{in}$ is the number of throws that are inside the circle, and $n_{tot}$ is the total number of throws.
using Plots
using Random
using LaTeXStrings
using Polynomials
using CurveFit
r = 1.0 # circle radius
Asquare = (2*r)^2 # area of square enclosing circle
n = 100000 # number of tries
rx = rand(n)*2*r .- r # -1 to 1; x locations
ry = rand(n)*2*r .- r # -1 to 1; y locations
rr = sqrt.(rx.^2 .+ ry.^2) # radial location
nhits = length(rr[rr .< r]) # number inside the circle
I = nhits/n * Asquare # integral of circle
Aexact = pi*r^2
println("A_mc = $I")
println("A_exact = $Aexact")
A_mc = 3.14628 A_exact = 3.141592653589793
i_n = zeros(n)
i_n[rr .< r] .= 1
itry = 1:n
Areas = cumsum(i_n)./itry*Asquare
errs = abs.((Areas.-Aexact)/Aexact)
plot(log10.(itry), log10.(errs), label="Relative Error")
plot!([0,log10(n)], [log10(errs[1]), log10(errs[1])-0.5*log10(n)], label="Slope=-1/2")
plot!(xlabel=L"\log_{10}(n_{tot})")
plot!(ylabel=L"\log_{10}(\epsilon_r)")
plot!(legend_foreground_color=nothing)
r = 1.0 # circle radius
As = (2*r).^2 # area of square enclosing circle
n = 10000 # number of tries
rx = rand(n)*2*r .- r # -1 to 1; x locations
ry = rand(n)*2*r .- r # -1 to 1; y locations
rr = sqrt.(rx.^2 + ry.^2); # radial location
scatter( rx[rr.<r], ry[rr.<r], markersize=1.5, markerstrokewidth=0)
scatter!(rx[rr.>r], ry[rr.>r], markersize=1.5, markerstrokewidth=0)
scatter!(aspect_ratio=:equal)
scatter!(legend=nothing)
we define a rectagle whose base is $[a,\,b]$, and whose height is larger than the maximum value of $f$ on $[a,\,b]$.
x = sort(rand(10))
y = rand(10)
p = poly_fit(x,y,5)
fp = Polynomial(p)
x = LinRange(0,1,1000)
y = fp.(x)
yshift = abs(minimum(y))
y = y .+ yshift
n = 10000
rx = rand(n)
ry = rand(n)*maximum(y)
fx = fp.(rx) .+ yshift
yhit = ry[ry .<= fx]
xhit = rx[ry .<= fx]
ymiss = ry[ry .> fx]
xmiss = rx[ry .> fx]
scatter( xhit, yhit, ms=1.5, markerstrokewidth=0)
scatter!(xmiss, ymiss, ms=1.5, markerstrokewidth=0)
plot!(x,y, lw=2, color="black")
plot!(legend=nothing)
plt=plot!(ylim=[0,1.1*maximum(y)])
display(plt)
Arectangle = 1*maximum(y)
I = Arectangle*length(xhit)/n
Ie = p[6]/6 + p[5]/5 + p[4]/4 + p[3]/3 + p[2]/2 + p[1] + yshift
println("I_mc = $I")
println("I_exact = $Ie")
I_mc = 23.735371322876222 I_exact = 23.61365454467787
Monte Carlo methods are particularly useful for integrating functions in high dimensions.
For example, in combustion systems we often deal with hundreds of chemical species.
$$\langle\omega_i(\vec{Y})\rangle = \int_{\vec{Y}}P(\vec{Y})\omega_i(\vec{Y})d\vec{Y}.$$
The integral is over the full multi-dimensional state space.
Such an integral cannot be computed using a 50 dimensional grid. Suppose we discretized each dimension with ten points. Then we would have $10^{50}$ total points to evaluate the function at, which is impossible.
Monte-Carlo methods reduce such calculations to manageable levels.