It looks like it's coming from Poly, since that is what's used to integrate polynomials. Poly ignores the precision on the coefficients and just uses 15.
The obvious workaround here is to not evalf until after you've integrated. Aaron Meurer On Sat, Nov 9, 2013 at 9:39 AM, Freddie Witherden <[email protected]> wrote: > Hi all, > > With SymPy 0.7.3 consider the following snippet: > > import sympy as sy > > p, q = sy.symbols('p q') > > p1 = sy.sqrt(2)/2 > p2 = 3*q/2 + sy.S(1)/2 > > p1f = p1.evalf(30) > p2f = p2.evalf(30) > > int1 = sy.integrate(sy.integrate(p1*p2, (q, -1, -p)), (p, -1, 1)) > int2 = sy.integrate(sy.integrate(p1f*p2f, (q, -1, -p)), (p, -1, 1)) > > print int1 > print int2 > > running this I find int1 is zero; as expected. However, int2 has a > residual of order 1e-16. But, given that both p1f and p2f are evaluated > to 30 decimal places (and I've verified both them and their product) > this implies that one of the integrate calls is truncating down to > floating point precision somewhere. > > Can anyone elucidate here? > > Regards, Freddie. > -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at http://groups.google.com/group/sympy. For more options, visit https://groups.google.com/groups/opt_out.
