ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Thu, 26 Aug 2010 09:48:47 +0200problems with symbolic integration and then numerical evaluatinghttps://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/can anyone explain this:
sage: integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1).n()
1.16508247725542e79
from approximation one know's that the legendre coefficients converge exponentially to zero and not to infinity!
and indeed with mpmath I get a better answer:
sage: import sage.libs.mpmath.all as mpmath
sage: mpmath.call(mpmath.quad,lambda x: mpmath.legendre(64,x)*mpmath.sin(pi/2*(x+1)),[-1,1])
-5.04684703543649e-25
Is there an overhead happening, when I numerically evaluate large rationals or something???
Thanks in advance,
maldun
Tue, 24 Aug 2010 19:08:05 +0200https://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/Answer by William Stein for <p>can anyone explain this:</p>
<pre><code>sage: integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1).n()
1.16508247725542e79
</code></pre>
<p>from approximation one know's that the legendre coefficients converge exponentially to zero and not to infinity!</p>
<p>and indeed with mpmath I get a better answer:</p>
<pre><code>sage: import sage.libs.mpmath.all as mpmath
sage: mpmath.call(mpmath.quad,lambda x: mpmath.legendre(64,x)*mpmath.sin(pi/2*(x+1)),[-1,1])
-5.04684703543649e-25
</code></pre>
<p>Is there an overhead happening, when I numerically evaluate large rationals or something???</p>
<p>Thanks in advance,
maldun</p>
https://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/?answer=11538#post-id-11538If you look at the output of
integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1)
you'll see that it is a massive expression with large rationals, powers of pi, etc. If `expr` is any Sage symbolic expression, then `expr.n(bits)` works by evaluating all the "leaves" of the expression to the given `bits` of precision, then doing the arithmetic to evaluate the expression. Thus large roundoff and cancellation can and sometimes will occur. Indeed, that is exactly what is happening here. You can get around this by increasing the precision sufficiently. If you want to be *sure* of the result, you can use interval arithmetic. Here's how:
sage: a = integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1)
sage: a.n(prec=300)
51516.6651093...60
sage: a.n(prec=500)
0
sage: a.n(prec=2000)
3.3204630...841332374315143e-97
sage: RealIntervalField(53)( a )
0.?e81
sage: RealIntervalField(200)( a )
0.?e37
sage: RealIntervalField(500)( a )
0.?e-54
sage: RealIntervalField(1000)( a )
3.32046...71175?e-97
For example, the lines
sage: RealIntervalField(500)( a )
0.?e-54
tell you that if evaluate the expression using 500 bits of precision, it's *definitely* 0.000000... (at least 53 zeros). Internally, this is done by evaluating with 500 bits of precision and rounding down and *also* evaluating with 500 bits and rounding up, and tracking the resulting "interval".
Tue, 24 Aug 2010 19:15:55 +0200https://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/?answer=11538#post-id-11538Comment by kcrisman for <p>If you look at the output of </p>
<pre><code>integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1)
</code></pre>
<p>you'll see that it is a massive expression with large rationals, powers of pi, etc. If <code>expr</code> is any Sage symbolic expression, then <code>expr.n(bits)</code> works by evaluating all the "leaves" of the expression to the given <code>bits</code> of precision, then doing the arithmetic to evaluate the expression. Thus large roundoff and cancellation can and sometimes will occur. Indeed, that is exactly what is happening here. You can get around this by increasing the precision sufficiently. If you want to be <em>sure</em> of the result, you can use interval arithmetic. Here's how:</p>
<pre><code>sage: a = integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1)
sage: a.n(prec=300)
51516.6651093...60
sage: a.n(prec=500)
0
sage: a.n(prec=2000)
3.3204630...841332374315143e-97
sage: RealIntervalField(53)( a )
0.?e81
sage: RealIntervalField(200)( a )
0.?e37
sage: RealIntervalField(500)( a )
0.?e-54
sage: RealIntervalField(1000)( a )
3.32046...71175?e-97
</code></pre>
<p>For example, the lines</p>
<pre><code>sage: RealIntervalField(500)( a )
0.?e-54
</code></pre>
<p>tell you that if evaluate the expression using 500 bits of precision, it's <em>definitely</em> 0.000000... (at least 53 zeros). Internally, this is done by evaluating with 500 bits of precision and rounding down and <em>also</em> evaluating with 500 bits and rounding up, and tracking the resulting "interval". </p>
https://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/?comment=22846#post-id-22846This is a great example that should make it into some doc/FAQ for those new to computer math software in general... which I will someday write...Thu, 26 Aug 2010 09:48:47 +0200https://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/?comment=22846#post-id-22846Comment by maldun for <p>If you look at the output of </p>
<pre><code>integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1)
</code></pre>
<p>you'll see that it is a massive expression with large rationals, powers of pi, etc. If <code>expr</code> is any Sage symbolic expression, then <code>expr.n(bits)</code> works by evaluating all the "leaves" of the expression to the given <code>bits</code> of precision, then doing the arithmetic to evaluate the expression. Thus large roundoff and cancellation can and sometimes will occur. Indeed, that is exactly what is happening here. You can get around this by increasing the precision sufficiently. If you want to be <em>sure</em> of the result, you can use interval arithmetic. Here's how:</p>
<pre><code>sage: a = integrate(legendre_P(64,x)*sin((1+x)*pi/2),x,-1,1)
sage: a.n(prec=300)
51516.6651093...60
sage: a.n(prec=500)
0
sage: a.n(prec=2000)
3.3204630...841332374315143e-97
sage: RealIntervalField(53)( a )
0.?e81
sage: RealIntervalField(200)( a )
0.?e37
sage: RealIntervalField(500)( a )
0.?e-54
sage: RealIntervalField(1000)( a )
3.32046...71175?e-97
</code></pre>
<p>For example, the lines</p>
<pre><code>sage: RealIntervalField(500)( a )
0.?e-54
</code></pre>
<p>tell you that if evaluate the expression using 500 bits of precision, it's <em>definitely</em> 0.000000... (at least 53 zeros). Internally, this is done by evaluating with 500 bits of precision and rounding down and <em>also</em> evaluating with 500 bits and rounding up, and tracking the resulting "interval". </p>
https://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/?comment=22896#post-id-22896that's an interesting behavior.
Thanks for your quick response!Tue, 24 Aug 2010 19:23:47 +0200https://ask.sagemath.org/question/7629/problems-with-symbolic-integration-and-then-numerical-evaluating/?comment=22896#post-id-22896