Issue 241 in mpmath: This integral will not work

[mpm...@googlecode.com]

Status: New
Owner: ----

New issue 241 by grant.wa...@gmail.com: This integral will not work
http://code.google.com/p/mpmath/issues/detail?id=241

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.



--
You received this message because this project is configured to send all
issue notifications to this address.
You may adjust your notification preferences at:
https://code.google.com/hosting/settings

[mpm...@googlecode.com]

Comment #1 on issue 241 by fredrik....@gmail.com: This integral will not
work
http://code.google.com/p/mpmath/issues/detail?id=241

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.

[mpm...@googlecode.com]

Comment #2 on issue 241 by grant.wa...@gmail.com: This integral will not
work
http://code.google.com/p/mpmath/issues/detail?id=241

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?

[mpm...@googlecode.com]

Comment #3 on issue 241 by fredrik....@gmail.com: This integral will not
work
http://code.google.com/p/mpmath/issues/detail?id=241

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).

[mpm...@googlecode.com]

Comment #4 on issue 241 by grant.wa...@gmail.com: This integral will not
work
http://code.google.com/p/mpmath/issues/detail?id=241

Ok, thank you! Everything makes sense now.