# Monte Carlo Double Integration Implementation

Am implementing a monte carlo integration routine to compute this double integral in eqn 0.3 of page 2 of this paper 'Mobius energy of knots and unknots', Annals of Mathematics, http://www.math.ucsb.edu/~zhenghwa/data/research/pub/Mobiusenergy-94.pdf in a case where both the integrals run over a circle.

$$E(\gamma)=\int\int \left( \frac{1}{|\gamma(v)-\gamma(u)|^2}-\frac{1}{D(\gamma(v),\gamma(u))^2} \right) |\dot{\gamma}(u)||\dot{\gamma}(v)|\text{d}u\text{d}v$$

For the sake of trivia, the first author of the paper in the above link, Michael H. Freedman is a Fields Medalist. $\gamma(.)$ denotes the parametric form of the circle and returns a 2d vector. $D(\gamma(u),\gamma(v))$ denotes the shortest arc-length on the circle between two points on it. The following is the code I wrote. If that works, then I'd be able to generalize it to other parametric curves. Can you help me do the below integration, properly? In the below code I've computed total length of the curve(circle) using the differential geometric notion of $\int_0^{2\pi}||\gamma'(t)||dt$.

Q1: In the case of a double-integral when u > v, R produces a NA or a negative result going by the rules of integration. In that case, how can a double integral be applied in a monte carlo routine, when both the integrals have the same limits of $0$ to $2\pi$ as there is always a chance of NA's occurring during the sample phase. Can you correct the implementation, so as to do the monte carlo double integration correctly?

Q2: Why is the monte carlo integration not producing a theoretical result of approximately 4, as suggested in the paper for the case of the circle?

library(cubature)
I.2d <- function(z) {
flag=0
FUN = function(x){
sqrt( (dxdt(x))^2 + (dydt(x))^2 + (dzdt(x))^2 )

}
dxdt <- function(x) {return(-1*sin(x))}
dydt <- function(x) {return(cos(x))}
dzdt <- function(x) {return(0)}
x = z[1]
y = z[2]
x1Coord = cos(x);y1Coord = sin(x);z1Coord = 0
x2Coord = cos(y);y2Coord = sin(y);z2Coord = 0
EuclideanDist = c(x1Coord,y1Coord,z1Coord) - c(x2Coord,y2Coord,z2Coord)
EuclideanDist = 1 / sqrt(sum(EuclideanDist^2 ))
totalLength = integrate(FUN, 0,2*pi)$value arcLength = integrate(FUN, x,y)$value
if(arcLength > (totalLength/2)){
GeodesicDistance = totalLength - arcLength
}
if(arcLength <= (totalLength/2)){
GeodesicDistance = arcLength
}
GeodesicDistance = 1 / (GeodesicDistance)^2
if(is.na(EuclideanDist-GeodesicDistance))
{
show("NA Found")
return(0)
}
if(!is.na(EuclideanDist-GeodesicDistance)){
return(EuclideanDist-GeodesicDistance)
}
}

adaptIntegrate(I.2d, c(0, 0), c(2*pi, 2*pi), maxEval=10000)

• To respond to your question: (1) No, you can't cross-post this exact question. If you change the emphasis significantly, you could post a significantly different version of this question on another Stack Exchange site. (2) Stack Overflow is a Stack Exchange site; all Stack Exchange sites are under the Stack Exchange umbrella. Jul 9 '15 at 21:08
• Thanks for clarifying on my question of whether I can post on multiple stackexchange / stackoverflow sites. Will have this question be just here in that case. Jul 10 '15 at 15:28

In this line:

EuclideanDist = 1 / sqrt(sum(EuclideanDist^2 ))


the sqrt should not be present: the integrand is one over the squared norm, not one over norm.

In this line:

arcLength = integrate(FUN, x,y)$value  this will give you a negative (incorrect) value for arc length whenever$x>y$, which confuses the rest of the code. The correct evaluation would be arcLength = abs(integrate(FUN, x,y)$value)


Finally, the reason for NAs is that the integration method will sometimes try to evaluate your functions at diagonal points with $x=y$, for which the integrand is undefined. One simple way to fix this is to replace the line

< y = z[2]
> y = z[2] + 1e-12


to break exact equality, or simply return zero when $x$ and $y$ are exactly equal.

The integrand FUN is $\sqrt{\cos^2x+\sin^2x} = 1$, so I'm not sure why you are integrating it numerically.

• Hey Kirill, I removed the sqrt and added 1e-12 as well. Upon integrating I get 9.18456 along with an error estimate of 7.704842. Why is the error bar so high-can it be improved? Also, it is still not close to 4 as mentioned in the paper, for the circle. What is still wrong-how can it be improved? Jul 10 '15 at 21:44
• @peow Did you make sure arcLength is positive by putting abs around it? I've got a modified working version of your code, so I've checked this answer is correct. Jul 10 '15 at 21:52
• Perfect! That bought it to 4.03442 with error bar of 0.2179383. Just add this comment into your above answer and the bounty is yours!. Will accept it. Jul 10 '15 at 21:54
• It says I may award bounty in 16 hours from now. Will give it then. Thanks again. Jul 10 '15 at 22:01