This integral will not work

[GoogleCodeExporter]

What steps will reproduce the problem?

>>> from mpmath import *
>>> mystring = 
"mpf('2')/(mpf('6')-mpf('2')*cos(y))**(mpf('1')/mpf('2'))*ellipk(mpf('2')/(mpf('
3')-cos(y)))"
>>> def func(y):
...     out = eval(mystring)
...     return out
... 
>>> mp.dps = 30
>>> myint = quad(func,[0,pi])
>>> print myint
+inf

What is the expected output? What do you see instead?

The output I expect should be similar to maple or scipy.
Maple answer: 6.34499346714920200531149676632
scipy answer: 6.3449934671492025
using Gauss-legendre: 6.34493993891195322762533970189

What version of the product are you using? On what operating system?
I'm using mpmath 0.17, on Mac OSX 10.6.8

Please provide any additional information below.

I tried some other things, such as trying different precisions (as high as 160, 
as low as 15), turning up maxdegree to 20.  I turned on the verbose option, and 
it outputs:
Integrating from 0 to 3.14159 (degree 1 of 6)
Integrating from 0 to 3.14159 (degree 2 of 6)
('Estimated error:', 'nan')
Integrating from 0 to 3.14159 (degree 3 of 6)

As an independent check, I compared the function "func" for various values of y 
in maple and as implemented above.  These agreed.

Original issue reported on code.google.com by grant.wa...@gmail.com on 24 Jun 2013 at 9:11

[GoogleCodeExporter]

2/(3-cos(y)) evaluates to 1 when y is close to 0, and then ellipk returns +inf

It should work if you wrap the integrand with, say, extraprec(100).

The reason this only happens with tanh-sinh quadrature is that it uses 
interpolation nodes much closer to the endpoints than the other algorithms.

Original comment by fredrik....@gmail.com on 26 Jun 2013 at 5:13

[GoogleCodeExporter]

Ah, okay I see, thank you!  Works for me now.  I apologize that this ended up 
not being an actual issue.

I guess I'm not really understanding something here.  If working at low 
precision, and evaluating at y close to 0, suppose we get something like 
1.0e500, which would count as infinity given a precision of 30, correct?

By using the extraprec decorator with sufficient n, func is able to resolve the 
integrand, but then this number gets passed out to the internals of the 
integrator, which are operating under the global precision of 30, no?  Or does 
the integrator tweak the precision as needed on its own?

Original comment by grant.wa...@gmail.com on 26 Jun 2013 at 6:58

[GoogleCodeExporter]

1.0e500 does not count as infinity. However, the numerical integration code in 
mpmath is based on absolute errors and not relative errors, so it might fail to 
converge if there are such large numbers.

The integrator does not tweak the precision except for adding a few guard bits 
against normal rounding error. If the precision is 100 bits, say, then quad() 
will set the precision to something like 110 bits and assume that this is 
sufficient for evaluating the integrand accurately to 100 bits.

If 200 bits are required to evaluate the integrand accurately to 100 bits, then 
user must increase the working precision appropriately inside the integrand.

The autoprec decorator can also be useful. In this case, it does not actually 
help since two consecutive precisions both happen to produce the same value 
(+inf).

Original comment by fredrik....@gmail.com on 26 Jun 2013 at 7:08

[GoogleCodeExporter]

Ok, thank you! Everything makes sense now.

Original comment by grant.wa...@gmail.com on 26 Jun 2013 at 7:21