Wrong result for digamma integral representation #26523
Comments
|
Looks like this comes from meijerg: In [11]: integrate((1 - t**z) / (1 - t), (t, 0, 1), meijerg=False)
Out[11]:
1
⌠
⎮ z
⎮ t - 1
⎮ ────── dt
⎮ t - 1
⌡
0
In [12]: integrate((1 - t**z) / (1 - t), (t, 0, 1), meijerg=True)
Out[12]:
⎛ 2⋅ⅈ⋅π ⎞ ⎛ 2⋅ⅈ⋅π ⎞
z⋅Φ⎝ℯ , 1, z + 1⎠⋅Γ(z + 1) Φ⎝ℯ , 1, z + 1⎠⋅Γ(z + 1)
- ────────────────────────────── - ──────────────────────────── + ∞
Γ(z + 2) Γ(z + 2)
In [13]: integrate((1 - t**z) / (1 - t), (t, 0, 1))
Out[13]: ∞ + ⅈ⋅π |
|
The result from meijerg=True looks like it would be correct if the In [58]: integrate((1 - t**z) / (1 - t), (t, 0, 1), meijerg=True).simplify()
Out[58]:
⎛ 2⋅ⅈ⋅π ⎞
- Φ⎝ℯ , 1, z + 1⎠ + ∞The |
|
@oscarbenjamin , I am not familiar with the Lerch transcendent, but evaluating the indefinite integral I get the following which results in |
That's coming from meijerg's indefinite integration: In [8]: integrate((1 - t**z) / (1 - t), t, meijerg=True).simplify()
Out[8]:
⎧ ⎛ z + 1 ⎛ 2⋅ⅈ⋅π ⎞ ⎞
⎪-⎝t ⋅Φ⎝t⋅ℯ , 1, z + 1⎠ + log(t - 1)⎠ for t > 1
⎨
⎪ ⎛ z + 1 ⎛ 2⋅ⅈ⋅π ⎞ ⎞
⎩-⎝t ⋅Φ⎝t⋅ℯ , 1, z + 1⎠ + log(1 - t)⎠ otherwise
In [9]: integrate((1 - t**z) / (1 - t), t).simplify()
Out[9]:
z + 1
- t ⋅Φ(t, 1, z + 1) - log(t - 1)
Maybe. Note that meijerg's definite integration is not really handled by computing the indefinite integral in the form shown. I think it would actually compute the definite integral using G-functions and then try to express the result in ordinary looking functions at the end. |
|
Ok. First things first, this is the first time I look into the sympy code for integration, so I might miss a lot of things (for example, I don't really know what meijerg's method is). import sympy as sym
t = sym.Symbol('t', nonnegative=True)
z = sym.Symbol('z')
f = (1 - t**z) / (1 - t)
I = sym.Integral((1 - t**z) / (1 - t), (t, 0, 1))
I.doit()and it seems that it is doing the following:
which evaluates to
but I don't think this is a problem since After that, it indeed tries to evaluate the primitive in the two limits in line
Then, for The problem then comes for 't = 1', which is evaluated at
This limit evaluates to infinite, which is the wrong part. The limit as Thus, in conclusion, it looks like the limit is the reason for the wrong value. |
What exactly is the expression whose limit is being computed? Can you reproduce the problem just using |
|
The expression is the primitive function: Indeed this can be reproduced with the following code: import sympy as sym
t = sym.Symbol('t', positive=True)
z = sym.Symbol('z')
F = (-t * t**z * z * sym.lerchphi(t, 1, z + 1) * sym.gamma(z + 1) / sym.gamma(z + 2)
- t * t**z * sym.lerchphi(t, 1, z + 1) * sym.gamma(z + 1) / sym.gamma(z + 2)
- sym.log(t - 1))
# An equivalent form
F = -t**(z+1)*sym.lerchphi(t, 1, z+1) - sym.log(t - 1)
L = sym.Limit(F, t, 1, '-')
L.doit()This return infinite, but as I justified, it should be digamma(z+1) + GammaEuler (- I*pi). |
|
So it looks like the bug is in limit... In [13]: L.args[0].subs(z, S(1)/2).subs(t, 0.9999).n()
Out[13]: 0.61365563825509 - 3.14159265358979⋅ⅈ
In [14]: L.args[0].subs(z, S(1)/2).limit(t, 1, '-')
Out[14]: ∞ |
|
Ok, I have looked into the limit implementation, and after going through the rabbit hole, I found the origin of the problem. import sympy as sym
a = sym.Symbol('a')
e = sym.Symbol('varepsilon', positive=True)
f = sym.lerchphi(1 - e, 1, a)
F = - (1 - e)**a * f - sym.log(e)
f._eval_as_leading_term(e, logx=None, cdir=1) # Out: lerchphi(1, 1, a)
F._eval_as_leading_term(e, logx=None, cdir=1) # Out: -lerchphi(1, 1, a) - log(varepsilon)The Lerch transcendent has no
If I add the following code to lerchphi (this is by no means a complete implementation), def _eval_as_leading_term(self, x, logx=None, cdir=0):
from sympy.functions.special.gamma_functions import digamma
z, s, a = self.args
z0, s0, a0 = (arg.subs(x, 0).cancel() for arg in self.args)
if s == 1 and z0 == 1 and not a.has(x):
lt = -log(-log(z))._eval_as_leading_term(x, logx=logx, cdir=cdir)
return lt - S.EulerGamma - digamma(a)
return super()._eval_as_leading_term(x, logx=logx, cdir=cdir)the previous code runs as expected: import sympy as sym
a = sym.Symbol('a')
e = sym.Symbol('varepsilon', positive=True)
f = sym.lerchphi(1 - e, 1, a)
F = - (1 - e)**a * f - sym.log(e)
f._eval_as_leading_term(e, logx=None, cdir=1) # Out: -log(varepsilon) - polygamma(0, a) - EulerGamma
F._eval_as_leading_term(e, logx=None, cdir=1) # Out: polygamma(0, a) + EulerGammaHowever, this still does not solve the original problem with the integration. The reason is that the function
But I guess that Taking the limit of the unsimplified primitive function, we arrive at the line
which returns
The previous expression fails this test. So the problem could also be solved by simply simplifying the coefficient, but again, I don't know if this might have a negative impact on performance. After that, it tries to apply the Gruntz algorithm, which still gives the wrong result: So I guess that trying to understand why the gruntz function returns the wrong result might be the best option. |
|
So looks like the root problem is a bug in |
|
Well, I wouldn't say that the root problem is in |
|
Okay, so two issue then:
To me the second problem is a bug but the first is just a missing feature. Of course it is still fine to fix the first but not the second if the second is not easy to fix but it does mean that the bug is still there. |
|
Ok, I was using the latest public version of sympy (1.12), where apparently there was a bug: whenever Working now with the master branch, the first problem I face is the following: import sympy as sym
t = sym.Symbol('t', positive=True)
z = sym.Symbol('z')
expr = -t**(z + 1)*sym.lerchphi(t, 1, z + 1) - sym.log(t - 1)
sym.gruntz(expr, t, 1, '-') # Out: ooThe problem is on these lines in
When doing the substitution, the output is By defining def _eval_is_finite(self):
z, s, a = self.args
if z == 1 and s == 1:
return Falsethe Gruntz algorithm no longer returns the wrong result; it fails with the following error: This should indeed be the expected result. However, the following code still gives a "wrong" result import sympy as sym
t = sym.Symbol('t', positive=True)
z = sym.Symbol('z')
expr = -t**(z + 1)*sym.lerchphi(t, 1, z + 1) - sym.log(t - 1)
L = sym.Limit(expr, t, 1, '-')
res = L.doit() # Out: -lerchphi(1, 1, z + 1) + ooOf course, one could argue that this is not strictly a wrong result since In summary:
|
The problem is an expression like this: In [1]: x = symbols('x')
In [2]: print(x.is_finite)
NoneI would say that it is not robust to depend on
I agree with all these points. Adding the missing methods is a useful thing to do regardless of anything else. |
Ok, as I said, I'm not sure if it is worth it to fix this problem since 99% of the time it should work fine and it would probably break a lot of limits. However, here is an option: We could add a condition like Then, regarding heuristics, I found the problem. It is essentially the same as the rest; when calling heuristics to a function, it simply computes the limits of each argument and then evaluates the function at the argument, which doesn't work if the function is infinite. Thus, I would add diff --git a/sympy/series/limits.py b/sympy/series/limits.py
index b3976e5512..e9553e58fc 100644
--- a/sympy/series/limits.py
+++ b/sympy/series/limits.py
@@ -115,7 +115,7 @@ def heuristics(e, z, z0, dir):
l = limit(e3, z, z0, dir)
rv = l * Mul(*r2)
- if rv is S.NaN:
+ if rv is S.NaN or e.is_Function and rv.is_finite is False:
try:
from sympy.simplify.ratsimp import ratsimp
rat_e = ratsimp(e)I don't know whether it is better to use Also, I realised there is another potential bug in heuristics, specifically in this part:
The idea is to simplify the expression and try the limit again, with an implicit assumption that if heuristics is called again, then e will be already simplified, and thus the check
This is problematic since, in some cases, import sympy as sym
from sympy.simplify import ratsimp, powsimp
t = sym.Symbol('t', positive=True)
z = sym.Symbol('z')
expr = t**(z+1)
rat_e = ratsimp(expr) # t*t**z
pow_rat_e = powsimp(rat_e) # t**(z + 1) == exprSo, even though I didn't find this problem for any real example (since the functions need to be more complicated for the limits to fail in a particular way to see this), this has the potential to fall into an infinite loop. With this, the limit of However, everything here still depends on the original expression being simplified. The original integral now evaluates to
This rewrites the expression
it returns
This would indeed fix the original problem. Of course, a full |
|
I'm finding it difficult to follow all of this so I'm not sure exactly what to suggest. Probably there are various things that could be fixed. I would say let's start with the most obvious fixes and get them in first. In general More generally it would be better if the limits etc code used some sort of canonical representation rather than lots of heuristic checks and simplification. |
I was trying to evaluate the following integral using sympy 1.12 (see complete list of packages below)
package listThe integration returns
oo + I*pi. The integral is a standard representation for the digamma function.sympy is able to do the integral correctly for the different concrete values of
zI have tried, but not for the general value ofz.The integral also gives the incorrect result even if I add the constrain that
zis positive.Also, sympy doesn't seem to be able to recognise the integral representation of the polygamma functions, though in that case it simply returns the unevaluated integral.
The text was updated successfully, but these errors were encountered: