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