Wrong result when integrating (weighed) Chebyshev polynomials

[GoogleCodeExporter]

I'm logging this here as this seems to be a mpmath problem. The polynomial
'f' defined below is of order 2, and so it does not contain the 3rd
Chebyshev polynomial of the first kind: taking the scalar product
(integrating with the weight function 1/sqrt(1-x**2)) should return 0.
Sympy gets this right symbolically, but when integrating using quadts
the result depends on the given precision. For high levels of precision
quadts returns +inf, which is completely off. 

Some test code to illustrate this behaviour:

>>> from sympy.interactive import *
>>> from sympy import mpmath

>>> f = x**2 + x - 2

>>> weight = diff(asin(x), x)
>>> def integrand(n, func):
...     return weight * func * chebyshevt(n, x)

>>> mpmath.mp.dps = 100
>>> mpmath.quadts(lambda t : integrand(3, f).subs(x, t).evalf(), [-1, 1])
mpf('+inf')
>>> mpmath.mp.dps = 10
>>> mpmath.quadts(lambda t : integrand(3, f).subs(x, t).evalf(), [-1, 1])
mpf('-1.781695653964e-7')
>>> integrate(integrand(3, f), (x, -1, 1))
0

Original issue reported on code.google.com by jorn.baa...@gmail.com on 18 Oct 2009 at 11:42

[GoogleCodeExporter]

Well the following works, so it rather seems to be an issue with evalf 
returning +inf
close to the singularity at -1 even when it shouldn't:

>>> from mpmath import *
>>> n = 3
>>> f = lambda x: diff(asin, x) * (x**2+x-2)*chebyt(3,x)
>>> mp.dps = 100
>>> quad(f, [-1,1])
mpf('-2.059865464227744827411977669273198632525765165715921865963072495169589619
204642003714978954969971380444e-53')

It's inaccurate due to the singularity which is of type 1/sqrt(x). About half 
the
precision is lost due to the amplification of the uncertainty in the quadrature
nodes. The following simpler examples also demonstrate this:

>>> mp.dps = 10; 2-quad(lambda x: 1/sqrt(x), [0,1])
mpf('1.201842678711e-7')
>>> mp.dps = 20; 2-quad(lambda x: 1/sqrt(x), [0,1])
mpf('1.2219366736851552598941e-12')
>>> mp.dps = 30; 2-quad(lambda x: 1/sqrt(x), [0,1])
mpf('1.41170839398550092442887420002314e-17')
>>> mp.dps = 40; 2-quad(lambda x: 1/sqrt(x), [0,1])
mpf('1.402406641150261557057869824475216789991847e-22')

A simple workaround is to just do the integrations at twice the precision.

This problem would be possible to solve more elegantly in mpmath, either by 
always
computing nodes to double precision (which would be somewhat wasteful) or 
adding an
option to quad to handle this type of singularity better.

Original comment by fredrik....@gmail.com on 18 Oct 2009 at 12:22