Approximate the value of an integral using simulation

Question

Calculate the following integral by generating $200$ of appropriate independent random variable.

$\int\int_{x_1 \geq 0, x_2\geq0, 1<x_1+x_2<6}x_1^2(1+x_2^2)^{-2}dx_1dx_2$

From what I currently know, the value of this integral is $\approx 42.1437$.

Here is what my code looks like so far in R:

set.seed(31)
n1=20000 #Number of vectors to generate
count=0
x1=rep(0,n1)
xt1=rep(0,n1)
xt2=rep(0,n1)
Y=rep(0,n1)
while (count<n1){
  x1=runif(1,0,1)*6
  x2=runif(1,0,1)*6
  if (x1+x2>1 & x1+x2<6)
  {count=count+1 
    xt1[count]=x1 
    xt2[count]=x2
  }

}
plot(xt1,xt2)
Y=xt1^2*(1+xt2^2)^(-2)
mean(Y)
##2.410251

I know something is wrong somewhere but I don't know where. If I look at the plot the points are in the correct region. I also decided to generate more vectors.

Answer 1

You previously asked this question and before you deleted it, one of my comments said "You need to multiply by the area of the quadrilateral bounded by $x_1\ge 0, x_2\ge 0, x_1+x_2 \gt 1, x_1+x_2\lt 6$".

The reason for this is that your simulation is giving an approximation of the average value of $x_1^2(1+x_2^2)^{-2}$ in that area, so to find what is in effect the volume under the function, you need to multiply by the area of integration.

That area is $\frac12 6^2-\frac121^2=17.5$ and $2.410251\times 17.5 \approx 42.18$ which is not far away from the exact value of $39\tan^{-1}(6)-\frac16\left(\log_e\left(\frac{37}{2}\right)+\pi +70\right)$

You can also simulate the area in your R code, since your random numbers are over an area of $6^2=36$, and your code will be quicker if you adjust it to avoid the loop. So perhaps something like

set.seed(31)
n1     <- 10^6 #Number of values to generate
x1     <- runif(n1, 0, 1) * 6
x2     <- runif(n1, 0, 1) * 6
inarea <- (x1+x2 > 1 & x1+x2 < 6)
mean(inarea) * 6^2
# 17.50741
xt1    <- x1[inarea]
xt2    <- x2[inarea]
Y      <- xt1^2 * (1+xt2^2)^(-2)
mean(Y)
# 2.400096
mean(Y) * mean(inarea) * 6^2
# 42.01947