Added aseries support for error functions

[0sidharth]

References to other Issues or PRs

Continuing #17303

Brief description of what is fixed or changed

Added series support for error functions

In [1]: Si(x).series(x, oo)
Out[1]: 
    ⎛  6    1    ⎛1        ⎞⎞          ⎛24   2         ⎛1        ⎞⎞       
    ⎜- ── + ─ + O⎜──; x → ∞⎟⎟⋅sin(x)   ⎜── - ── + 1 + O⎜──; x → ∞⎟⎟⋅cos(x)
    ⎜   3   x    ⎜ 5       ⎟⎟          ⎜ 4    2        ⎜ 5       ⎟⎟       
π   ⎝  x         ⎝x        ⎠⎠          ⎝x    x         ⎝x        ⎠⎠       
─ - ──────────────────────────────── - ───────────────────────────────────
2                  x                                    x                 

Other comments

Release Notes

• functions
  • Added aseries expansion to error functions

[sympy-bot]

✅

Hi, I am the SymPy bot (v161). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• functions
  • Added aseries expansion to error functions (#20888 by @0sidharth)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.8.

Click here to see the pull request description that was parsed.

<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234" (see
https://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Continuing https://github.com/sympy/sympy/pull/17303

#### Brief description of what is fixed or changed
Added series support for error functions
```py
In [1]: Si(x).series(x, oo)
Out[1]: 
    ⎛  6    1    ⎛1        ⎞⎞          ⎛24   2         ⎛1        ⎞⎞       
    ⎜- ── + ─ + O⎜──; x → ∞⎟⎟⋅sin(x)   ⎜── - ── + 1 + O⎜──; x → ∞⎟⎟⋅cos(x)
    ⎜   3   x    ⎜ 5       ⎟⎟          ⎜ 4    2        ⎜ 5       ⎟⎟       
π   ⎝  x         ⎝x        ⎠⎠          ⎝x    x         ⎝x        ⎠⎠       
─ - ──────────────────────────────── - ───────────────────────────────────
2                  x                                    x                 

```
#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* functions
    * Added `aseries` expansion to error functions
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[0sidharth]

EDIT: So I think the issue is mainly with the symbolic form, and how sympy handles Order, and not related to the series code.

                                      2    
⎛   3      1     1    ⎛1        ⎞⎞  -x     
⎜- ──── + ──── - ─ + O⎜──; x → ∞⎟⎟⋅ℯ       
⎜     5      3   x    ⎜ 6       ⎟⎟         
⎝  4⋅x    2⋅x         ⎝x        ⎠⎠         

Is indeed O(1/x**6), if e^-x**2 is distributed inside. So the simplification is indeed correct, and I don't think a better expression can be obtained. But then the question is whether the series expansion should be implemented in its current mathematically correct way, or It can be simplified to return 1 + O(1/z**n, x) (Perhaps just S.One)

Current status:

In [1]: erf(x).series(x, oo)
Out[1]: 
     ⎛1        ⎞
1 + O⎜──; x → ∞⎟
     ⎜ 6       ⎟
     ⎝x        ⎠

Reposting from gitter -
When _eval_nseries calls _eval_aseries (L697 in function.py), adding a print statement says the asymptotic series being returned is

-(_x - _x**3/2 + 3*_x**5/4 + O(_x**6))*exp(-1/_x**2)/sqrt(pi) + 1

which is what I am expecting, but on returning this expression changes to

1 + O(x**(-6), (x, oo))

I think the answer sympy is finally returning is wrongly simplified and the original expression is better. This conversion is happening in _eval_nseries of Expr (L3024 in expr.py), where collect(s1, x) outputs the original expression but collect(s1, x) + o is outputting the changed expression. (where o is O(_x**6))

I have tried a few things to make sympy output the desired expressions instead, but none of them worked and I am stuck. I would appreciate any help. If this is indeed an undesirable simplification, is it a bug not in the implementation but on the current master?

[jksuom]

I think the answer sympy is finally returning is wrongly simplified and the original expression is better.

I agree. Asymptotic series are currently not handled in a satisfactory way. (The docstring of Expr.aserieseven says that 'This is equivalent to self.series(x, oo, n).' which is not true as asymptotic expansions can appear anywhere and expansions at oo can be convergent Laurent series.)

It seems to me that the series expansion code in expr.py needs to be improved but I have no detailed ideas about it. As a provisional fix I would try to apply this code

sympy/sympy/core/expr.py

Lines 3029 to 3032 in b6a1a3f

	try:
	return collect(s1, x) + o
	except NotImplementedError:
	return s1 + o

only if there s1 has no (embedded) Order term.

[jksuom]

Looks good, thanks! I have enabled auto-merge.

[0sidharth]

Looks like a test from ode took a long time on the CI -

sympy/solvers/ode/tests/test_single.py::test_factorable - Took 2742.024 seconds

However, I tried running this entire file again locally and it finished in the expected time -

(sympy-dev) (error_aseries) ~/sympy $ bin/test sympy/solvers/ode/tests/test_single.py
==================================================================== test process starts =====================================================================
executable:         /home/sidharth/anaconda3/envs/sympy-dev/bin/python  (3.8.5-final-0) [CPython]
architecture:       64-bit
cache:              yes
ground types:       gmpy 2.0.8
numpy:              1.19.4
random seed:        23845405
hash randomization: on (PYTHONHASHSEED=1622638637)

sympy/solvers/ode/tests/test_single.py[34] ..w...w......w.ww..www............                                                                             [OK]

_______________________________________________________________________ slowest tests ________________________________________________________________________
sympy/solvers/ode/tests/test_single.py::test_Liouville_ODE - Took 11.154 seconds
sympy/solvers/ode/tests/test_single.py::test_nth_order_linear_euler_eq_nonhomogeneous_variation_of_parameters - Took 11.595 seconds
sympy/solvers/ode/tests/test_single.py::test_nth_linear_constant_coeff_var_of_parameters - Took 11.849 seconds
sympy/solvers/ode/tests/test_single.py::test_1st_exact - Took 12.669 seconds
sympy/solvers/ode/tests/test_single.py::test_nth_order_linear_euler_eq_nonhomogeneous_undetermined_coefficients - Took 13.157 seconds
sympy/solvers/ode/tests/test_single.py::test_nth_linear_constant_coeff_homogeneous - Took 16.183 seconds
sympy/solvers/ode/tests/test_single.py::test_2nd_linear_bessel_equation - Took 20.801 seconds
sympy/solvers/ode/tests/test_single.py::test_nth_algebraic - Took 22.992 seconds
sympy/solvers/ode/tests/test_single.py::test_separable_reduced - Took 28.637 seconds
sympy/solvers/ode/tests/test_single.py::test_2nd_2F1_hypergeometric - Took 46.420 seconds
sympy/solvers/ode/tests/test_single.py::test_lie_group - Took 49.506 seconds
sympy/solvers/ode/tests/test_single.py::test_separable - Took 67.436 seconds
sympy/solvers/ode/tests/test_single.py::test_2nd_nonlinear_autonomous_conserved - Took 84.460 seconds
sympy/solvers/ode/tests/test_single.py::test_factorable - Took 139.178 seconds
================================================== tests finished: 26 passed, 8 skipped, in 584.51 seconds ===================================================

I don't think it is an issue with the PR because as I said I tried it again (although master has not been merged), so is this a known thing? What could be the possible cause of this?

[jksuom]

I got 106 sec on master and 103 sec on this branch, which I think is normal variation. I cannot see any connection with this PR.