Monte Carlo Integration

Question

I need to find the center of mass of a bar of length 12 using basic Monte Carlo integration. The mass distribution $ m(x) = .06x^2-3x+36$ and the definition of the center of mass is

$x_{com} = \frac{\int_0^{12}xm(x)dx}{\int_0^{12}m(x)dx}$

Therefore, in the Wolfram Language(non Monte Carlo) we have

m[x_]:=.06x^2-3x+36
Integrate[x m[x],{x,0,12}]/Integrate[m[x],{x,0,12}]
 = 4.68966

The Monte Carlo Method employs the fact that $\lim_{n\rightarrow\infty}\frac{1}{n}\left(\sum_{i=1}^nf(x_i)\right) = $ The average of the function on the domain $[a,b]$ where $x_i$ is a uniformly distributed random real such that $a<x_i<b$.

$f_{avg} = \frac{1}{b-a}\int_a^bf(x)dx$ and so we can draw the following conclusion:

$\int_a^bf(x)dx\approx\frac{b-a}{n}\left(\sum_{i=1}^nf(x_i)\right) $

Therefore, I write

n=100000;
a=0;
b=12;
xi:=RandomReal[{a,b}]
mcNumerator = (1/n)Sum[xi m[xi],{i,1,n}];
mcDenominator = (1/n)Sum[m[xi],{i,1,n}];
MC = mcNUmerator/mcDenominator

The output of this little program converges to something near 6 for large $n$. What am I missing here?

Edit I see a potential problem in the mcNumerator where each of those $x_i$ are different.

Answer 1

The problem in your original code is that it generates different random values whenever you access the symbol xi (which is a function without arguments), so there are two different xi values in the xi * m[xi] - thus you actually computed two separated integrals for the numerator; more precisely, instead of computing

Integrate[x*m[x],{x,a,b}]

you computed

Integrate[x,{x,a,b}] * Integrate[m[x],{x,a,b}] / (b-a)

The numerator/denominator was then

Integrate[x,{x,0,12}] / 12

which is indeed 6.

The solution you proposed may be slightly inefficient - by appending one by one thousands (if not millions) of elements to the list there might be several reallocations and moving of data.

Still, your original version can be corrected by minor alterations:

m[x_]:=.06x^2-3x+36;
n=10000;
a=0;
b=12;
xi=RandomReal[{a,b}, n];
mcNumerator = ((b-a)/n)Total[xi* m[xi]];
mcDenominator = ((b-a)/n)Total[m[xi]];
MC = mcNumerator/mcDenominator

Of course, you may drop the (b-a)/n in front the two variables - I used them to check if each integral is computed correctly.

Answer 2

For something different:

In the following 2 rectangles that contain the regions of interest in numerator and denominator are randomly populated and fraction ‘under the respective curves’ is calculated and scaled:

m[x_] := 0.06  x^2 - 3  x + 36
r1 = ReIm@RandomComplex[{0, 12 + 130  I}, 100000];
r2 = ReIm@RandomComplex[{0, 12 + 36  I}, 100000];
a1 = (True /. CountsBy[r1, 0 < #[[2]] < #[[1]] m[#[[1]]] &])  12*130/
    100000.;
a2 = (True /. CountsBy[r2, 0 < #[[2]] < m[#[[1]]] &])  12*36/100000.;
a1/a2

This yielded 4.687

Visualising:

ListPlot[GroupBy[r1, 0<#[[2]]< #[[1]] m[#[[1]]] &]]
ListPlot[GroupBy[r2, 0<#[[2]]<m[#[[1]]] &]]

Correction 2024-04-06

@Mariusz Iwanuik correctly identified that my original code did not run. My original post contained multiple typographical errors (related to OCR errors pasting from image of code on cloud). I have corrected and confirmed as screenshot demonstrates.

Answer 3

I got this sorted out using a Do loop of all things:

m[x_] := .06 x^2 - 3 x + 36
data = {};
n=10000;
Do[{x = RandomReal[{0, 12}], AppendTo[data, {m[x], x m[x]}]},n]
t = Total[data]
t[[2]]/t[[1]]

This program converges to 4.68966

Answer 4

m[x_] := 0.06 x^2 - 3 x + 36
NIntegrate[x m[x], {x, 0, 12}]/NIntegrate[m[x], {x, 0, 12}]
n = 50000;
a = 0;
b = 12;
xi = Table[RandomVariate[UniformDistribution[{a, b}]], {i, 0, n}];
mcNumerator = Sum[xi[[i]] m[ xi[[i]] ], {i, 1, n}];
mcDenominator = Sum[m[xi[[i]]], {i, 1, n}];
MC = mcNumerator/mcDenominator

Answer 5

I apologize for not knowing how to enter Latex code to without triggering the "that's code; it must be indented" nagger.

m=Function[x, 0.06x^2-3 x+36]
x=RandomReal[12,10^4];
arm=Total[Map[# m[#]&,x]] ==> 
mass=Total[Map[m,x]]
arm/mass ==> 4.70536
ClearAll[x]
\[Integral]x m[x]\[DifferentialD]x ==> 18 x^2-x^3+(3 x^4)/200
\[Integral] m[x] \[DifferentialD]x ==> 36 x-(3 x^2)/2+x^3/50
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(12\)]\(x\ m[x] \[DifferentialD]x\)\) ==> 1175.04
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(12\)]\ \(m[x] \[DifferentialD]x\)\) --> 250.56
1175.04/250.56 ==> 4.68966