-
Notifications
You must be signed in to change notification settings - Fork 0
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 PowerSeriesThis 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.
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)
[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 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)
[a[0], a[1], a[2], a[3], a[4], a[5], a[6], ...]
>>> 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
FalseBy 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!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*ySo 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)
[a[0], 0, 2*a[2], a[3], 0, 2*a[5], a[6], ...]
>>> c[11]
2*a[11]factorization
>>> a = Sequence('a')
>>> b = Sequence(formula=(k, S(1)/factorial(k)))
>>> c = SeqMulEW(a, b)
>>> pprint(c)
a[2] a[3] a[4] a[5] a[6]
[a[0], a[1], ----, ----, ----, ----, ----, ...]
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]**2Note, 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))
[a[2], a[3], a[4], a[5], a[6], a[7], a[8], ...]
>>> pprint(a.shiftright(2))
[0, ..., a[0], a[1], a[2], a[3], a[4], ...]
>>> pprint(a.shiftright(2).shiftleft(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.
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] + ...>>> 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 + ...>>> 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]**2For the nested functions, see Faà di Bruno's Formula.
This is similar to the formal power series, but the coefficients are placed near
x**n/n! terms.
>>> from sympy.series import PowerESeries
>>> seq = Sequence((1, oo), formula = (k, S(1)/k))
>>> PowerESeries(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 = PowerESeries(x, periodical=(1, 0))
>>> tcosh
1 + x**2/2 + x**4/24 + ...
>>> tsinh = PowerESeries(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 >>> tcos = PowerESeries(x, periodical=(1, 0, -1, 0))
>>> tcos**(-1)
1 + x**2/2 + 5*x**4/24 + 61*x**6/720 + 277*x**8/8064 + ...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 = PowerESeries(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 = PowerESeries(x, periodical=(1, 0, -1, 0))
>>> texp = PowerESeries(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 + ...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 = PowerESeries(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 + ...Lets define some PowerESeries:
>>> ts = PowerESeries(x, periodical=(0, 1))
>>> ts
x + x**3/6 + x**5/120 + x**7/5040 + ...
>>> type(ts)
<class 'sympy.series.power_e.PowerESeries'>
>>> 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.PowerSeries'>
>>> pprint(ps.sequence)
[0, 1, 0, 1/6, 0, 1/120, 0, ...]- 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