problems with symbolic integration and then numerical evaluating

i like this post (click again to cancel)
1
i dont like this post (click again to cancel)

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

asked Aug 24 '10

maldun gravatar image maldun
81 3 4 9

updated Aug 24 '10

William Stein gravatar image William Stein
1414 5 20 46
http://wstein.org/
i like this answer (click again to cancel)
3
i dont like this answer (click again to cancel)

If 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".

link

posted Aug 24 '10

William Stein gravatar image William Stein
1414 5 20 46
http://wstein.org/
that's an interesting behavior. Thanks for your quick response! maldun (Aug 24 '10)
This 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... kcrisman (Aug 26 '10)

Your answer

Please start posting your answer anonymously - your answer will be saved within the current session and published after you log in or create a new account. Please try to give a substantial answer, for discussions, please use comments and please do remember to vote (after you log in)!
Login/Signup to Post

Question tools

Tags:

Stats:

Asked: Aug 24 '10

Seen: 553 times

Last updated: Aug 24 '10

powered by ASKBOT version 0.7.22
Copyright Sage, 2010. Some rights reserved under creative commons license.