RiccatiSolve

first order right factors of a linear Mahler operator

Calling sequence:

$\operatorname{RiccatiSolve}(L, x, M, b)$

Parameters:

• $L$, a linear Mahler operator
  
• $x$, the variable
  
• $M$, the name of the Mahler operator
  
• $b$, the radix
  
• $\mathit method$, a keyword parameter to choose the method of computation
  
• $\mathit field\_extension$, a keyword parameter which defines a field extension
  
• $\mathit free$, a keyword parameter which provides the name to parametrize the solutions
  
• $\mathit sigma\_min$, $\mathit sigma\_max$, $\mathit sigma\_ratio$, $\mathit sigma\_shift$, some optional keyword parameters when '$\mathit method$' = '$\mathit hermite\_pade\_approximants$' is chosen.

Description:

• This procedure will be deleted from future versions. It is intended to illustrate the future cited below article by Chyzak et alii. The standard way to compute the solutions of the Riccati equation associated to $L$ uses the procedure $LMOpSolve$ with the option '$\mathit output$' = '$\mathit riccati$'.
  
• The procedure computes the ramified rational solutions of the Riccati equation associated to the operator $L$, that is the ramified rational functions $u$ such that $M - u$ is a right factor of $L$. In other words, if the operator writes $$ L= \ell_0(x) + \ell_1(x) M + \dots + \ell_r(x) M^r,$$ then the Riccati equation writes $$\ell_0(x) + \ell_1(x) u(x) + \dots + \ell_r(x) u(x) u(x^b)\dots u(x^{b^{r-1}}) = 0$$ and the procedure provides the set of ramified rational functions solutions $u(x)$ of it.
  
• More precisely the solutions are given as a set of parametrized ramified rational functions expressions and when the tuples of parameters, say $(K_1 : K_2 : \ldots : K_s)$, vary in the projective space ${\mathbb P}^{s-1}({\mathbb L})$, if $\mathbb L$ is the chosen field extension, the obtained solutions are all distinct and no two parametrizations can give the same solution.
  
• It is assumed that $L$ is a polynomial in $M$ with rational functions coefficients in $x$.
  
• It is not assumed that the valuation of $L$ wrt $M$ is $0$.
  
• The requested solutions have their coefficients in a field extension.
  • It may be defined à la Maple, with '$\mathit field\_extension$' = $\alpha$ or '$\mathit field\_extension$' = $\{\alpha, \beta, \ldots\}$, where $\alpha$, $\beta\ldots$ are expressed by $\operatorname{RootOf}$'s and NOT by radicals (but $\operatorname{RootOf}$'s and radicals can be employed with the procedure $\operatorname{LMOpSolve}$, which is here the standard for solving). With '$\mathit field\_extension$' = $\{\}$, the field of coefficients of $L$ is used.
    
  • The user can also use the algebraic closure of the ground field with '$\mathit field\_extension$' = '$\mathit splitting\_field$'.
    
• The choice '$\mathit field\_extension$' = '$\mathit splitting\_field$'
  • reinforces the exponential complexity of the Petkošek method by increasing the number of factors to be taken into account;
    
  • guarantees the termination and correction of the Hermite-Padé approximation approach.
    
• The default value of the name used to parametrize the solutions is $\mathit \_C$. Any keyword parameter '$\mathit free$' = '$K$' with $K$ a name or a symbol is valid. Then the solutions are expressed with the parameters $K_1$, $K_2\ldots$
  
• The algorithm used to compute the solutions can be specified by the optional keyword parameter method:
  • '$\mathit method$' = '$\mathit basic\_petkovsek$',
    
  • '$\mathit method$' = '$\mathit improved\_petkovsek$',
    
  • '$\mathit method$' = '$\mathit hermite\_pade\_approximants$'. The last one is the default.
    
• In case '$\mathit method$' = '$\mathit hermite\_pade\_approximants$' is used there are supplementary optional keyword parameters (see the cited work by Chyzak et alii for a deeper explanation):
  • The optional keyword parameter $\mathit sigma\_min$ permits to choose the initial accuracy $\sigma$ of the formal series expansions in the search for approximate syzygies between formal power series solutions of the ancillary equations. The default is the parameter $\nu$ in the referenced article about linear Mahler equations (Chyzak et alii, 2018). Whatever the value of $\mathit sigma\_min$, the minimum value used for $\sigma$ is $\geq \nu$.
    
  • The optional keyword parameter $\mathit sigma\_ratio$ gives the multiplier for the value of $\sigma$ from one iteration to the next. By default, the multiplier is the golden ratio $\varphi = 1.618\ldots$ introduced by '$\mathit sigma\_ratio$' = '$\mathit golden\_ratio$'. The other possibility is a real number $q > 1$ which can be numerically evaluated. With '$\mathit sigma\_ratio$' = $q$, the tested values of $\sigma$ are a geometric sequence $\mathit sigma\_min$ * $q^k$, $k \geq 0$.
    
  • The optional keyword parameter $\mathit sigma\_shift$ takes precedence over the optional keyword parameter $\mathit sigma\_ratio$. With '$\mathit sigma\_shift$' = $q$, a positive integer, the tested values of $\sigma$ are an arithmetic sequence $\mathit sigma\_min$ + $q\times k$, $k \geq 0$.
    
  • The optional keyword parameter $\mathit sigma\_max$ is an upper bound for $\sigma$. If the value of $\sigma$ is larger than $\mathit sigma\_max$, the computation process stops and $\mathit FAIL$ is returned.
    
• In case '$\mathit hermite\_pade\_approximants$' is used, the procedure may fail, and then returns $\mathit FAIL$, for three reasons:
  • the internal computation of syzygies may fail because Maple is not able to deal with algebraic numbers;
    
  • the option '$\mathit field\_extension$' = '$\mathit splitting\_field$' is not used and the number of iterations on $\sigma$ becomes too large;
    
  • $\sigma$ becomes larger than $\mathit sigma\_max$.
    
• In case '$\mathit hermite\_pade\_approximants$' is used, the procedure is currently based on the procedure absPrimdecGTZ of Singular, which has to be installed.

References:

• Frédéric Chyzak, Thomas Dreyfus, Philippe Dumas, and Marc Mezzarobba (2018). Computing solutions of linear Mahler equations. Mathematics of Computation 87.314, pp. 2977–3021.
  
• Roques Julien (2018). On the algebraic relations between Mahler functions. Trans. Amer. Math. Soc. 370.1, pp. 321--355.
  
• Frédéric Chyzak, Thomas Dreyfus, Philippe Dumas, and Marc Mezzarobba. First-order factors of linear Mahler operators. In preparation.

Example:

libname := libname, FileTools:-JoinPath(["maple","lib","dcfun.mla"],base=homedir):

 with(dcfun):

 b := 3;

$$3$$

We tinker an example using a least common left multiple of first order operators.

 L1 := RamRatPolyToLMOp(1/(1 + 2*x), x, M, b);

$$1+2 x -\left(2 x^{3}+1\right) M$$

 L2 := RamRatPolyToLMOp(x^(1/2)/(1 + x), x, M, b);

$$\left(-x^{2}+x -1\right) M +x$$

 L := LMOpLCLM([L1, L2], x, M, b);

$$\left(4 x^{18}+2 x^{16}-2 x^{15}-2 x^{13}+2 x^{12}+2 x^{10}+4 x^{9}+x^{7}-x^{6}-x^{4}+x^{3}+x +1\right) M^{2}+\left(-4 x^{15}-4 x^{14}-4 x^{13}-6 x^{12}-2 x^{11}-4 x^{10}-2 x^{9}-2 x^{8}-x^{7}-4 x^{6}-3 x^{5}-4 x^{4}-4 x^{3}-x^{2}-2 x -1\right) M +4 x^{13}+6 x^{12}+6 x^{11}+2 x^{10}+2 x^{7}+3 x^{6}+3 x^{5}+3 x^{4}+3 x^{3}+3 x^{2}+x$$

 L := collect(L, M, factor);

$$\left(2 x^{3}+x +1\right) \left(2 x^{9}+1\right) \left(x^{6}-x^{3}+1\right) M^{2}-\left(x^{2}+1\right) \left(2 x^{3}+1\right) \left(2 x^{10}+2 x^{9}+x^{5}+2 x +1\right) M +x \left(1+2 x \right) \left(2 x^{9}+x^{3}+1\right) \left(x^{2}+x +1\right)$$

 r, d := degree(L, M), degree(L, x);

$$2, 18$$

We use each of the available methods in turn.

 RiccatiSolve(L, x, M, b, 'free' = K, 'method' = 'basic_petkovsek');

$$\left\{\frac{\left(2 K_{2} x^{\frac{9}{2}}+2 x^{3} K_{1}+K_{2} x^{\frac{3}{2}}+2 K_{1}\right) \left(1+2 x \right)}{\left(x^{2}-x +1\right) \left(2 x^{3}+1\right) \left(2 K_{2} x^{\frac{3}{2}}+2 K_{1} x +K_{2} \sqrt{x}+2 K_{1}\right)}\right\}$$

 RiccatiSolve(L, x, M, b, 'free' = K, 'method' = 'improved_petkovsek');

$$\left\{\frac{\left(2 x^{\frac{9}{2}} K_{1}+K_{2} x^{3}+x^{\frac{3}{2}} K_{1}+K_{2}\right) \left(1+2 x \right)}{2 \left(x^{2}-x +1\right) \left(2 x^{\frac{3}{2}} K_{1}+\sqrt{x}\, K_{1}+K_{2} \left(1+x \right)\right) \left(x^{3}+\frac{1}{2}\right)}\right\}$$

 RiccatiSolve(L, x, M, b, 'free' = K, 'method' = 'hermite_pade_approximants');

$$\left\{\frac{\left(1+2 x \right) \left(2 K_{2} x^{\frac{9}{2}}+x^{3} K_{1}+K_{2} x^{\frac{3}{2}}+K_{1}\right)}{\left(2 x^{3}+1\right) \left(x^{2}-x +1\right) \left(2 K_{2} x^{\frac{3}{2}}+K_{1} x +K_{2} \sqrt{x}+K_{1}\right)}\right\}$$

The computation with the basic method à la Petkovšek is a little bit slow, perhaps $10$ s, because there are $16 \times 8 = 128$ pairs $(A, B)$ of divisors respectively of the coefficients $\ell_0$ and $\ell_r$ to consider. The improved method runs slightly better, $5$ s say, and the Hermite-Padé approach is seriously better on this example.

The result shows that the set of solutions is the image of a projective line.