Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

integrate can return a meijerg object that evaluates to nan #13681

Open
ethankward opened this issue Dec 6, 2017 · 4 comments
Open

integrate can return a meijerg object that evaluates to nan #13681

ethankward opened this issue Dec 6, 2017 · 4 comments

Comments

@ethankward
Copy link
Contributor

ethankward commented Dec 6, 2017

>>> from sympy import *
>>> x =  Symbol('x')
>>> r = integrate(log(1 + x**2) / (1 + x**2), [x, 0, oo])
>>> r
meijerg(((1/2, 1, 1), ()), ((1/2, 1), (0,)), 1)/2
>>> N(r)
nan

The meijerg result might be correct and just be evaluating incorrectly (I don't know), but the actual value of the integral is pi log(2).

@asmeurer
Copy link
Member

asmeurer commented Dec 6, 2017

It seems the value is correct. There are a lot of issues here.

If you replace 1 with z, you get

>>> pprint(hyperexpand(meijerg(((S(1)/2, 1, 1), ()), ((S(1)/2, 1), (0,)), z)/2))
     ⎛⎧            ⅈ⋅π             ⎞
     ⎜⎪acoth(√z) + ───             ⎟     ⎛⎧log(z - 1) - ⅈ⋅π  for │z│ > 1⎞
     ⎜⎪             2              ⎟   π⋅⎜⎨                             ⎟
     ⎜⎪───────────────  for │z│ > 1⎟     ⎝⎩  log(-z + 1)      otherwise ⎠
π⋅√z⋅⎜⎨       √z                   ⎟ + ──────────────────────────────────
     ⎜⎪                            ⎟                   2
     ⎜⎪   atanh(√z)                ⎟
     ⎜⎪   ─────────      otherwise ⎟
     ⎝⎩       √z                   ⎠

Now limit can't handle this directly (issue 1):

>>> limit(hyperexpand(meijerg(((S(1)/2, 1, 1), ()), ((S(1)/2, 1), (0,)), z)/2), z, 1)
Traceback (most recent call last):
  File "./sympy/series/limits.py", line 209, in doit
    r = gruntz(e, z, z0, dir)
  File "./sympy/series/gruntz.py", line 658, in gruntz
    r = limitinf(e0, z)
  File "./sympy/series/gruntz.py", line 428, in limitinf
    c0, e0 = mrv_leadterm(e, x)
  File "./sympy/series/gruntz.py", line 513, in mrv_leadterm
    series = calculate_series(f, w, logx=logw)
  File "./sympy/series/gruntz.py", line 466, in calculate_series
    for t in e.lseries(x, logx=logx):
  File "./sympy/core/expr.py", line 2696, in yield_lseries
    for si in s:
  File "./sympy/core/expr.py", line 2761, in _eval_lseries
    series = self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/add.py", line 386, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/add.py", line 386, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/expr.py", line 2848, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/mul.py", line 1606, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/mul.py", line 1606, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/expr.py", line 2848, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/add.py", line 386, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/add.py", line 386, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/expr.py", line 2848, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/function.py", line 669, in _eval_nseries
    raise PoleError("Cannot expand %s around 0" % (self))
sympy.core.function.PoleError: Cannot expand acoth(sqrt(_w + 1)) around 0

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "./sympy/series/limits.py", line 209, in doit
    r = gruntz(e, z, z0, dir)
  File "./sympy/series/gruntz.py", line 658, in gruntz
    r = limitinf(e0, z)
  File "./sympy/series/gruntz.py", line 428, in limitinf
    c0, e0 = mrv_leadterm(e, x)
  File "./sympy/series/gruntz.py", line 513, in mrv_leadterm
    series = calculate_series(f, w, logx=logw)
  File "./sympy/series/gruntz.py", line 466, in calculate_series
    for t in e.lseries(x, logx=logx):
  File "./sympy/core/expr.py", line 2696, in yield_lseries
    for si in s:
  File "./sympy/core/expr.py", line 2761, in _eval_lseries
    series = self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/mul.py", line 1606, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/mul.py", line 1606, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/expr.py", line 2848, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/add.py", line 386, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/add.py", line 386, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/expr.py", line 2848, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/function.py", line 669, in _eval_nseries
    raise PoleError("Cannot expand %s around 0" % (self))
sympy.core.function.PoleError: Cannot expand acoth(sqrt(_w + 1)) around 0

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "./sympy/series/limits.py", line 209, in doit
    r = gruntz(e, z, z0, dir)
  File "./sympy/series/gruntz.py", line 658, in gruntz
    r = limitinf(e0, z)
  File "./sympy/series/gruntz.py", line 428, in limitinf
    c0, e0 = mrv_leadterm(e, x)
  File "./sympy/series/gruntz.py", line 513, in mrv_leadterm
    series = calculate_series(f, w, logx=logw)
  File "./sympy/series/gruntz.py", line 466, in calculate_series
    for t in e.lseries(x, logx=logx):
  File "./sympy/core/expr.py", line 2696, in yield_lseries
    for si in s:
  File "./sympy/core/expr.py", line 2761, in _eval_lseries
    series = self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/mul.py", line 1606, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/mul.py", line 1606, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/expr.py", line 2848, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/add.py", line 386, in _eval_nseries
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/add.py", line 386, in <listcomp>
    terms = [t.nseries(x, n=n, logx=logx) for t in self.args]
  File "./sympy/core/expr.py", line 2848, in nseries
    return self._eval_nseries(x, n=n, logx=logx)
  File "./sympy/core/function.py", line 669, in _eval_nseries
    raise PoleError("Cannot expand %s around 0" % (self))
sympy.core.function.PoleError: Cannot expand acoth(sqrt(_w + 1)) around 0

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "./sympy/series/limits.py", line 209, in doit
    r = gruntz(e, z, z0, dir)
  File "./sympy/series/gruntz.py", line 658, in gruntz
    r = limitinf(e0, z)
  File "./sympy/series/gruntz.py", line 428, in limitinf
    c0, e0 = mrv_leadterm(e, x)
  File "./sympy/series/gruntz.py", line 487, in mrv_leadterm
    Omega, exps = mrv(e, x)
  File "./sympy/series/gruntz.py", line 309, in mrv
    "Don't know how to calculate the mrv of '%s'" % e)
NotImplementedError: Don't know how to calculate the mrv of '((acoth(sqrt(1 + 1/_p)) + I*pi/2)/sqrt(1 + 1/_p), True)'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "./sympy/series/limits.py", line 68, in limit
    return Limit(e, z, z0, dir).doit(deep=False)
  File "./sympy/series/limits.py", line 213, in doit
    r = heuristics(e, z, z0, dir)
  File "./sympy/series/limits.py", line 80, in heuristics
    l = limit(a, z, z0, dir)
  File "./sympy/series/limits.py", line 68, in limit
    return Limit(e, z, z0, dir).doit(deep=False)
  File "./sympy/series/limits.py", line 213, in doit
    r = heuristics(e, z, z0, dir)
  File "./sympy/series/limits.py", line 80, in heuristics
    l = limit(a, z, z0, dir)
  File "./sympy/series/limits.py", line 68, in limit
    return Limit(e, z, z0, dir).doit(deep=False)
  File "./sympy/series/limits.py", line 213, in doit
    r = heuristics(e, z, z0, dir)
  File "./sympy/series/limits.py", line 80, in heuristics
    l = limit(a, z, z0, dir)
  File "./sympy/series/limits.py", line 68, in limit
    return Limit(e, z, z0, dir).doit(deep=False)
  File "./sympy/series/limits.py", line 224, in doit
    raise NotImplementedError()
NotImplementedError

I'm guessing limits of Piecewise not working is already a known issue.

If you extract the z = 1 part, you get

>>> pprint(pi*atanh(sqrt(z)) + pi*log(1 - z)/2)
π⋅log(-z + 1)
───────────── + π⋅atanh(√z)
      2

But that still doesn't work with limit (issue 2):

>>> limit(pi*atanh(sqrt(z)) + pi*log(1 - z)/2, z, 1)
Limit(pi*log(-z + 1)/2 + pi*atanh(sqrt(z)), z, 1)

If you rewrite the atanh as log, though, it works, and gives the answer you said:

>>> limit(pi*atanh(sqrt(z)).rewrite(log) + pi*log(1 - z)/2, z, 1)
pi*log(2)

So there are a couple of issues with limit. And there are also two issues with meijerg:

  • hyperexpand doesn't work on the expression unless it is generalized (1 -> z)
  • evalf gives the wrong answer (nan)

@asmeurer
Copy link
Member

asmeurer commented Dec 6, 2017

Also want to say that in an ideal world even if SymPy can't compute the limit, hyperexpand would return Limit(pi*log(-z + 1)/2 + pi*atanh(sqrt(z)), z, 1).

@skirpichev
Copy link
Contributor

evalf issue sounds like mpmath's bug

@skirpichev
Copy link
Contributor

mpmath/mpmath#378

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants