GSoC 2014 Application Avichal Dayal Series Expansion - sympy/sympy GitHub Wiki

Personal Details:

My name is Avichal Dayal and I am a second year Computer Science student at International Institute of Information Technology, Hyderabad.

Email : [email protected]

IRC : avichal on freenode

Skype : avichaldayal

Github : avichaldayal

Programming and Mathematical background

I have done coding in various high-level languages (C, C++, Python) with projects in Python and Python-related technologies(web2py). These include Data Visualization Comparator(analyzing images), signal analysis and web-applications.

I was introduced to the open source community by my college's open source development group and I am amazed by the open source community. I have used svn before and I have become familiar with git by participating in SymPy's workflow. I have been using Linux for two years now.

Math has been a strong point for me. I have taken relevant courses in Math like Linear and Abstract algebra and I am quite familiar with the concepts.

Contribution:

Unmerged:

  • #2905 This involves implementing asymptotic expansion of Fresnel integrals. Task is completed and I think reviewing is almost done.

  • #2706 Enhanced degree function in polynomial to work for constants without generators.

Merged:

  • #2812 Passed logx parameter in leadterm to avoid expanding log(x) term

  • #2717 Fixed bug by changing force=False to True when combining log

  • #2690 Removed duplicate output from solvers by passing it through set

Closed:

  • #2899 This aimed to improve residue but I realized the need to implement asymptotic expansion in SymPy. I discuss this in my ideas.

The Project:

  • Formal Power Series:

Other CAS are capable of returning an infinite series in terms of a generating function or a sequence. I feel we should have this. Currently series() returns something like 1 + x + x**2/2 + O(x**3) while Mathematica and others are also capable of returning something like summation(x**n/factorial(n), (n, 0, oo)).

Examples:

In[1] := FormalSeries(exp(x), x, 0).as_summation()
Out[1] := Summation(x**k/k!, (k, 0, oo))

In[2] := FormalSeries(cos(x), x, 0).as_summation()
Out[2] := Summation((-1)**k * x**(2*k) / (2*k)!, (k, 0, oo))
  • Asymptotic series:

Most of the routines are implemented in gruntz.py. However still work is to be done to implement it. Series is extensively used by gruntz and asymptotic expansion is not implemented in SymPy. This would improve gruntz greatly.

Examples:

In[1] := series(sin(1/x+exp(-x)) – sin(1/x), x, oo)
Out[1] := (1 – 1/(2*x**2) +1/(24*x**4) + O(1/x**6)) * exp(-x)

In[2] := series((exp(exp(x)) * (exp(sin(l/x+l/exp( exp(x))))-exp(sin(1/x))), x=oo, 3)
Out[2] := 1 + 1/x - 1/(2*x**3) + O(1/x**4)

Overall aim of the project:

  • Infinite series from sequence, differential equation, recurrence equation or a function
  • Convert function to summation (as a by product from above point)
  • Play with infinite series (addition, multiplication, invert, composition etc.)
  • Implement asymptotic series (MrvAsympt algorithm)
  • Special functions series implementation

Timeline: (tentative)

  • Week 1: Basic class structure for FormalSeries and Sequence
  • Week 2: Implement periodical, formula and recurrence sequences
  • Week 3: Add test-cases for FormalSeries with sequences
  • Week 4: Implement SimpleDE and DE-to-RE algorithms
  • Week 5: Implement solveRE and add sophistcated test-cases
  • Week 6: Implement rational algorithm and start implementing operations on FPS.
  • Week 7: Implement operations on FPS – addition, multiplication, convolution, inversion, composition
  • Week 8: Make FormalSeries work at points other than zero. Solve bugs if present and make necessary changes. Add test-cases for operations on FormalSeries
  • Week 9: Implement MrvAsymptotic algorithm
  • Week 10: Implement asymptotic expansion for special functions
  • Week 11: Implement rewrite_as_tractable functions where-ever possible
  • Week 12: Add test cases for asymptotic series. This should cover examples in gruntz as well as series
  • Week 13: Documentation, docstrings and stylistic changes where necessary

Formal Power Series

Basic class structure: <--FormalSeries-->

Attributes:

  • Attribute 1: Sequence
In[1] := FormalSeries(seq=(1, 2, 3), x)
Out[1] := 1 + 2*x + 3*x**2 + x**3 + 2*x**4 + 3*x**5 + x**6 + ...
  • Attribute 2: Formula
In[1] := FormalSeries(formula=(1/k, k), x)
Out[1] := 1 + x + x**2/2! + x**3/3! + x**4/4! + ...
  • Attribute 3: Recurrence equation
  • Attribute 4: Differential equation
  • Attribute 5: Function
Algorithm described below mentions how we can convert DE, RE or a function to a formula

Methods:

  • Method 1: simpleDE()
Converts a function to a simple differential equation
  • Method 2: DE-to-RE()
Converts a differential equation to recurrence equation
  • Method 4: solveRE()
Solves recurrence equation to find a formula
  • Method 5: as_summation()
Writes the infinite series in summation form
  • Method 6: Add
  • Method 7: Multiply
  • Method 8: Convulate
  • Method 9: Invert
  • Method 10: Shift
  • Method 11: Composition
Various operations on series

Representation of infinite series (important):

Reference[1] talks about how we can represent infinite series and how operations can be done on it. Infinite series are usually represented by a generating function or a formula of some sort. That's the advantage compared to truncated series. There is no need to store the coefficients, they can be generated on the fly.

Keeping in mind multi-variable series and also the latest development in series module i.e. Lpoly1, I plan to store it in the form c*X**a where X**a is a monomial in several variables. See the wikipedia link which describes the same method to represent the generating functions: Wikipedia link With this there is no other special structure introduced like generating function, sequence etc.

1) Finding DE of degree k (simpleDE)

We need to find DE of the following form: summation(A(j) * f(x).diff(x, j)), (j, 0, k)) where A(j) is a rational function in x.

(a) Keep an upper bound (papers mention k=4 is safe option)
(b) Set n=1 and form the DE.
(c) Solve for A(j)'s in the DE. If all are rational then break, else increase n by 1

Example(step by step):-

>>> simpleDE(exp(x) * sin(x))
… For n = 1, A(0) = 1+cot(x) which is not rational in x
… Try n = 2, DE = 2*exp(x)*cos(x) + A(1)*(exp(x)*(sin x + cos x)) + A(0)*exp(x)*sin(x )
… Solving we get A(0) = 2, A(1) = -2 which are rational in x

Example(from references):-

In[1] := simpleDE(((1+x)/(1-x))**n, x)
Out[1] := 2*n*f(x) + (x-1)*(x+1)*f(x).diff(x, 1)

In[2] := simpleDE(asin(x), x)
Out[2] := x*f(x).diff(x, 1) + (x-1)*f(x).diff(x, 2)

2) Convert DE to RE

This is a rather easy step. In the DE, we just need to substitute

f.diff(x, k) * x**j -> pochhammer(n+1-j, k) * a(n+k-j)

This is proved in Reference[3]. I didn't understand the proof completely but I don't think it will affect my work.

Examples(from references):

In[1] := DEtoRE(2*n*f(x) + (x-1)*(x+1)*f(x).diff(x, 1), x)
Out[1] := k*a(k) + 2*n*a(k+1) - (k+2)*a(k+2)

3) Solving RE

First part is to recognize the recurrence equation.

(a) If RE has constant coefficients then subsituting b(n) = a(n)*n! does the work. The corresponding equation can then be solved with rsolve
(b) If RE contains only one summand on the right hand then it's of hypergeometric type which has known hypergeometric series corresponding to it. (reference [2] page-2)

4) Rational Algorithm

This is a kind of heuristic which finds a generating function and an infinite series for a rational function f.

Given f, we apply a partial fractional decomposition to separate into terms of form c/(x-a)**n. These terms can then be expanded separately using binomial theorem. SymPy already has the apart function which makes this task easy.

5) Asymptotic

At points other than zero: FormalSeries at points other than zero can be done by shifting similar to what series does. Also asymptotic expansion at oo or singularity points can be done using the same algorithm. Cases where singularity exist give an error.

6) Special functions

Also I was thinking of having formalSeries made separately for special functions similar to what eval_aseries and eval_nseries do. Although this algorithm can find series for error functions, exponential integral functions etc, as well, it will be better to implement for them separately.

7) Operations on FPS

I am aware of the difficulties in doing much with formal power series as of now. Our summation algorithms are not that great which limits a few of the following functionalities. However we can still implement it although the final answer will not be in the most simplified form.

Reference[6] i.e. the wikipedia article on formal power series mentions the formula for various operations on FPS.

  • Addition: Addition and subtraction is term-wise

Examples:

>>> ps1 = FormalSeries(seq=(1, 2, 3), x)
>>> ps2 = FormalSeries(seq=(2, 3, 4), x)
>>> ps1 + ps2
3 + 5*x + 7*x**2 + 3*x**3 + 5*x**4 + 7*x**5 + ...
  • Multiplication: If the final form is of Summation(c(n)*x**n, (n, 0, oo)) then c(n) is given as:
c(n) = summation(a(i)*b(n-i), (i, 0, oo))
where a(i) and b(i) are coefficients of ps1 and ps2
  • Inverse
>>> ps1 = Summation(a(n)*x**n, (n, 0, oo))
>>> ps1.invert()
Summation(b(n)*x**n, (n, 0, oo))
where b(n) = -Summation(a(i)*b(n-i), (i, 0, oo))/a(0)
  • Shift
  • Differentiation and Integration: This is term-wise differentiation and integration
>>> ps1 = Summation(a(n)*x**n, (n, 0, oo))
>>> ps1.diff(x, 1)
Summation(a(n)*n*x**(n-1), (n, 0, oo))
  • Composition
  • Convolution

There exists a closed formula for composition of two power series and also to find the coefficients. This link mentions them clearly: http://en.wikipedia.org/wiki/Fa%C3%A0_di_Bruno%27s_formula

These are the ones that the current summation algorithms won't be able to simplify but it's better to implement them nonetheless in my opinion.

Sample run:

>>> formalseries(exp(x)*sin(x), x, 0)
… procedure simpleDE(exp(x)*sin(x)):
… 	Trying n = 1
… 	A*exp(x)*sin(x) + (exp(x)*cos(x)+exp(x)*sin(x)) {does not break as A is not rational in x}
… 	Trying n = 2
… 	returns A*exp(x)*sin(x) + B*(exp(x)*cos(x) + exp(x)*sin(x)) + 2*exp(x)*cos(x) {breaks as A=2, B=-2 are rational in x}
… procedure DE-to-RE(f(x).diff(x, 2) – 2*f(x).diff(x) + 2f)
… 	DE has constant coefficients hence b(n) → a(n)/n!
…	returns b(n+2) – 2*b(n+1) + 2*b(n)
… procedure solveRE(r=b(n+2) – 2*b(n+1) + 2*b(n))
… 	a(n) = rsolve(r, {b(0):0, b(1):1})
…	return a(n)
>>> ans = summation(a(n)*x**n, (n, 0, oo))

Asymptotic Expansion

Relevant issue: Issue 13

Most of the routines to find asymptotic expansion are implemented in gruntz.py which makes the implementation easier.

Example:-
>>> series(sin(1/x+exp(-x)) – sin(1/x)
>>> O(1/x**6) {currently}
Better would be:-
>>> (1 – 1/(2*x**2) +1/(24*x**4) + O(1/x**6)) * exp(-x)

The MrvAsympt algorithm described in [1] page-90 solves this

1) Mrv, rewriting and recursion

Finding the most-rapidly varying term and rewriting expression in terms of omega is already implemented. After rewriting the expression, coefficients are recursively expanded in terms of their most rapidly varying subexpression.

Algorithm is terminated when coefficients become constant. Special handling is done in case recursion never terminates by keeping an upper bound or Heirarchical series.

(a) Determine set of mrv subexpressions
(b) Choose w from above set and rewrite other expressions as g(x) = h(x)*w**p
(c) Compute series around w = 0

Reference [1] and [4] give a detailed explanation and it won't be much trouble implementing this.

2) Heirarchical series

There are cases like exp(1/x + exp(-x**2) * (exp(ax) - exp(bx)) - exp(1/x) whose asymptotic expansion using the above algorithm will result in a complicated expansion. In such cases, it is better to terminate early in the recursion process. Coefficients need not be expanded completely. Also sometimes not expanding coefficients results in a more informative expansion.

In[1] := exp(1/x + exp(-x**2) * (exp(a*x) - exp(b*x))  - exp(1/x)
Out[1] := exp(1/x)*(exp(a*x) - exp(b*x))*exp(-x*x) + exp(1/x)/2*(exp(a*x) - exp(b*x))**2*exp(-x*x)**2 + O(exp(-x*x)**3)

3) Special functions

Relevant issue: Issue 2958

Quite a few special functions don't have their asymptotic as well as eval_nseries method implemented. I already implemented it for fresnel integrals and it won't take long for other functions. Also after this I think various test-cases will need to be added both to series and gruntz. In case expansion does not exist at infinity, functions need to be rewritten in terms of exponential functions.

TODO part of each file describes the not-implemented part thoroughly. Many functions don't have their _eval_nseries and _eval_aseries method implemented. Even if they are, they can be represented in a nicer way. This should be possible to do as already there is enough literature available on the net. I would refer various resources along with [5] for information regarding this.

Notes

I have no other plans in summer and I would dedicate 40 hours a week to complete my work according to my timeline.

Most of the references have pseudo-code given in their paper and I would refer them for my help. I plan to be an active SymPy contributor after GSoC period and continue work on the series module in particular.

References:

[1] On computing limits in a symbolic manipulation world – Dominik Gruntz

[2] "Formal Power Series" by Dominik Gruntz and Wolfram Koepf

[3] Power Series in Computer Algebra - Wolfram Koepf

[4] "A New Algorithm Computing for Asymptotic Series" by Dominik Gruntz

[5] http://mathworld.wolfram.com/

[6] Wikipedia - formal power series

Relevant group discussions and PR