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.