Series expansion of erf() fails

[petschge]

I ran into the following problem with version 1.13.1 and 1.13.3. I did not other check other version.

Consider the following code:

import sympy

x = sympy.symbols("x", positive=True)

f1 = sympy.erf(x)
print(f1)
s1 = sympy.series(f1, x, x0=sympy.oo, n=3)
print(s1)

f2 = 2*sympy.erf(x)
print(f2)
s2 = sympy.series(f2, x, x0=sympy.oo, n=3)
print(s2)

For the first part it outputs

erf(x)
-(1/x + O(x**(-3), (x, oo)))*exp(-x**2)/sqrt(pi) + 1

which is reasonable. But for the second part it outputs

2*erf(x)
Traceback (most recent call last):
  File ".local/lib/python3.11/site-packages/sympy/core/expr.py", line 2994, in series
    s = self.subs(x, cdir/x).series(x, n=n, dir='+', cdir=1)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/expr.py", line 3023, in series
    s1 = self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/mul.py", line 1997, in _eval_nseries
    ords3 = [coeff_exp(term, x) for term in fac]
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/mul.py", line 1997, in <listcomp>
    ords3 = [coeff_exp(term, x) for term in fac]
             ^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/mul.py", line 1953, in coeff_exp
    lt = term.leadterm(x)
         ^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/expr.py", line 3527, in leadterm
    l = self.as_leading_term(x, logx=logx, cdir=cdir)
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/cache.py", line 72, in wrapper
    retval = cfunc(*args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/expr.py", line 3490, in as_leading_term
    obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/mul.py", line 2036, in _eval_as_leading_term
    return self.func(*[t.as_leading_term(x, logx=logx, cdir=cdir) for t in self.args])
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/mul.py", line 2036, in <listcomp>
    return self.func(*[t.as_leading_term(x, logx=logx, cdir=cdir) for t in self.args])
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/cache.py", line 72, in wrapper
    retval = cfunc(*args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/expr.py", line 3490, in as_leading_term
    obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/mul.py", line 2036, in _eval_as_leading_term
    return self.func(*[t.as_leading_term(x, logx=logx, cdir=cdir) for t in self.args])
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/mul.py", line 2036, in <listcomp>
    return self.func(*[t.as_leading_term(x, logx=logx, cdir=cdir) for t in self.args])
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/cache.py", line 72, in wrapper
    retval = cfunc(*args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/core/expr.py", line 3490, in as_leading_term
    obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".local/lib/python3.11/site-packages/sympy/functions/elementary/exponential.py", line 556, in _eval_as_leading_term
    raise PoleError("Cannot expand %s around 0" % (self))
sympy.core.function.PoleError: Cannot expand exp(-1/_k**2) around 0

instead of the expected

-(2/x + O(x**(-3), (x, oo)))*exp(-x**2)/sqrt(pi) + 2

[AlnesrAltair]

I believe that for your f1, the expression is an instance of the error function class and the series function used to define s1 ends up calling the series methods within the erf class - which are of course designed to handle the error function. However in the case of f2, this is an instance of the Mul class (ie multiplication of 2 and erf(x)) and we end up in the series methods as implemented in Mul, which it seems are not able to manage the error function.

More explicitly, we go into Mul._eval_nseries , which determines a Taylor series for each of the two multiplicative factors (2 and erf(x)) then attempts to zip the two series together. To do this it tries to inspect each term to determine the "leading term" and then uses these to multiply the two series term by term. In this case, we end up trying to find the leading term of exp(-1/x**2) which takes us through the local function coeff_exp within Mul._eval_nseries to exp.eval_as_leading_term which raises the PoleError.

I'm new to sympy so please let me know if I missed something.

I'm not sure what the fix should be. coeff_exp calls leadterm from within a try block, which is specifically looking to catch ValueErrors. Potentially this could be expanded to include PoleErrors as well - but I am not too familiar with the ways of sympy and am not sure what side-effects this would have.

[AlnesrAltair]

I've looked at my suggested fix again and for

f2 = 2*sympy.erf(x)
s2 = sympy.series(f2, x, x0=sympy.oo, n=3)
print(s2) 

I get the output:

2 + O(x**(-3), (x, 00))

rather than

-(2/x + O(x**(-3), (x, oo)))*exp(-x**2)/sqrt(pi) + 2

Note this means that application of series' toerf` would not be stable, in that:

f1 = sympy.erf(x)
s1 = sympy.series(f1, x, x0=sympy.oo, n=3)
print(s1)
s2 = sympy.series(s1, x, x0=sympy.oo, n=3)
print(s2)

would output:

-(1/x + O(x**(-3), (x, oo)))*exp(-x**2)/sqrt(pi) + 1
1 + O(x**(-3), (x, oo))

I think this is because the _eval_nseries method in the Mul class is set up to output power series (it explicitly constructs the output in terms of powers of x). What we get from the erf class is a more general asymptotic series. The output of series seems to be an Add class and is not aware of the type of series it is etc for this to be passed on to the next call of series.