Monte Carlo Double Integration Implementation

Question

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)

Answer 1

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.