Mfrhtyj

Consider a linear, homogenous 2nd order ODE:
begin{equation}
Lleft[y(x)right] = left[frac{d^2}{d x^2} + a_1(x) frac{d}{d x} + a_0(x)right] y(x)=0
end{equation}
In https://arxiv.org/pdf/1606.01576.pdf the authors present an interesting heuristic algorithm for finding closed form solutions of those ODEs if such solutions exist. The algorithm can be summarized as follows:

1. For all regular singular points of the ODE in question find exponents and power series expansions of solutions to a given order. If the indicial equation is degenerate (which means the second solution cannot be found) find that second solution via the Wronskian method http://mathworld.wolfram.com/Second-OrderOrdinaryDifferentialEquationSecondSolution.html.
   Note that the solutions in question will be of the form:
   begin{equation}
   y_p(x):= x^{nu_p} sumlimits_{i=0}^infty left( a_i^{(0)} + a_i^{(1)} log(x) right)cdot x^i
   end{equation}
   for $p=1,2$.




2. Having found the solutions in question (meaning of course series expansions being truncated to some given order) find a pair of differential operators $({mathcal B}_0, {mathcal B}_1)$ that are globally integral. A precise definition of that notion is given in the paper above and in references therein however for our purposes it amounts to point out the following intuitive meaning. An operator is globally integral if its action on any of the solutions yields a series expansion that comprises nonnegative powers of both $x$ and $log(x)$ only. It turns out (see the paper for details) that such operators have the following form:
   begin{equation}
   {mathcal B}_p= r_0(x) + r_1(x) cdot frac{d}{d x}
   end{equation}
   for $p=1,2$.




3. Choose one of such operators, say ${mathcal B}_1$ and then construct another differential operator ${tilde L}$ such that
   begin{equation}
   {tilde L}left[{mathcal B}_1 y_p(x)right] = 0
   end{equation}
   Note that this new operator ${tilde L}$ is also w linear 2nd order operator.
   Now one would naively expect that this new operator will be more complicated than the original operator $L$. This is indeed the case in general. However if it just happens that $L$ has closed form solutions in terms of special functions then the new operator will be much simpler and will allow us to find those closed form solutions.

Now, I have tested this algorithm on the family of ODEs specified in the first example in my answer to Gauge transformation of differential equations I . That family is a five parameter family . I have therefore generated a random five tuple of such parameters by sampling from rational numbers then for such five tuple I have generated the operator $L$ and then found its global integral basis which I then used to find the new operator $tilde{L}$ .
To my surprise it turns out that the new operator ${tilde L}$ was always much simpler than the original one and I could guess that it corresponds to an abscissa transformed Gaussian hypergeometric operator. Below I enclose a code snippet that I used:

In[929]:= x =.; Dx =.; A =.; B =.; CC =.;

order = 4;



AA = RandomInteger[{1, 10}];

{A2, A1, B2} = Rationalize[RandomReal[{-5, 5}, 3], 1/10];

{A, B, CC} = {1/7, 1/3, 1/2};

{L, mysubst} = GenerateODE[A2, A1, B2, A, B, CC, AA]

g = GlobalIntegralBasis[L];

g1 = NormalizeAtInfinity[L, g]

myStandardForm[L, g1[[1]]]



L00 = ChangeOfVars[{1, (CC - (A + B + 1) x)/(

    x (1 - x)), -((A B)/(x (1 - x)))}, (AA x)^2];

Factor[L00/L00[[1]]]





Out[934]= {{1, -((

   4 (37632 - 52920 x - 3488051 x^2 - 1014720 x^3 + 28451380 x^4 + 

      12779520 x^5))/(

   21 x (-1 + 4 x) (1 + 4 x) (-3584 + 3360 x + 165415 x^2 + 

      49920 x^3))), (

  2 (10752 - 15120 x - 1066847 x^2 - 577920 x^3 + 6640000 x^4 + 

     4792320 x^5))/(

  3 x^2 (-1 + 4 x) (1 + 4 x) (-3584 + 3360 x + 165415 x^2 + 

     49920 x^3))}, 

 b^n$34032_ :> (14 b^(-3 + n$34032))/195 - (7 b^(-2 + n$34032))/

104 - (33083 b^(-1 + n$34032))/9984 /; n$34032 >= 3}



Out[936]= {{-((

   112 ((80 - 15 x - 1792 x^2) f[x] + 

      48 x (-1 + 16 x^2) Derivative[1][f][x]))/(

   x (-3584 + 3360 x + 165415 x^2 + 49920 x^3))), ((1792 + 560 x - 

      11763 x^2 + 16640 x^3) f[x] + 

   1792 x (-1 + 16 x^2) Derivative[1][f][x])/(

  1792 x^2 (-3584 + 3360 x + 165415 x^2 + 49920 x^3))}, {0, 1}}



Out[937]= {1, (656 x)/(21 (-1 + 4 x) (1 + 4 x)), 64/(

 21 (-1 + 4 x) (1 + 4 x))}



Out[939]= {1, (656 x)/(21 (-1 + 4 x) (1 + 4 x)), 64/(

 21 (-1 + 4 x) (1 + 4 x))}

As I said before for the randomly sampled parameter five tuple the code produces firstly the operator $L$ (Out[934]) then it produces the global integral basis (Out[936]) and finally it produces the new operator ${tilde L}$ (Out[937]) along with an abscissa transformed Gaussian hypergeometric operator (Out[939]). I have played with this code and it was always that the last two lines of the output were the same, meaning Out[937] was the same as Out[939] . This means that the algorithm is always able to solve the ODE in question provided of course it has a closed form solution which in this cases it always has. I have also run the algorithm on other ODEs which do not have a closed form solution and in there the algo has failed as expected.

Now having said all this my question would be as follows. Why does this algorithm work? Could anyone enlighten me why an integral basis should be useful as a means for transforming a complicated ODE into a simple one?
Can we prove that this algorithm is going to work if a closed form solution exists?