我正在实施一个蒙特卡洛积分例程来计算本文第 2 页的 eqn 0.3 中的双重积分，“结和未结的莫比乌斯能量”，数学年鉴，http://www.math.ucsb.edu/~zhenghwa/data /research/pub/Mobiusenergy-94.pdf在两个积分都在一个圆上运行的情况下。

$$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$$

为方便起见，上述链接中论文的第一作者 Michael H. Freedman 是菲尔兹奖得主。$\gamma(.)$表示圆的参数形式并返回一个二维向量。$D(\gamma(u),\gamma(v))$表示圆上两点之间的最短弧长。以下是我写的代码。如果可行，那么我可以将其推广到其他参数曲线。你能帮我正确地做下面的整合吗？在下面的代码中，我使用微分几何概念计算了曲线（圆）的总长度$\int_0^{2\pi}||\gamma'(t)||dt$.

Q1：在双积分的情况下，当 u > v 时，根据积分规则，R 产生 NA 或负结果。在那种情况下，当两个积分具有相同的极限时，如何在蒙特卡罗例程中应用双积分$0$到$2\pi$因为在采样阶段总是有可能发生 NA。您能否更正实施，以便正确进行蒙特卡洛双积分？

Q2：为什么蒙特卡洛积分没有产生大约 4 的理论结果，正如论文中针对圆的情况所建议的那样？

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)