Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

problem computing a numeric double integral

So, I'd like to compute

$$ \int_{0}^{2\pi} \int_{0}^{2\pi} \left( \frac{1}{\left| \gamma(u) - \gamma(v)\right|^2} - \frac{1}{D(\gamma(u), \gamma(v))^2}\right) \left| \gamma'(u) \right| \left| \gamma'(v)\right| du dv $$ where $\gamma : [0,2\pi] \to \mathbb{R}^3$ is a curve, and $D(\gamma(u), \gamma(v))$ is the length of the shortest path (on the curve) between $\gamma(u)$ and $\gamma(v)$ (the geodesic arclength).

Here's my code:

def speed(gamma):
"""returns norm{gamma'}, the speed of the curve gamma(t)"""
return gamma.diff(t).norm().function(t)

def arcLength(gamma,a,b):
    """Returns the arc-length of a curve, gamma = f(t), from t=a to t=b """
    return abs(numerical_integral(speed(gamma), a, b)[0])

def geodesicArcLength(gamma,a,b):
    """finds the smallest distance between points t=a and t=b on a curve, gamma(t), t \in [0,2*pi]"""
    totalLength = arcLength(gamma,0,2*pi)
    if a <= b:
        lengthOfPath1 = arcLength(gamma,a,b)
    else:
        lengthOfPath1 = arcLength(gamma,b,a)

    if lengthOfPath1 <= (totalLength/2.0):
        return lengthOfPath1
    else:
        return totalLength - lengthOfPath1

def euclideanDistance(gamma,a,b):
"""finds the euclidean distance on gamma(t), a curve, between gamma(a) and gamma(b)"""
return norm(gamma(a) - gamma(b))

def mobiusEnergyIntegrand(gamma,u,v):
return speed(gamma)(u)*speed(gamma)(v)/(euclideanDistance(gamma,u,v)^(-2) - geodesicArcLength(gamma,u,v)^(-2))

def calcMobiusEnergy(gamma):
return numerical_integral( lambda v: numerical_integral( lambda u: mobiusEnergyIntegrand(gamma,u,v) , (u,0,2*pi) )[0] , (v,0,2*pi) )

However, when I try and test this on something simple like

sage: circle(t) = (cos(t), sin(t), 0)
sage: calcMobiusEnergy(circle)

I get

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)

/home/steven/sage-4.7/local/lib/python2.6/site-packages/sage/all_cmdline.pyc in <module>()

/home/steven/sage-4.7/local/lib/python2.6/site-packages/sage/all_cmdline.pyc in calcMobiusEnergy(gamma)

/home/steven/sage-4.7/local/lib/python2.6/site-packages/sage/gsl/integration.so in sage.gsl.integration.numerical_integral (sage/gsl/integration.c:2432)()

/home/steven/sage-4.7/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.__float__ (sage/symbolic/expression.cpp:5543)()

TypeError: unable to simplify to float approximation

Is my entire strategy just not right? Does anyone have any ideas on how to make this work? I used the response to this old question as a model for the function "calcMobiusEnergy".

Thanks,

Steven