UD Sequences and formal power series prototype

UD - Sequences and formal power series - prototype

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
>>> from sympy.abc import x, k
>>> from sympy.series.sequences import Sequence
>>> from sympy.series.sequencesexpr 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.series.sequences import Sequence
>>> from sympy.printing.pretty.pretty import pprint

>>> a = Sequence((0, oo), 'a')
>>> pprint(a)
[a[0], a[1], a[2], a[3], a[4], a[5], a[6], ...]

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, ..., a[3], a[4], a[5], a[6]]

>>> b = Sequence((5, 9), 'b')
>>> pprint(b)
[0, ..., b[5], b[6], b[7], b[8], b[9]]

>>> (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((0, oo), 'a')
>>> a
a

>>> pprint(a)
[a[0], a[1], a[2], a[3], a[4], a[5], a[6], ...]

>>> a[2]
a[2]

Periodical sequence:

>>> seq = Sequence((0, oo), 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.series.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((0, oo), 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 scalar expression:

>>> 2*a
2*SeqPer([0, oo), (1, 0))

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

>>> (3*a).coefficient
3

>>> 2*(3*a).coefficient
6

So we can use the linear combination of them:

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

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

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

>>> pprint(c)
[a[0], 0, 2*a[2], a[3], 0, 2*a[5], a[6], ...]

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

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 main method for binary '*' operation:

>>> a = Sequence((0, oo), 'a')
>>> b = Sequence((0, oo), '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 as the number of Cauchy product by itself.

>>> a = Sequence((0, oo), '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((0, oo), 'a')
>>> pprint(a.toleft(2))
[a[2], a[3], a[4], a[5], a[6], a[7], a[8], ...]

>>> pprint(a.toright(2))
[0, 0, a[0], a[1], a[2], a[3], a[4], ...]

>>> pprint(a.toright(2).toleft(2))
[a[0], a[1], a[2], a[3], a[4], a[5], a[6], ...]

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((0, oo), periodical=(1, 0)))
>>> a
1 + x**2 + x**4 + x**6 + x**8 + ...

Summation of to series is obviously element-wise:

>>> a = PowerSeries(x, sequence=Sequence((0, oo), periodical=(1, 0)))
>>> b = PowerSeries(x, sequence=Sequence((0, oo), 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, sequence=Sequence((0, oo), 'a'))
>>> b = PowerSeries(x, sequence=Sequence((0, oo), '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 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 respresent 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.first_nonzero_n
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 Taylor series

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

>>> from sympy.series import Sequence, TaylorSeries

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

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

>>> a = Sequence((0, oo), periodical = (1, 0))
>>> tcosh = TaylorSeries(x, sequence=a)
>>> tcosh
1 + x**2/2 + x**4/24 + ...

>>> b = Sequence((0, oo), periodical = (0, 1))
>>> tsinh = TaylorSeries(x, sequence=b)
>>> 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 and restricted by stack. This is because of that the code is rough, some calculated straightly and with no optimization.

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

Slow:
>>> c[100]      # doctest: +SKIP
0

Restrictions

Critics of the implementation

Provide feedback

Saved searches

Use saved searches to filter your results more quickly