GSoC 2019 Report Arighna Chakrabarty : Improving Series Expansions

The page summarizes the work that I have carried out during the course of the summer, and the work which is still left to be done in the future. The following wiki page contains the links to all my PRs in chronological order. People can kindly refer to my blog posts for more details.

About Me

My name is Arighna Chakrabarty, and I am a current final year student at the Indian Institute of Technology, Guwahati, and I am completing my degree in the field of Electronics and Communication Engineering.

Merged Pull Requests

The various Pull Requests enacted out in the course of the summer, and have been subsequently merged are --

1. Corrected the symbolic XFAIL tests in test_formal.py - There were some pre-existing symbolic functions for which the fps of the functions could not be calculated.

For example :-

x, n = symbols('x n')
f = x**n*sin(x**2)
assert fps(f,x).truncate(8) = x**2*x**n - x**6*x**n/6 + O(x**(n + 8), x)

The above test would fail previously. The fps of the above function would return a NotImplemented Error. This PR introduced the changes to correctly execute the symbolic fps tests. The polynomial function in the FormalPowerSeries class was tweaked to handle the symbolic cases.

2. sympy.fps return a formal power series object for polynomial functions - This was a minor issue which came up during my pre-GSoC days, and my then published PR had some errors and merge conflicts in it. So, I decided to come up with a new PR, whose working is described below --

Formerly, polynomial functions, when referenced through fps, did not return a Formal Power Series object, instead the function itself got returned.

>>> from sympy import Symbol, fps
>>> x = Symbol('x')
>>> p = fps(x ** 2)
>>> p
x**2
>>> type(p)
<class 'sympy.core.power.Pow'>
>>> p[0]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'Pow' object does not support indexing

I implemented the Coeff class according to the lines of Aaron's comment here. Further integration of class methods and functions might be required into the class. Coefficient sequence of the formal power series of polynomial functions contains a Coeff object now.

Now, A Formal Power Series object gets returned.

>>> from sympy import Symbol, fps
>>> x = Symbol('x')
>>> p = fps(x ** 2)
>>> p
FormalPowerSeries(x**2 + x + 1, x, 0, 1, (SeqFormula(Coeff(x**2 + x + 1, x, _k), (_k, 1, oo)), SeqFormula(x**_k, (_k, 0, oo)), 1))
>>> p.truncate()
1 + x + x**2 + O(x**6)
>>> p[0]
1

3. Implemented convolution operation of two formal power series - This PR introduces the convolution of two FormalPowerSeries objects. Presently, it takes in a order term, and prints all the terms in the convoluted resultant fps upto order. No FormalPowerSeries object is being returned as of now.

For example --

>>> f1, f2 = fps(sin(x)), fps(exp(x))
>>> f1.convolute(f2, x, order=6)
x + x**2 + x**3/3 - x**5/12 + O(x**6)

This PR contains the implementation of linear convolution. It only contains the code logic, thus not returning any FormalPowerSeries class object. The class logics have been integrated later.

4. Implementation of composition and inversion of formal power series - This PR introduces the omposition and inversion operations of two FormalPowerSeries objects. It takes in a order n term, and prints all the terms in the composed resultant fps upto order. No FormalPowerSeries object is being returned as of now.

For composition, there exists a closed form expression for the resultant power series, where the resultant coefficient sequence ak involves Bell polynomials, as follows --

For inversion, there exists a similar Bell-polynomial closed form expression. You can check the tests for this PR, for getting more ideas as to how fps(exp(sin(x))) might look like. This PR only contains the code logic, not any class wrappers as of now.

5. Introduction of FormalPowerSeries sub-classes, which is being returned by product, compose and inverse functions - This was the main class PR, one which binds all the operations PR. I created a FiniteFormalPowerSeries sub-class of FormalPowerSeries class, which contains some common code for the next three sub classes. All the three functions mentioned above have three classes of their own, namely FormalPowerSeriesProduct, FormalPowerSeriesCompose, and FormalPowerSeriesInverse. You can refer to the implementation in the above PR.

Let me show you the implementation of FormalPowerSeriesCompose.

class FormalPowerSeriesCompose(FiniteFormalPowerSeries):
    """Represents the composed formal power series of two functions.
    No computation is performed. Terms are calculated using a term by term logic,
    instead of a point by point logic.
    There are two differences between a `FormalPowerSeries` object and a
    `FormalPowerSeriesCompose` object. The first argument contains the outer
    function and the inner function involved in the omposition. Also, the
    coefficient sequence contains the generic sequence which is to be multiplied
    by a custom `bell_seq` finite sequence. The finite terms will then be added up to
    get the final terms.
    See Also
    ========
    sympy.series.formal.FormalPowerSeries
    sympy.series.formal.FiniteFormalPowerSeries
    """

    @property
    def function(self):
        """Function for the composed formal power series."""
        f, g, x = self.f, self.g, self.ffps.x
        return f.subs(x, g)

    def _eval_terms(self, n):
        """
        Returns the first `n` terms of the composed formal power series.
        Term by term logic is implemented here.
        The coefficient sequence of the `FormalPowerSeriesCompose` object is the generic sequence.
        It is multiplied by `bell_seq` to get a sequence, whose terms are added up to get
        the final terms for the polynomial.
        Examples
        ========
        >>> from sympy import fps, sin, exp, bell
        >>> from sympy.abc import x
        >>> f1 = fps(exp(x))
        >>> f2 = fps(sin(x))
        >>> fcomp = f1.compose(f2, x)
        >>> fcomp._eval_terms(6)
        -x**5/15 - x**4/8 + x**2/2 + x + 1
        >>> fcomp._eval_terms(8)
        x**7/90 - x**6/240 - x**5/15 - x**4/8 + x**2/2 + x + 1
        See Also
        ========
        sympy.series.formal.FormalPowerSeries.compose
        sympy.series.formal.FormalPowerSeries.coeff_bell
        """
        f, g = self.f, self.g

        ffps, gfps = self.ffps, self.gfps
        terms = [ffps.zero_coeff()]

        for i in range(1, n):
            bell_seq = gfps.coeff_bell(i)
            seq = (ffps.bell_coeff_seq * bell_seq)
            terms.append(Add(*(seq[:i])) / ffps.fact_seq[i-1] * ffps.xk.coeff(i))

        return Add(*terms)

Every subclass consists of an _eval_terms() custom functions, which actually contains the logic of the generating series. FOr more details about the wrapper classes and their method functions, visit the blog posts and the subsequent PR.

6. Added aseries function - This PR introduces the implementation of asymptotic expansion of series. This is equivalent to self.series(x, oo, n).

BEFORE

>>> from sympy import *
>>> x=symbols('x')
>>> e = sin(1/x + exp(-x)) - sin(1/x)
>>> e.series(x,oo)
O(x**(-6), (x, oo))

NOW

>>> from sympy import *
>>> x=symbols('x')
>>> e = sin(1/x + exp(-x)) - sin(1/x)
>>> e.aseries(x)
(1/(24*x**4) - 1/(2*x**2) + 1 + O(x**(-6), (x, oo)))*exp(-x)

Open Pull Requests

1. Added aseries for special functions - This PR includes the implementation of asymptotic and series expansions for these special functions --

    erf: gauss error function
    erfc: complementary error function
    erfi: imaginary error function
    Ei: exponential integral
    expint: generalised exponential integral
    li: logarithmic integral
    Li: offset logarithmic integral
    Si: sine integral
    Ci: cosine integral
    bessely
    besselj
    lerchphi function
    Riemann zeta function

There are some bugs to be solved, for example aseries of erf functions still gives you 1 + O(x**6), and not the defined asymptotic series.

2. Corrected some limit XFAIL tests - There are some tests, which currently fail under the coverage of the present limits module. I have noted down the issues in Issue #17380. Previously,

>>> from sympy import *
>>> x = symbols('x')
>>> limit(x*(((x + 1)**2 + 1)/(x**2 + 1) - 1), x, oo)
0

Now,

>>> from sympy import *
>>> x = symbols('x')
>>> limit(x*(((x + 1)**2 + 1)/(x**2 + 1) - 1), x, oo)
2

Future Work

A lot of things have been covered in the span of the summer. Firstly, I believe that I achieved my intended goals in the field of FormalPowerSeries. Operations have been implemented, and the class structure is also very friendly. But there are some work left to be done in the asymptotic expansions and in the limits module. I hope that I can continue these things post GSoC, and I also hope that this report will be helpful for anyone who wishes to contribute further in the field of Series Expansions. Following is a list of ideas, and leftovers that can be extended in this domain --

1. Implementation and extension of aseries functions for special functions. Almost all the work has been done in the PR mentioned above. Some bugs needs to be fixed, and then it will be ready to go.
2. Adding some more operations in FormalPowerSeries like reversion, or coefficient sequence of power of fps. The class structure is ready and flexible, so adding newer operations should be easy.
3. Correcting some more limit and convergent tests.

There are some points which were there in my proposal, but were left unimplemented. They are --

1. Improving ring_series module.
2. Unifying series expansion, and assigning one global class for all, for e.g. series, fps, ring_series should all come under one unified class hierarchy.

Conclusion

This summer has been a great learning experience. I plan to actively contribute to SymPy, specifically to this project. I am grateful to my mentor, Sartaj for reviewing my work, giving me valuable suggestions, and being readily available for discussions.