How can I use Monte Carlo integration for non-rectangular regions?

[steven187]

Dear friends

I have been trying to write a code in Mathematica for solving a double integral with a non rectangular region through monte carlo integration, I have written a code that will solve for rectangular regions but i have no idea how to change it to work for non rectangular regions

for example let the function be fun[x_,y_] := 2x
what do i need to do to make this code to solve it for the non rectangular region {(x,y):y<=x<=4y,0<=y<=2}

fun[x_,y_] := 2x
a =
b =
c =
d =
Monte2D[n_]:=(b-a)*(d-c)/n*Sum[fun[Random[Real,{a,b}],Random[Real,{c,d}]],{n}]

nsamples=1000; nrepetitions=100;
mc=Table[Monte2D[nsamples],{nrepetitions}];

Print["For n=",nsamples," with ",nrepetitions," repetitions, \
mean=",Mean[mc]," st.dev.=",StandardDeviation[mc]


Please help

 

[juvenal]

You calculate the integral in the rectangular region, as you have done, but for the smaller domain, you need to check to see whether each (x,y) point is inside the smaller domain.

The way I would do it would be to keep two running totals: the sum in the rectangular region(=sum), and the number of (x,y) points generated that are inside the smaller domain(=npass).

The answer will then be: sum*npass/nsamples

I'm not sure if this can be done without revamping your code. But I'm not a expert on the intracacies of Mathematica.

 

[saltydog]

Thanks guys. It was very interesting and I learned a lot. This is what I came up with. I checked it against the integral of the indicated function over the specified domain. The integral is 40. The routine below returns 40.062.

fun[x_,y_]:=2 x;
MonteTest[fun_,nsamples_,nreps_,a_,b_,c_,d_]:=Module[{sum,tlist,xval,yval},
      tlist=Table[0,{nreps}];
      For[j=1,j<=nreps,j++,
        npass=0;
        sum=0;
        For[i=1,i<=nsamples,i++,
          xval=Random[Real,{a,b}];
          yval=Random[Real,{c,d}];
          If[yval<=xval && yval>=xval/4,
            sum+=fun[xval,yval];
            ];
          tlist[[j]]=sum*(b-a)*(d-c)/nsamples;
          ];
        ];
      Return[Mean[tlist]];
      ];
MonteTest[fun,1000,100,0,8,0,2]

Edit: Oh God. Should I not be doing this for Steve? You know Steve, I'm not sure this it right. Needs to be checked with other functions, other areas, perhaps optimized. It ran kind of slow. Don't use this for homework Ok?

Edit2: Ok, the tlist needs to go outside of the first For loop. Nevermind Steve.