ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sat, 13 Aug 2011 10:54:33 -0500problem computing a numeric double integralhttp://ask.sagemath.org/question/8270/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](http://ask.sagemath.org/question/68/double-integral) as a model for the function "calcMobiusEnergy".
Thanks,
StevenFri, 12 Aug 2011 14:42:58 -0500http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/Answer by niles for <p>So, I'd like to compute </p>
<p>$$
\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).</p>
<p>Here's my code:</p>
<pre><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) )
</code></pre>
<p>However, when I try and test this on something simple like</p>
<pre><code>sage: circle(t) = (cos(t), sin(t), 0)
sage: calcMobiusEnergy(circle)
</code></pre>
<p>I get</p>
<pre><code>---------------------------------------------------------------------------
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
</code></pre>
<p>Is my entire strategy just not right? Does anyone have any ideas on how to make this work? I used the response to <a href="http://ask.sagemath.org/question/68/double-integral">this old question</a> as a model for the function "calcMobiusEnergy". </p>
<p>Thanks,</p>
<p>Steven</p>
http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?answer=12573#post-id-12573When 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?
Sat, 13 Aug 2011 05:08:20 -0500http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?answer=12573#post-id-12573Comment by StevenPollack for <p>When I ran this code, I got an error about global variable <code>v</code> not being defined. I changed the last function to </p>
<pre><code>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)
</code></pre>
<p>which, I think, lets sage work a little faster. Then I got warnings (not errors) about division by zero. I changed <code>mobiusEnergyIntegrand</code> to:</p>
<pre><code>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))
</code></pre>
<p>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?</p>
<p>UPDATE:</p>
<p>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 <code>f1</code>, which I determined by inserting extra print statements). For the circle, the answer I get is <code>(0.0,0.0)</code> -- can this possibly be right?</p>
http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21380#post-id-21380So 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?Sat, 13 Aug 2011 10:54:33 -0500http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21380#post-id-21380Comment by StevenPollack for <p>When I ran this code, I got an error about global variable <code>v</code> not being defined. I changed the last function to </p>
<pre><code>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)
</code></pre>
<p>which, I think, lets sage work a little faster. Then I got warnings (not errors) about division by zero. I changed <code>mobiusEnergyIntegrand</code> to:</p>
<pre><code>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))
</code></pre>
<p>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?</p>
<p>UPDATE:</p>
<p>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 <code>f1</code>, which I determined by inserting extra print statements). For the circle, the answer I get is <code>(0.0,0.0)</code> -- can this possibly be right?</p>
http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21385#post-id-21385Thanks 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.Sat, 13 Aug 2011 05:21:44 -0500http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21385#post-id-21385Comment by StevenPollack for <p>When I ran this code, I got an error about global variable <code>v</code> not being defined. I changed the last function to </p>
<pre><code>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)
</code></pre>
<p>which, I think, lets sage work a little faster. Then I got warnings (not errors) about division by zero. I changed <code>mobiusEnergyIntegrand</code> to:</p>
<pre><code>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))
</code></pre>
<p>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?</p>
<p>UPDATE:</p>
<p>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 <code>f1</code>, which I determined by inserting extra print statements). For the circle, the answer I get is <code>(0.0,0.0)</code> -- can this possibly be right?</p>
http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21384#post-id-21384So, 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?Sat, 13 Aug 2011 06:38:16 -0500http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21384#post-id-21384Comment by niles for <p>When I ran this code, I got an error about global variable <code>v</code> not being defined. I changed the last function to </p>
<pre><code>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)
</code></pre>
<p>which, I think, lets sage work a little faster. Then I got warnings (not errors) about division by zero. I changed <code>mobiusEnergyIntegrand</code> to:</p>
<pre><code>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))
</code></pre>
<p>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?</p>
<p>UPDATE:</p>
<p>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 <code>f1</code>, which I determined by inserting extra print statements). For the circle, the answer I get is <code>(0.0,0.0)</code> -- can this possibly be right?</p>
http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21383#post-id-21383Great! 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.htmlSat, 13 Aug 2011 06:47:48 -0500http://ask.sagemath.org/question/8270/problem-computing-a-numeric-double-integral/?comment=21383#post-id-21383