GSoC Week 11

This week I primarily worked on some bugs and added singular initial conditions to the result of expr_to_holonomic().

Firstly I fixed a bug in .unify(). There were some errors being raised in it when one of the Holonomic Function had ground domain RR or an extension of it. This now works fine.

In [9]: expr_to_holonomic(1.4*x)*expr_to_holonomic(a*x, x)
Out[9]: HolonomicFunction((-2.0) + (1.0*x)*Dx, x, 0, {2: [1.4*a]})

In [10]: _9.to_expr()
Out[10]:
       2
1.4⋅a⋅x

Later I fixed a bug in converting the two types of initial condition into one another. Apparently I forgot to add the factorial term while converting.

After that I added singular initial conditions to the result of expr_to_holonomic()  whenever the Indicial equation have one root. For example:

In [14]: expr_to_holonomic(x*exp(x))
Out[14]: HolonomicFunction((-x - 1) + (x)*Dx, x, 0, {1: [1]})

In [15]: _.to_expr()
Out[15]:
   x
x⋅ℯ

I also changed printing of the class HolonomicFunction to include the initial conditions inside the call, so as to make it proper python. These are implemented in #11480 and #11451.

Right now we are trying to find a way to include convergence conditions in the result while converting Holonomic Functions to expressions.

Advertisement

GSoC Week 10

Started off this week by continuing my work on #11422. I added support for singular initial conditions in multiplication and made some amendments in addition too. They now can return a Holonomic function with singular initial condition. The input functions can have singular initial condition both or one of them can have singular one and the other one with ordinary initial condition.

# one function have singular initial condition and the other have ordinary.
In [4]: expr_to_holonomic(x) + expr_to_holonomic(sqrt(x))
Out[4]: HolonomicFunction((1/2) + (-x/2)Dx + (x**2)Dx**2, x), {1/2: [1], 1: [1]}

In [5]: _4.to_expr()
Out[5]: √x + x

In [6]: expr_to_holonomic(x) * expr_to_holonomic(sqrt(x))
Out[6]: HolonomicFunction((-3/2) + (x)Dx, x), {3/2: [1]}

In [7]: _6.to_expr()
Out[7]:
 3/2
x

# both have singular initial conditions.
In [9]: expr_to_holonomic((x)**(S(1)/3)) + expr_to_holonomic(sqrt(x))
Out[9]: HolonomicFunction((1/6) + (x/6)Dx + (x**2)Dx**2, x), {1/3: [1], 1/2: [1]}

In [10]: _9.to_expr()
Out[10]:
3 ___
╲╱ x  + √x
In [11]: expr_to_holonomic((x)**(S(1)/3))*expr_to_holonomic(sqrt(x))
Out[11]: HolonomicFunction((-5/6) + (x)Dx, x), {5/6: [1]}

In [12]: _11.to_expr()
Out[12]:
 5/6
x

I found some problems in coding because of storing these initial conditions in two different attributes. So I merged them to a single attribute and instead added methods to check which one is stored and refactored the existing code using it.

I opened a new PR #11451 majorly focused on adding singular initial conditions to the result of .expr_to_holonomic() when necessary. At first I added it in converting polynomials. Here are some examples:

In [14]: expr_to_holonomic(3*x**3+4*x**2)
Out[14]: HolonomicFunction((-9*x - 8) + (3*x**2 + 4*x)Dx, x), {2: [4, 3]}

In [15]: _14.to_expr()
Out[15]:
 2
x ⋅(3⋅x + 4)

In [16]: expr_to_holonomic(x)
Out[16]: HolonomicFunction((-1) + (x)Dx, x), {1: [1]}

In [17]: _16.to_expr()
Out[17]: x

I also a found a bug in .to_hyper() when the recurrence relation has order 0. Added its fix too. Earlier the output also considered negative exponents which weren’t needed.

# previously
In [18]: expr_to_holonomic(y*x, x).integrate(x).to_expr()
Out[18]:
 2
x ⋅y
──── + C₁
 2
# after fix
In [19]: expr_to_holonomic(y*x, x).integrate(x).to_expr()
Out[19]:
 2
x ⋅y
────
 2  

What Next:

I hope to add singular initial conditions to more types of functions in .expr_to_holonomic().
 

GSoC Week 8 and 9

I couldn’t write a blog post last week so including progress of week 8 and 9 both here.

Week 8

I continued working on the PR #11360. We added functionality to store a different type of initial condition for regular singular points other than the usual [y(0), y'(0), ...]. The exact format is described here in master, though it is changed to a more elegant form in #11422. This type of initial condition provides more information at regular singular points and is helpful in converting to expressions. Examples on how to use it:

In [22]: expr_to_holonomic(sin(x)/x**2, singular_ics={-1: [1, 0, -1]}).to_expr()
Out[22]:
sin(x)
──────
   2
  x

I also added method to compute this type of initial condition for algebraic functions of the form P^r, for some Polynomial P and a Rational Number R.

In [25]: expr_to_holonomic(sqrt(x**2+x))
Out[25]: HolonomicFunction((-x - 1/2) + (x**2 + x)Dx, x), {1/2: [1]}

In [26]: _25.to_expr()
Out[26]:
     _______
√x⋅╲╱ x + 1

After that I made some changes in `to_meijerg()` to return the polynomial itself if the `meijerg` function represents a polynomial instead of raising `NotImplementedError`.

In [40]: expr_to_holonomic(4*x**3 + 2*x**2, lenics=3).to_meijerg().expand()
Out[40]:
   3      2
4⋅x  + 2⋅x

I also added code to return the general solution in `_frobenius()` if none of the roots of indicial equation differ by an integer.

Week 9

I wasn’t able to do much this week because my college started. I travelled back and had some college related stuff to do.

I opened #11422 and first added a basic method to determine the domain for polynomial coefficients in the differential equation.

In [77]: expr_to_holonomic(sqrt(y*x+z), x=x, lenics=2).to_expr()
Out[77]:
  _________
╲╱ x⋅y + z 

In [78]: expr_to_holonomic(1.1329138213*x)
Out[78]: HolonomicFunction((-1.1329138213) + (1.1329138213*x)Dx, x), f(0) = 0

Then I added support for the new type of initial condition on regular singular points in .integrate().

In [83]: expr_to_holonomic(sin(x)/x**3, singular_ics={-2: [1, 0, -1]}).integrate(x).to_expr()
Out[83]:
 ⎛ 2                          ⎞
-⎝x ⋅Si(x) + x⋅cos(x) + sin(x)⎠
────────────────────────────────
                 2
              2⋅x               

Also added support for the same in addition.

In [6]: expr_to_holonomic(sqrt(x)) + expr_to_holonomic(sqrt(2*x))
Out[6]: HolonomicFunction((-1/2) + (x)Dx, x), {1/2: [1 + sqrt(2)]}

In [7]: _6.to_expr()
Out[7]: √x⋅(1 + √2)

I plan to continue my work on this PR and add more support for this initial condition.

GSoC Week 7

I started working on the support for additional symbolic parameters in the module. So I added an extra argument in the method expr_to_holonomic for using custom domains. This will be helpful in integrations. For instance:

In [5]: expr_to_holonomic(sin(x*y), x=x, domain=QQ[y]).integrate(x).to_expr()
Out[5]:
-cos(x⋅y) + 1
─────────────
     y
In [15]: expr_to_holonomic(log(x*y), x=x, domain=QQ[y]).integrate((x, 1, x)).to_expr()
Out[15]: x⋅log(x) + x⋅log(y) - x - log(y) + 1

Code upto here was merged at #11330.

After that I implemented the Frobenius method to support functions where the series may have negative or fractional exponents. Although there are limitations in the algorithm. First one being that we don’t have an algorithm to compute the general solution. There doesn’t seem to be a direct algorithm if roots of the indicial equation differ by an integer. Second one is the initial conditions. Initial conditions of the function (most of the times they won’t even exist) can’t be used if the exponents are fractional or negative.

In [23]: expr_to_holonomic(sqrt(x**2-x)).series()
Out[23]:
            3/2       5/2       7/2         9/2         11/2
        C₀⋅x      C₀⋅x      C₀⋅x      5⋅C₀⋅x      7⋅C₀⋅x        ⎛ 6⎞
C₀⋅√x - ─────── - ─────── - ─────── - ───────── - ────────── + O⎝x ⎠
           2         8         16        128         256            

In [28]: expr_to_holonomic(cos(x)**2/x**2, initcond=False).series()
Out[28]:
         2         4                       3         5
     C₂⋅x    2⋅C₂⋅x    C₁   2⋅C₁⋅x   2⋅C₁⋅x    4⋅C₁⋅x    C₀    ⎛ 6⎞
C₂ - ───── + ─────── + ── - ────── + ─────── - ─────── + ── + O⎝x ⎠
       3        45     x      3         15       315      2
                                                         x

I will be working more on this and hope to add more functionality in this method.

Right now I am working on an algorithm to convert a given Holonomic Function to Meijer G-function. This also doesn’t look so straightforward. Kalevi helped me with the theory and we wrote a basic that works for some cases. We hope to add more things in it.

In [17]: hyperexpand(expr_to_holonomic(exp(x)).to_meijerg())
Out[17]:
 x
ℯ
In [19]: hyperexpand(expr_to_holonomic(log(x)).to_meijerg()).simplify()
Out[19]: log(x)

These are at #11360.

GSoC: Week 6

This week was majorly focused on fixing issues. Thanks to Ondrej, Aaron and Kalevi for testing the module and pointing them out. Some of the issues required a new functionality to be added and some were just plain technical bugs.

I started by writing a method which allows the user to change the point x0 where the initial conditions are stored. Firstly it tries converting to SymPy and then convert back to holonomic i.e. from_sympy(self.to_sympy(), x0=b). If the process fails somewhere then the method numerically integrates the differential equation to the point b. This method was then used in adding support for initial conditions at different points in addition and multiplication. Issues are  #11292 and #11293.

After that I implemented differentiation for a holonomic function. The algorithm is described here. A limitation is that sometimes it’s not able to compute sufficient initial conditions (mainly if the point x0 is singular).

There was also a bug in to_sequence() method. Apparently, whenever there weren’t sufficient initial conditions for the recurrence relation,  to_sympy() would fail. So I changed it to make to_recurrence() return unknown symbols C_j , where C_j representing the coefficient of x^j in the power series. So now even if we don’t have enough initial conditions, to_sympy() will return an answer with arbitrary symbols in it. This is also useful in power series expansion.

Earlier the method to_sympy() would only work if initial conditions are stored at 0 which is a big limitation. Thanks to Kalevi for the solution we now can find the hyper() representation of a holonomic function for any point x0 .(of course only if it exists.) This is further used by to_sympy() to convert to expressions. There is also something interesting Kalevi said “The best points to search for a hypergeometric series are the regular singular points”. I observed this turned out to be true a lot of times. We should keep this in mind when using to_sympy().

Later I added functionality to convert an algebraic function of the form p^(m/n) , for a polynomial p . Thanks to Aaron and Ondrej for the solution. Also added custom Exception to use in the holonomic module. Some small bugs #11319, #11316 and 11318 were also fixed.

Currently I am trying to make from_sympy() work if additional symbols are given in the expression and also support other types of fields (floats, rationals (default right now), integers  ). This involves extending the ground domain of the Polynomials used internally. Hopefully this should be done in the next couple of days.

The unmerged fixes are at #11330.

 

GSoC Week 5

Earlier this week, I continued my work on the method to convert a given Holonomic Function to Hypergeometric.

There were some major bugs that produced wrong results because of the fact that the recurrence relation doesn’t holds for all non-negative n. So I changed it to also store the value, say n0, so the recurrence holds for all n >=n0. I changed the method .to_hyper() accordingly and it works pretty fine now.

I also added the method to convert to symbolic expressions to_sympy() . It uses hyperexpand , however will only work for cases when a hyper representation is possible.

These were implemented at [#11246].

After that I have been majorly working on some issues Ondrej pointed out. These are #11292, #11291, #11287, #11285, #11284. The fixes are implemented at [#11297].

This week I will be fixing some more bugs, the ones listed at [Label: Holonomic Functions] and also a bug discussed with Kalevi on Gitter.

Cheers!

GSoC Week 4

I have been working on two things since my last blog post and they are: converting Symbolic Functions/Expressions to Holonomic and converting Holonomic to Hypergeometric. Let’s discuss the former first:

I wrote a method for converting Expressions earlier too, but that was very trivial, so this week began by me writing a Look-Up table technique to convert expressions. It is similar to the way meijerint._rewrite1 converts to meijerg in the integrals module. The method recursively searches for chunks in the expression that can be converted to Holonomic (using tables) and then uses the operation +, *, ^ defined on Holonomic Functions to convert the whole expression. #11222

Then I started working on to write a method that converts Holonomic Functions to Hypergeometric or a linear combination of them. (if it’s possible, because every Holonomic Functions is not Hypergeometric.) First I asked questions I had to Kalevi on Gitter about how to do it, and then finally wrote the code here #11246.

Let me also give you an overview of things we have in the module so far.

• Differential Operators and Recurrence Operators.
• Holonomic Functions with following operations: addition, multiplication, power, integration and composition).
• Computing Recurrence relation in Taylor coefficients.
• Series expansion using the recurrence. (Some know issues are there if origin is a regular singular point.)
• Numerical Computation (Runge-Kutta 4th Order and Euler’s method.)
• Converting to Hypergeometric. (Currently working on this.)
• Converting Hypergeometric to Holonomic.
• Converting Meijer G-function to Holonomic.
• Converting expressions to Holonomic. (There doesn’t seem to be a direct algorithm for converting algebraic functions which are not polynomials or rationals, so they are not supported right now.)

This week I will be continuing my work on converting to hyper and then further extending it to convert to expressions using hyperexpand.

Good luck to all fellow GSoCers for Midterm Evaluation.

Goodbye!

GSoC: Week 3

As I wrote in my blog post last week, Our plan was to implement numerical methods to compute values of a Holonomic function at given points and convert symbolic function/expressions. Let’s look at what I implemented during this time one by one.

I started by working on to implement the most basic and popular Euler's method for numerical integration of holonomic ode’s. Given a set a points in the complex plane, the algorithm computes numerical values of Holonomic function at each point using the Euler’s method. #11180

Next I moved on to write a method for converting Meijer G-function to a Holonomic Function. In the module meijerint we have method meijerint._rewrite1 to convert expressions/functions to meijerg . Now this method can be used first to convert to meijerg and then converting this result finally to Holonomic Functions. I wrote methods from_meijerg and from_sympy respectively for this purposes. #11187

We are also familiar with the fact that Euler's method doesn’t give much accurate results and for better accuracy one needs higher order numerical methods e.g. Runge-Kutta 4th Order method A.K.A RK4 method. I wrote this method and made it the default method for numerical computation. As expected, the results are much more accurate using RK4 and enough for our present purposes. #11195

I am now looking for bugs inside everything implemented so far in the module, will be working on solving them and also will make things optimal wherever there is a possibility. More on what to implement next is yet to be discussed.

Thank You!

 

GSoC: Week 2 Begins

First, I am going to talk about what has been done since last week. I implemented methods to find a recurrence relation in the coefficients of Power series expansion of a Holonomic Function, and then method to find the Power series.

One of the many interesting properties of Holonomic Functions is that the coefficients in Power series expansion of a Holonomic function is a Holonomic sequence, i.e. these coefficients satisfy a recurrence relation having polynomial coefficients.

This recurrence relation can be used to efficiently calculate the coefficients of Power series. The PR 11153 had these methods and got merged today. Here are a couple of examples of the implementation.

In []: p = HolonomicFunction(Dx – 1, x, 0, [1])  # exp(x)
In []: p.to_sequence()

Out []: HolonomicSequence((-1) + (n + 1)Sn, n), u(0) = 1

In []: p.series()

Out []: 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)

For this week it is planned first to implement a method to numerically integrate differential equations of holonomic type from any point x=a to x=b in the complex plane. The next thing to implement after this would probably be converting symbolic functions/expressions to Holonomic Functions.

Cheers and Happy Coding.

GSoC: First Week of Coding Phase

Since my last blog post I have opened two new PR’s for the Project. We’ll see what they do and discuss the goals of this week.

A major issue in the first PR was the slow performing algorithms, a consequence of using recursive expressions internally. Most of the algorithm implemented used matrices module to solve the linear system which doesn’t support DMP and DMF objects.

So I defined a new class subclassed from MutableDenseMatrix and changed some methods to make it work with Polynomials and Fractions and used this to use Polynomials internally. Thanks to Kalevi for this idea. It works much more robust now. I have also added methods to find composition of Holonomic Functions and converting a Hypergeometric Function to Holonomic. These things are added in this PR. I hope the PR gets merged in a couple of days.

A new PR was opened for features relating to recurrence relations in coefficients of Power Series expansion of Holonomic Functions. The first thing I did was defined a class RecurrenceOperator parallel to DifferentialOperator to store the recurrence relation.

Goals of the Week:

In this week, I have planned to define a function to find the Recurrence Relation of series coefficients and then go for numerical computation of Holonomic Functions. Let me know If anything else should be implemented first as I haven’t discussed this with mentors yet.

The chronology might be different from what I wrote in the Proposal but we are quite ahead of that.

Cheers Everyone.