# 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))

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:

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) ) , (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

edit retag close merge delete

Sort by » oldest newest most voted

When I ran this code, I got an error about global variable v not being defined. I changed the last function to

def calcMobiusEnergy(gamma):
def f2(v):
def f1(u):
return mobiusEnergyIntegrand(gamma,u,v)
return numerical_integral(f1,0,2*pi)
return numerical_integral(f2,0,2*pi)


which, I think, lets sage work a little faster. Then I got warnings (not errors) about division by zero. I changed mobiusEnergyIntegrand to:

def mobiusEnergyIntegrand(gamma,u,v):
ed = euclideanDistance(gamma,u,v)
gal = geodesicArcLength(gamma,u,v)
if gal == 0:
return 0
else:
return speed(gamma)(u)*speed(gamma)(v)*(ed)^(2) - 1/(gal^(2))


and now the code has been running, with no warnings or errors (or output), for a few minutes. How long does it take for you to get your error message?

UPDATE:

I got bored, so I quit the calculation added %cython to the top of the (notebook) cell containing definitions, to see what would happen. The calculation finished surprisingly fast, but only after printing 't' many many times (once for each call to f1, which I determined by inserting extra print statements). For the circle, the answer I get is (0.0,0.0) -- can this possibly be right?

more

Thanks for the reply. I'll have to play around with your code. As far as I can tell, the circle cannot have mobius energy = 0, since you're summing non-negative quantities, and this would suggest that the integrand is constantly zero (something we know isn't true). P.S. how to do I do that %cython thing when not using the notebook? / Is there any documentation on that I could see? P.P.S. ed should be raised to the -2 inside your final return statement in mobiusEnergyIntegrand.

So, I just finished running calcMobiusEnergy on (cos(t), sin(t), 0), and found an energy equal to 400.8287... The number is a bit unsettling (it should probably be less that 6*pi + 4), but what's more troublesome is the speed in which the calculation is completed. Running sage from my computer (not the notebook) it took well over 10 minutes to calculate this. Any way to speed this up?

Great! Cython is the best way to speed up. Just telling sage to compile with cython gave substantial improvements for me. You can make further optimizations by telling it what type of objects the in/outputs of your functions will be. The docs are at http://docs.cython.org/src/quickstart/build.html

So I looked at the link and converted my code to http://sagenb.org/home/TeamSkynet/6/cells/1/__sagenb_servers_sage_notebook_sagenb_sagenb_home_TeamSkynet_6_code_sage112_spyx.html however when I run calcMobiusEnergy, I now get the error "Exception TypeError: 'unable to simplify to float approximation' in '_sagenb_servers_sage_notebook_sagenb_sagenb_home_TeamSkynet_6_code_sage\ 112_spyx_0.arcLength' ignored"... What gives?