Monte-Carlo Method

1. Feb 19, 2013

unscientific

1. The problem statement, all variables and given/known data

Hi everyone, I'm supposed to find the mass, centre of mass and moment of inertia of a toroid. The basic idea I have is:

a)Take a sample of N set of points (x,y,z)
b)Assign a value of 1 to each correct set of (x,y,z) that lie in the toroid (by seeing whether it satisfies the equation)
c)Sum up all the number of correct sets
d)Ratio R is defined as the sum/N
e)Volume of shape is R*(volume of region, which is usually a box)

Assuming p = 1 as in the question, I've managed to find the mass and centre of mass.

3. The attempt at a solution

The centre of mass Ixx is defined as mixx. I basically found the mass by simply finding the volume of the shape * 1 as ρ = constant = 1. For the centre of mass, I took the Ʃ(xcorrect * 1)/(number of correct sets) which simply gives the average x..

But for mixx i have this idea:
1. The number of sets of Ʃ(y^2 + z^2) you can take before you finish summing up the entire volume depends on how many partitions the volume is split up into.
2. Thus the number of sets taken is simply the number of partitions = Volume/number of correct sets

This gives a value for mixx to be around 82, miyy 140, mizz, 210

Code (Text):

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main()

{

int i=1;
double x;
double y;
double z;
int s;
double p;
double v;
double cmx1;
double cmx2;
double cmy1;
double cmy2;
double cmz1;
double cmz2;
double mixx1;
double mixx2;

srand48(time(NULL));

while (i<=100000)

{
x = 1 + 3 * (drand48());

y = 4 - 7 * (drand48());

z = 1 - 2 * (drand48());

if ( x*x + y*y + z*z - 6*sqrt(x*x + y*y) <= -8)
{s = 1;}
else
{s = 0;}

p = p + s;

v = (p/i)*42;

cmx1 = cmx1 + x*s;
cmy1 = cmy1 + y*s;
cmz1 = cmz1 + z*s;

cmx2 = cmx1/p;
cmy2 = cmy1/p;
cmz2 = cmz1/p;

mixx1 = mixx1 + (pow(y,2) + pow(z,2))*s;
mixx2 = (mixx1/p)*v;

i++;
printf("%d \t %.0f \t %f \t %f \t %f \t %f \t %f \t %f\n", s, p, v, cmx2, cmy2, cmz2, mixx1, mixx2);

}

}

Last edited: Feb 19, 2013
2. Feb 20, 2013

unscientific

bumpp?

3. Feb 22, 2013

unscientific

Nevermind, problem solved!