5
$\begingroup$

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.

$\endgroup$
3
  • 2
    $\begingroup$ xi generates a new random number everytime it appears. So xi m[xi] actually generates two random numbers. But Monte-Carlo integrationswants to have the same random number both times. Also, you need the factor 1/n in both the numerator and denominator, no? $\endgroup$ Feb 2, 2023 at 18:20
  • 1
    $\begingroup$ I saw the multivaluedness in mcNumerator after I posted and 1/n was a typo. Fixed. Thanks. I will rewrite for get xi the same in xi m[xi] and test. $\endgroup$
    – JEM
    Feb 2, 2023 at 19:01
  • 1
    $\begingroup$ Just in case, surely you are aware of the Method -> "MonteCarlo" option on NIntegrate? $\endgroup$
    – kirma
    Feb 3, 2023 at 11:39

5 Answers 5

3
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ Very nice answer. Your explanation of why the original code converged to 6 was very helpful. $\endgroup$
    – JEM
    Feb 3, 2023 at 19:50
4
$\begingroup$

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];
al = (True / . CountsBy[r1, 0<#[[2]]<[[1]]m[#[[1]1] &]) 12 × 130 / 100000.
a2 = (True / . CountsBy[r2, 0<#[[2]] <m[#[ [1]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]]] &]]

enter image description here

$\endgroup$
2
$\begingroup$

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

$\endgroup$
5
  • $\begingroup$ Functional way to write this is something like With[{m = .06 #^2 - 3 # + 36 &}, Total[#2]/Total[#1] & @@ Transpose[{m[#], # m[#]} & /@ RandomReal[{0, 12}, 100000]]]. Maybe the use of anonymous functions (... &) a bit excessive on this, but using Do (especially in combination with AppendTo) is something most people working with Mathematica rather avoid, using Map, Table, Array and such instead. $\endgroup$
    – kirma
    Feb 3, 2023 at 12:14
  • $\begingroup$ ... or even With[{m = .06 #^2 - 3 # + 36 &}, Divide @@ Total[{# m[#], m[#]} & /@ RandomReal[{0, 12}, 100000]]]. $\endgroup$
    – kirma
    Feb 3, 2023 at 12:21
  • $\begingroup$ Or even Divide @@ Total[{#, 1} (.06 #^2 - 3 # + 36) & /@ RandomReal[{0, 12}, 100000]]. Okay, enough obfuscation is enough. :) $\endgroup$
    – kirma
    Feb 3, 2023 at 12:34
  • 1
    $\begingroup$ @ kirma I suspected the existence of an efficient functional approach to this but couldn't quite get there myself. Thank you. In your opinion, what is the best resource available that will take a traditional programmer and aid them in developing the kind of functional programming intuition you have displayed here? $\endgroup$
    – JEM
    Feb 3, 2023 at 19:42
  • $\begingroup$ It's an interesting question, I believe some of the old, and highest-voted answers here on Mathematica SE might dive more into it. I'm a bit biased on the subject since I learned to appreciate functional programming almost three decades ago on the introductory university CS course, and it stuck. I must admit that I was very much an imperative programming guy before that (but so was the hobbyist microcomputer world...). $\endgroup$
    – kirma
    Feb 4, 2023 at 3:36
2
$\begingroup$
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
$\endgroup$
-1
$\begingroup$

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
$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.