Integral of real function has complex result

[eendebakpt]

A minimal example:

import sympy
from sympy import symbols, expand, simplify, integrate, S
t = symbols('t', positive=True)
pdf = ((1+t)/t**3)*sympy.exp(-1/t)
print(pdf)

E= integrate(pdf*t, (t, 0, 1e9))
print(E.evalf()) # seems fine

E= integrate(pdf*t, (t, 0, S.Infinity))
print(E.evalf()) # complex number!

has output

(t + 1)*exp(-1/t)/t**3
21.1460501720449
0.422784335098467 + 3.14159265358979*I

[oscarbenjamin]

Possibly related to #23337

[asmeurer]

It is coming from the computation of the FTC:

>>> t = symbols('t', positive=True)
>>> pdf = ((1+t)/t**3)*sympy.exp(-1/t)
>>> expr = integrate(pdf*t)
>>> expr
-Ei(-1/t) + exp(-1/t)
>>> limit(expr, t, 0)
0
>>> limit(expr, t, oo)
-EulerGamma + 1 + I*pi

I don't know if the limit is wrong or if the indefinite integral is simply inappropriate for naively computing the definite integral.

[eendebakpt]

@asmeurer @skirpichev I was unaware of the diofant package. It does a little better, so it might be worthwhile to look at changes there. But both sympy and diofant fail at the following

import diofant
from diofant import symbols, expand, simplify, integrate, S
t = symbols('t', positive=True)
pdf = ((1+t)/t**3)*diofant.exp(-1/t)

E= integrate(pdf*t, (t, 0, S.Infinity))
print(E.evalf()) # inf, whereas sympy fails

E= integrate(pdf*t, (t, 1, S.Infinity))
print(E.evalf()) # complex number!

[eendebakpt]

There is a problem with the limits:

import sympy
from sympy import Ei, symbols
from sympy.core.singleton import S

t = symbols('t', positive=True)

r=sympy.limit(Ei(1/t), t, S.Infinity, '+')
print(f'Ei(1/t) +: {r}')
r=sympy.limit(Ei(1/t), t, S.Infinity, '-')
print(f'Ei(1/t) -: {r}')

r=sympy.limit(Ei(-1/t), t, S.Infinity, '+')
print(f'Ei(-1/t) +: {r}')
r=sympy.limit(Ei(-1/t), t, S.Infinity, '-')
print(f'Ei(-1/t) -: {r}')

has output

Ei(1/t) +: EulerGamma
Ei(1/t) -: EulerGamma
Ei(-1/t) +: EulerGamma - I*pi
Ei(-1/t) -: EulerGamma - I*pi

[eendebakpt]

This term is already false:

t = symbols('t', real=True)
r=sympy.limit(Ei(t), t, 0 )

The series expansion of Ei(t) around zero should have a term log(t).

[anutosh491]

Not fully sure if there is an integration related error here but the series error here is quite closely related with another issue I plan to tackle during my GSoC project . So I tried approaching this and it seems fix for both issues is same. The following diff is what needs to be implemented here

diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index 1b0f58826f..e8f1f7bbc7 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -1212,10 +1212,15 @@ def _eval_rewrite_as_tractable(self, z, limitvar=None, **kwargs):
         return exp(z) * _eis(z)

     def _eval_as_leading_term(self, x, logx=None, cdir=0):
+        from sympy import re
         x0 = self.args[0].limit(x, 0)
+        arg = self.args[0].as_leading_term(x, cdir=cdir)
+        cdir = arg.dir(x, cdir)
         if x0.is_zero:
-            f = self._eval_rewrite_as_Si(*self.args)
-            return f._eval_as_leading_term(x, logx=logx, cdir=cdir)
+            if logx is not None:
+                return logx + S.EulerGamma - (
+                    S.ImaginaryUnit*pi if re(cdir).is_negative else S.Zero)
+            return log(-arg) + S.EulerGamma if re(cdir).is_negative else log(arg)
         return super()._eval_as_leading_term(x, logx=logx, cdir=cdir)

     def _eval_nseries(self, x, n, logx, cdir=0):
@@ -2016,7 +2021,9 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if arg0 is S.NaN:
             arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
         if arg0.is_zero:
-            return S.EulerGamma
+            if logx is not None:
+                return S.EulerGamma + logx
+            return log(x)
         elif arg0.is_finite:
             return self.func(arg0)
         else:
@@ -2251,6 +2258,8 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
         if arg0 is S.NaN:
             arg0 = arg.limit(x, 0, dir='-' if re(cdir).is_negative else '+')
         if arg0.is_zero:
+            if logx is not None:
+                return S.EulerGamma + logx
             return S.EulerGamma
         elif arg0.is_finite:
             return self.func(arg0)

The changes are the following

>>> integrate(((1 + x)/x**2)*exp(-1/x), (x, 0, oo))
1 - EulerGamma
>>> # After applying the diff we have
>>> integrate(((1 + x)/x**2)*exp(-1/x), (x, 0, oo))
oo

The diff would also solve the limit errors you've pointed out above @eendebakpt

The series expansion of Ei(t) around zero should have a term log(t).

>>> limit(Ei(t), t, 0 )
-oo

EDIT: I also see this above diff playing a part in one of the failing tests in my #23592 , so I'll fix the series based errors here once I figure the error I am having (related to failing limit in Ei function)

[eendebakpt]

@anutosh491 The def _eval_as_leading_term in Chi is indeed incorrect. Can you make a PR out of your patch?

[anutosh491]

@anutosh491 The def _eval_as_leading_term in Chi is indeed incorrect. Can you make a PR out of your patch?

Yeah will address this soon ! I feel this would fix the issue .