problem computing a numeric double integral
So, I'd like to compute
∫2π0∫2π0(1|γ(u)−γ(v)|2−1D(γ(u),γ(v))2)|γ′(u)||γ′(v)|dudv where γ:[0,2π]→R3 is a curve, and D(γ(u),γ(v)) is the length of the shortest path (on the curve) between γ(u) and γ(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