UD Sequences and formal power series prototype

Alexey U. Gudchenko edited this page Feb 27, 2014 · 14 revisions

UD - Sequences and formal power series - prototype

 

Note: these are just the ideas of one person. It may not necessarily coincide with the ideas of others in the community.

Note: For GSoC, it doesn't have to be this way. On the contrary, we welcome the participants to describe their own ideas.

 

In this page we will introduce the working prototype which is in the sequences branch. Mostly all of the examples are given from the docstrings of the appropriated classes in this branch.

>>> from sympy import S, oo, Symbol, factorial
>>> from sympy.abc import x, y, k
>>> from sympy.sequences import Sequence, abstract_sequences, SeqFormula
>>> from sympy.sequences.expr import SeqMulEW
>>> from sympy.series.power import PowerSeries

Introduction

This approach is based on that we introduce the Sequence expression. And then use this sequence to define the coefficients for the constructing of formal power series. E.g. take periodical sequence [1, 2, 3, 1, 2, 3, ...]:

>>> seq = Sequence((0, oo), periodical=(1, 2, 3))
>>> ps = PowerSeries(x, sequence=seq)
>>> ps
1 + 2*x + 3*x**2 + x**3 + 2*x**4 + 3*x**5 + x**6 + 2*x**7 + 3*x**8 + ...

>>> _[11:]
3*x**11 + x**12 + 2*x**13 + 3*x**14 + x**15 + ...

The reasons for it are:

  • we work only with the coefficients of the series, and do not parse the expression, and represent series as the sum of sub-expressions and powers of formal variable only when print series.
  • some methods with formal power series are similar to the methods of sequences.
  • we can introduce various kind of sequences. So the implementation for specific kind of the sequences will be easer (and sometimes faster) rather then in common case.

Sequences

The main content of the Sequences is its coefficients:

>>> from sympy.sequences import Sequence
>>> from sympy.printing.pretty.pretty import pprint

>>> a = Sequence('a')
>>> pprint(a)
[a0, a1, a2, a3, a4, a5, a6, ...]

Another component is an interval, in which the index go through:

>>> a.interval
[0, oo)

The interval is helpful to do not use the elements which are obviously zero.

>>> a = Sequence((3, 6), 'a')
>>> pprint(a)
[0, ..., a3, a4, a5, a6]

>>> b = Sequence((5, 9), 'b')
>>> pprint(b)
[0, ..., b5, b6, b7, b8, b9]

>>> (a + b).interval    # element-wise summation
[3, 9]

>>> (a * b).interval    # Couchy product
[8, 15]

>>> SeqMulEW(a, b).interval         # element-wise multiplication
[5, 6]

Specific types of sequences

Specific types of sequences are used, as it is mentioned above, to implement specified algorithms (e.g. for periodical sequences). Or for the specific representation (e.g. abstract sequence). So there are implemented appropriated SumPy expressions.

It is possible to use specific constructor for them. Another way is to use the universal constructor with appropriated parameters.

Symbolic (abstract) sequence:

>>> a = Sequence('a')
>>> a
{a}

>>> pprint(a)
 [a0, a1, a2, a3, a4, a5, a6, ...]

>>> a[2]
a[2]

Periodical sequence:

>>> seq = Sequence(periodical=(1, 2, 3))
>>> seq
SeqPer([0, oo), (1, 2, 3))

>>> pprint(seq)
[1, 2, 3, 1, 2, 3, 1, ...]

Another way to construct is to use the SeqPer class directly:

>>> from sympy.sequences import SeqPer
>>> seq = SeqPer((4, oo), (1, 2, 3))
>>> pprint(seq)
[0, ..., 1, 2, 3, 1, 2, 3, 1, ...]

The same is for other kinds.

Finite sequence:

>>> seq = Sequence((0, 2), finitlist=(1, 2, 3))
>>> seq
SeqList([0, 2], (1, 2, 3))

>>> seq.is_infinite
False

By formula:

>>> seq = Sequence((3, oo), formula=(k, S(1)/k))
>>> seq
SeqFormula([3, oo), k, 1/k)

>>> pprint(seq)
[0, ..., 1/3, 1/4, 1/5, 1/6, 1/7, 1/8, 1/9, ...]

By function:

>>> f = lambda k: S(1)/k
>>> seq = Sequence((1, oo), function = f)
>>> seq
SeqFunc([1, oo), <function <lambda> at ...


>>> pprint(seq)
[0, 1, 1/2, 1/3, 1/4, 1/5, ...]

By recursion:

>>> from sympy import Function
>>> from sympy.abc import n
>>> y = Function("y")
>>> eq = (n-1)*y(n+2) - (n**2+3*n-2)*y(n+1) + 2*n*(n+1)*y(n)

>>> seq = Sequence((2, oo), recurr=(eq, y(n), { y(0):0, y(1):3 }) )
>>> pprint(seq)
[0, ..., 6, 6, -24, -264, -1968, ...]

>>> seq.formula
3*2**n - 3*n!

Operation with sequences

The summation is obviously element-wise:

>>> a = Sequence(periodical = (1, 0))
>>> b = Sequence((4, oo), periodical = (0, 2))

>>> a + b
SeqPer([0, oo), (1, 0)) + SeqPer([4, oo), (0, 2))

>>> pprint(a + b)
[1, 0, 1, 0, 1, 2, 1, ...]

We can use also the multiplication of the sequence by the scalar expression:

>>> s = Sequence(periodical = (1, 0))
>>> y = Symbol('y')
>>> 2*s
2*SeqPer([0, oo), (1, 0))

>>> pprint(3*s)
[3, 0, 3, 0, 3, 0, 3, ...]

>>> (3*s).coefficient
3

>>> (3*s*y*2).coefficient
6*y

So we can use the linear combination of them:

>>> pprint(3*s - b)
[3, 0, 3, 0, 3, -2, 3, ...]

The multiplication of two sequences can be done in two ways. One way is an element-wise multiplication:

>>> a = Sequence('a')
>>> b = Sequence(periodical=(1, 0, 2))
>>> c = SeqMulEW(a, b)
>>> c
SeqMulEW({a}, SeqPer([0, oo), (1, 0, 2)))

>>> pprint(c)
[a0, 0, 2*a2, a3, 0, 2*a5, a6, ...]

>>> c[11]
2*a[11]

factorization

>>> a = Sequence('a')
>>> b = Sequence(formula=(k, S(1)/factorial(k)))
>>> c = SeqMulEW(a, b)
>>> pprint(c)
             a2  a3  a4   a5   a6      
    [a0, a1, --, --, --, ---, ---, ...]
             2   6   24  120  720   

Another way is the Cauchy product or discrete convolution of sequences.

c_n=\sum_{k=0}^n a_k b_{n-k}.

This way is used as the main method for the binary operation '*':

>>> a, b = abstract_sequences('a,b')
>>> c = a*b
>>> c
{a}*{b}

>>> c[0]
a[0]*b[0]

>>> c[1]
a[0]*b[1] + a[1]*b[0]

>>> c[5]
a[0]*b[5] + a[1]*b[4] + a[2]*b[3] + a[3]*b[2] + a[4]*b[1] + a[5]*b[0]

We can define the power of sequence (exponentiation) to an natural number n as the n times of Cauchy product by itself.

>>> a = Sequence('a')
>>> a**3
{a}**3

>>> c = a**10

>>> c[0]
a[0]**10

>>> c[1]
10*a[0]**9*a[1]

>>> c[2]
10*a[0]**9*a[2] + 45*a[0]**8*a[1]**2

Note, that we can use the exponentiation not only with natural exponent, but also for with rational and negative exponent.

For this definition we can consider the result to be valid for the Formal Power Series algorithm.

Shifting of sequence:

>>> a = Sequence('a')
>>> pprint(a.shiftleft(2))
[a2, a3, a4, a5, a6, a7, a8, ...]

>>> pprint(a.shiftright(2))
[0, ..., a0, a1, a2, a3, a4, a5, a6, ...]

>>> pprint(a.shiftright(2).shiftleft(2))
[a0, a1, a2, a3, a4, a5, a6, ...]

It is used later for the Formal Power Series, to align the sequence to be with non-zero first element.

Formal power series

The natural construction of the formal power series is based on some sequence (coefficients) and the formal variable.

>>> a = PowerSeries(x, sequence=Sequence((1, oo), formula=(k, S(1)/k)))
>>> a
x + x**2/2 + x**3/3 + x**4/4 + x**5/5 + x**6/6 + x**7/7 + x**8/8 + ...

>>> a[9]
x**9/9

>>> a.coeff(9)
1/9

>>> a = PowerSeries(x, sequence=Sequence(periodical=(1, 0)))
>>> a
1 + x**2 + x**4 + x**6 + x**8 + ...

or shortly implicitly

>>> a = PowerSeries(x, periodical=(1, 0))
>>> a
1 + x**2 + x**4 + x**6 + x**8 + ...

>>> a = PowerSeries(x, 'a')
>>> a
a[0] + x*a[1] + x**2*a[2] + x**3*a[3] + x**4*a[4] + ...

Summation of the series is obviously element-wise:

>>> a = PowerSeries(x, periodical=(1, 0))
>>> b = PowerSeries(x, periodical=(0, 2))

>>> a + b
1 + 2*x + x**2 + 2*x**3 + x**4 + 2*x**5 + x**6 + 2*x**7 + x**8 + ...

The formula of multiplication is similar Cauchy product of sequenses:

>>> a = PowerSeries(x, 'a')
>>> b = PowerSeries(x, 'b')

>>> c = a*b

>>> c[0]
a[0]*b[0]

>>> c[1]
x*(a[0]*b[1] + a[1]*b[0])

>>> c[2]
x**2*(a[0]*b[2] + a[1]*b[1] + a[2]*b[0])

>>> c.coeff(2)
a[0]*b[2] + a[1]*b[1] + a[2]*b[0]

For the multiplication algorithm see:

Donald E. "Knuth Art of Computer Programming, Volume 2: Seminumerical Algorithms",
3rd ed., sec 4.7 "Manipulation of power series", p 525.

The exponentiation of the series is more interesting. See

Donald E. "Knuth Art of Computer Programming, Volume 2: Seminumerical Algorithms",
3rd ed., sec 4.7 "Manipulation of power series", p 526.

The main method (that can be used and for other functions) is firstly to align and represent the formal power series to the from with non-zero first element.

E.g.

x**4/4 + x**6/6 + x**8/8 + ... ==> x**4/4 * (1 + 2*x**2/3 + x**8/2 + ...)

(The shifting operation of the sequence to the left is used.)

Then we exponent the common term x**4/4 and exponent aligned sequence.

For those transformation the helper methods are present:

>>> a = PowerSeries(x, sequence=Sequence((2, oo), 'a'))
>>> a
x**2*a[2] + x**3*a[3] + x**4*a[4] + x**5*a[5] + ...

>>> a.sequence.order
2
>>> c = a**2
>>> c [0]
0

>>> c[4]
x**4*a[2]**2

>>> c[5]
2*x**5*a[2]*a[3]

>>> c.coeff(6)
2*a[2]*a[4] + a[3]**2

For the nested functions, see Faà di Bruno's Formula.

Formal PowerE series

This is similar to the formal power series, but the coefficients are placed near x**n/n! terms.

>>> from sympy.series import PowerSeries_0E


>>> seq = Sequence((1, oo), formula = (k, S(1)/k))
>>> PowerSeries_0E(x, sequence=seq)
x + x**2/4 + x**3/18 + x**4/96 + x**5/600 + ...

Let's construct the series for hyperbolic sin and cosine as examples and look to the operations with them

>>> tcosh = PowerSeries_0E(x, periodical=(1, 0))
>>> tcosh
1 + x**2/2 + x**4/24 + ...

>>> tsinh = PowerSeries_0E(x, periodical=(0, 1))
>>> tsinh
x + x**3/6 + x**5/120 + x**7/5040 + ...

>>> (tcosh + tsinh)**1000
1 + 1000*x + 500000*x**2 + 500000000*x**3/3 + 125000000000*x**4/3 + ...


>>> (tcosh**2 - tsinh**2)**1000
1 + ...

The first elements are calculated fast because there is no needed to calculate and expand the (sympy expression)**1000.

But the later elements now are slow. This is because of that the code is rough, some calculated straightly and with no optimization.

>>> c = tcosh**2 - tsinh**2
>>> c
1 + ...

Slow (~ 2 sec):
>>> c[100]  #doctest: +SKIP
0

>>> c[300] #doctest: +SKIP
0

Negative power

    >>> tcos = PowerSeries_0E(x, periodical=(1, 0, -1, 0))
    >>> tcos**(-1)
    1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + 277*x**8/8064 + ...

Nested series (Composition of the sequences)

What we must obtain:

>>> from sympy import sin, cos, exp

>>> sin(sin(x)).series(n=12) # doctest: +SKIP
x - x**3/3 + x**5/10 - 8*x**7/315 + 13*x**9/2520 - 47*x**11/49896 + O(x**12)

>>> cos(sin(x)).series(n=12) # doctest: +SKIP
1 - x**2/2 + 5*x**4/24 - 37*x**6/720 + 457*x**8/40320 - 389*x**10/172800 + O(x**12)

>>> exp(sin(x)).series(n=12)  # doctest: +SKIP
1 + x + x**2/2 - x**4/8 - x**5/15 - x**6/240 + x**7/90 + 31*x**8/5760 + ...
>>> tsin = PowerSeries_0E(x, periodical=(0, 1, 0, -1))
>>> ns = tsin.compose(tsin)
>>> ns
x - x**3/3 + x**5/10 - 8*x**7/315 + ...

>>> ns[11]
-47*x**11/49896

>>> ns[11:13]
-47*x**11/49896 + 15481*x**13/97297200 + ...

>>> tcos = PowerSeries_0E(x, periodical=(1, 0, -1, 0))
>>> texp = PowerSeries_0E(x, periodical=(1))

>>> tcos.compose(tsin)
1 - x**2/2 + 5*x**4/24 - 37*x**6/720 + ...

>>> texp.compose(tsin)
1 + x + x**2/2 - x**4/8 - x**5/15 - x**6/240 + x**7/90 + ...

Reversion (Composition inversion) of power series

It is inversion operation of compose method. That is we can obtain inverse function. E.g. arcsin from sin:

Target result:

>>> asin(x).series(n=10))  # doctest: +SKIP
x + x**3/6 + 3*x**5/40 + 5*x**7/112 + 35*x**9/1152
>>> t = PowerSeries_0E(x, periodical=(0, 1, 0, -1))
>>> t.reverse()
x + x**3/6 + 3*x**5/40 + 5*x**7/112 + ...

>>> t.reverse().compose(t)  # arcsin(sin(x)) == x
x + ...

The inverse function of the inverse function is the original function:

>>> t.reverse().reverse()  # sin(x)
x - x**3/6 + x**5/120 - x**7/5040 + ...

Converting PowerSeries_0E to PowerSeries and vice versa.

Lets define some PowerSeries_0E:

>>> ts = PowerSeries_0E(x, periodical=(0, 1))
>>> ts
x + x**3/6 + x**5/120 + x**7/5040 + ...

>>> type(ts)
<class 'sympy.series.power_0e.PowerSeries_0E'>

>>> pprint(ts.sequence)
 [0, 1, 0, 1, 0, 1, 0, ...]

and convert to PowerSeries object.

>>> ps = ts.to_power_series()

This object has another type, another bases and internal sequence, but the calculated terms (series) are the same:

>>> ps
x + x**3/6 + x**5/120 + x**7/5040 + ...

>>> type(ps)
<class 'sympy.series.power_0.PowerSeries0'>


>>> pprint(ps.sequence)
[0, 1, 0, 1/6, 0, 1/120, 0, ...]

Some more examples

Some examples (pdf) of the formal power series with abstract sequences, abstract functions and diff operators.

Restrictions

  • Only 1D (one variable) for series expansion can be used with 1D sequences. To make the expansion for tw variables it is need to use another structure. E.g. dictionary of terms or Mario Pernici Lpoly See also wikipedia about Power series in several variables

Critics of the implementation

Further road map

Clone this wiki locally
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.
Press h to open a hovercard with more details.