Ask Your Question
1

problem computing a numeric double integral

asked 2011-08-12 21:42:58 +0100

StevenPollack gravatar image

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

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
1

answered 2011-08-13 12:08:20 +0100

niles gravatar image

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)[0]
    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?

edit flag offensive delete link more

Comments

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.

StevenPollack gravatar imageStevenPollack ( 2011-08-13 12:21:44 +0100 )edit

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?

StevenPollack gravatar imageStevenPollack ( 2011-08-13 13:38:16 +0100 )edit

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

niles gravatar imageniles ( 2011-08-13 13:47:48 +0100 )edit

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?

StevenPollack gravatar imageStevenPollack ( 2011-08-13 17:54:33 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2011-08-12 21:42:58 +0100

Seen: 1,360 times

Last updated: Aug 13 '11