    # Computing Hypergeometric Solutions of Second Order Linear Differential Equations using Quotients of Formal Solutions and Integral Bases

We present two algorithms for computing hypergeometric solutions of second order linear differential operators with rational function coefficients. Our first algorithm searches for solutions of the form (∫ r dx)·_2F_1(a_1,a_2;b_1;f) where r,f ∈Q(x), and a_1,a_2,b_1 ∈Q. It uses modular reduction and Hensel lifting. Our second algorithm tries to find solutions in the form (∫ r dx)·( r_0 ·_2F_1(a_1,a_2;b_1;f) + r_1 ·_2F_1'(a_1,a_2;b_1;f) ) where r_0, r_1 ∈Q(x), as follows: It tries to transform the input equation to another equation with solutions of the first type, and then uses the first algorithm.

## Authors

04/02/2020

### Automated solving of constant-coefficients second-order linear PDEs using Fourier analysis

After discussing the limitations of current commercial software, we prov...
07/07/2017

### Measured Multiseries and Integration

A paper by Bruno Salvy and the author introduced measured multiseries an...
11/30/2021

### Numerical solution of several second-order ordinary differential equations containing logistic maps as nonlinear coefficients

This work is devoted to find the numerical solutions of several one dime...
07/27/2017

### Dealing with Rational Second Order Ordinary Differential Equations where both Darboux and Lie Find It Difficult: The S-function Method

Here we present a new approach to search for first order invariants (fir...
01/25/2019

### Symbolic integration of hyperexponential 1-forms

Let H be a hyperexponential function in n variables x=(x_1,…,x_n) with c...
01/23/2020

### Canonical form of modular hyperbolas with an application to integer factorization

For a composite n and an odd c with c not dividing n, the number of solu...
10/17/2020

### Learning second order coupled differential equations that are subject to non-conservative forces

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

A second order homogeneous linear differential equation with rational function coefficients

 A2y′′+A1y′+A0y=0 (3)

corresponds to the differential operator

 L=A2∂2+A1∂+A0∈Q(x)[∂]

where . Another representation of (3) is

 L(y)=0.

This paper gives two heuristic (see Remarks

1.2 and 3.4) algorithms to find a hypergeometric solution of (3) in the form (1) and (2). The form (2) is more general than in prior works Fang and van Hoeij (2011); Kunwar and van Hoeij (2013); van Hoeij and Vidunas (2015); Kunwar (2014). Papers van Hoeij and Vidunas (2015) and Kunwar (2014) were restricted to a specific number of singularities (4 in van Hoeij and Vidunas (2015) and 5 in Kunwar (2014)). Papers Kunwar and van Hoeij (2013) and Fang and van Hoeij (2011) were restricted to specific degrees (degree 3 in Kunwar and van Hoeij (2013) and a degree-2 decomposition in Fang and van Hoeij (2011)). Our algorithms are not restricted to a specific number of singularities or a specific degree. Moreover, our algorithms can find algebraic functions in (1) and (2).

Our first algorithm, Algorithm 3.1, tries to find solutions of (3) in the form of (1). Our second algorithm, Algorithm 4.5, tries to reduce equations with solutions of form (2) to equations and then calls our first algorithm.

We assume that (3) has no Liouvillian solutions and hence irreducible. Otherwise one can solve (3) with Kovacic’s algorithm Kovacic (1986).

Let ( input) be a second order linear differential operator, regular singular (details in Section 2.1) and without Liouvillian solutions:

• Goal 1: (Algorithm 3.1) Find a solution of form (1) if it exists.

• Goal 2: (Algorithm 4.5) Try to transform to a simpler operator (which hopefully has a solution in form (1)).

The crucial steps for Goal 1 are to find (candidates for) , , and the pullback function . Finding the parameters , , is the combinatorial part; Theorem 3.2 helps us to eliminate the vast majority of cases. Given (or equivalently a base operator ), if we know the value of a certain constant , by comparing quotients of formal solutions of and , we can compute . We have no direct formula for ; to obtain it with a finite computation, we take a prime number . Then, for each we try to compute modulo . If this succeeds, then we lift modulo a power of , and try rational number reconstruction.

Goal 2 is to find a transformation to convert to a simpler operator . The key idea is to follow the strategy of the POLRED algorithm in Cohen and Diaz Y Diaz (1991). It takes as input a polynomial and finds an element whose minimal polynomial is close to optimal. It works as follows:

1. “Finite points”: Compute an integral basis.

2. “Valuations at infinity”: Find an integral element with (near) optimal absolute values.

Our key idea is to apply POLRED’s strategy to . First compute an integral basis (introduced in Kauers and Koutschan (2015)). Then normalize at infinity (following (Trager, 1984, Section 2.3) where this is done for function fields) so we can select an element with minimal valuations at infinity.

###### Example 1.1 (Rational Pullback Function).

The differential operator

 Linp=147x(x−1)(x+1)∂2+(266x2−42x−98)∂+20x−5

has a -type solution in the form of (1), which is

 Y(x)=exp(∫rdx)⋅2F1(542,1142;23;f)

where

 exp(∫rdx)=(x+1)−521andf=4x(x+1)2 (4)

Section 3.4 shows how to find the parameters . Then is computed with the quotient method:

###### Remark 1.1 (The Quotient Method).

The hypergeometric function

 2F1(542,1142;23;x)

is a solution of the Gauss hypergeometric differential operator

 LB=∂2+(29x−14)21x(x−1)∂+551764x(x−1).

has two solutions at . They are

 y1(x) =2F1(541,1142;23,x)=1+551176x+…, y2(x) =x13(1+4752352x+194132519361664x2+…).

The so-called exponents of at are the exponents of in the dominant terms of the solutions and , so the exponents are and . The minimal operator for has the following solutions at :

 y1(f) =1+55294x−493986436x2+16135823304946208x3+…, y2(f) =c⋅x13(1+83588x+68051210104x2+…)

for some constant that depends on . The exponents are again and , because is a root of with multiplicity (Theorem 2.1). Let

 Y1(x) =exp(∫rdx)y1(f)=1−598x+4399604x2+…, (5) Y2(x) =exp(∫rdx)y2(f)=c⋅x13(1−19196x+…). (6)

(5) and (6) form a basis of solutions of . Here is the same as in (4). Denote the quotients of the formal solutions of and by

 q =y1(x)y2(x) Q =Y1(x)Y2(x)=y1(f)y2(f)=q(f)

respectively. It follows that gives a series expansion of at . Given enough terms we can compute with rational function reconstruction. This Quotient Method was already used in (van Hoeij and Vidunas, 2015, Section 5.1). In order to turn this into an algorithm for solving differential equations we need to answer the following questions:

1. [label=Q0.]

2. How many terms are needed to reconstruct ? This is equivalent to finding a degree bound for .

3. How to find the parameters , , ? This is the combinatorial part of our algorithm.

4. The exponents , of at only determine up to a constant factor (see Remark 2.2 in Section 2.3). This means the quotient is only known up to a constant . How to find this constant?

5. What if has logarithmic solutions at ?

6. What if is an algebraic function?

7. What if does not have solutions in the form of (1), but has solutions in the form of (2)?

Remark 3.2, Section 3.4, Section 3.6, and Section 3.5.2 provide answers to Q1, Q2, Q3, and Q4 respectively. Section 3.6 answers Q5. The parts Q1,…,Q5 are already in our ISSAC 2015 paper Imamoglu and van Hoeij (2015). Algorithm 4.5, which finds solutions of form (2), is new compared to Imamoglu and van Hoeij (2015). So Q6 is the main new part in this paper. It will be discussed in Section 4. Example 1.3 will illustrate to Q6.

###### Remark 1.2.

Both algorithms are very effective in practice but they are not proven. For completeness for Goal 1 we still need a theorem for good prime numbers. A good prime is a prime for which reconstruction will work.

###### Example 1.2 (Algebraic Pullback Function).

The differential operator

 Linp=∂2+14x4−44x3+1206x2−44x+1(x2−34x+1)2x2

has a -type solution in the form of (1), which is

 Y(x)=exp(−12∫rdx)⋅2F1(13,23;1;f)

where

 −x5+22x4−55x3−343x2+58x−1+6x(x2−7x+1)√x2−34x+1x(x4−41x3+240x2−41x+1)(x+1)

and

 f=121+30x−24x2+x3−(x2−7x+1)√x2−34x+11+3x+3x2+x3.

Here the pullback function is an algebraic function: is an algebraic extension of of degree ( is if and only if is a rational function, as in Example 1.1). Algorithm 3.1 can find this solution.

###### Example 1.3 (Finding Solutions in the form of (2) using an Integral Basis).

Consider the differential operator222Prof. Jean-Marie Maillard sent us this differential operator.

 Linp=∂2 −512x5+384x4−64x3−88x2−10x−1x(4x−1)(4x+1)(16x3+24x2+5x+1)∂ +512x5+64x4−128x3−60x2−8x−1x2(4x−1)(4x+1)(16x3+24x2+5x+1).

Algorithm 3.1 can not solve . We try to transform to simpler operator . First we compute an integral basis. Then we normalize the basis at infinity and obtain where

 B0= 16x4−x2(16x3+24x2+5x+1)x∂ +−34359738400x3−51539607556x2−10737418241x−2147483648(16x3+24x2+5x+1)x

and

 B1=16x3−x(16x3+24x2+5x+1)x∂+−32x2−4x−1(16x3+24x2+5x+1)x.

We try to find a suitable . It should be a combination of and . For this example, we take . This is called a gauge transformation. It maps solutions of to solutions of

 ~Linp=∂2 +48x2−1x(16x2−1)∂+1616x2−1.

has a solution in the form of (1),

 y(x)=2F1(12,12;1;16x2)

which is easy to find with Algorithm 3.1. Then we apply the inverse gauge transformation and obtain a solution of in the form of (2), which is

## 2 Preliminaries

This section recalls the concepts needed in later sections.

### 2.1 Differential Operators, Singularities, Formal Solutions

We start with some classical definitions which can be also found in Ince (1926); van der Put and Singer (2003); Kunwar (2014); Fang (2012); Debeerst (2007); Yuan (2012).

###### Definition 2.1.

Let be an operator of order .

1. [label=()]

2. A point is called a singularity of if it is a zero of the leading coefficient of or a pole of any other coefficient of . The point is called a singularity if is a singularity of . Here is the differential operator obtained from via a change of variables (note that sends to .

3. If is not a singularity, then it is called a regular point of .

4. A singularity is called a regular singularity if is analytic at for . The point is a regular singularity if is a regular singularity of .

5. is regular singular if all its singularities are regular singular.

### 2.2 Gauss Hypergeometric Differential Operator and 2f1 Function

Definitions also can be found in Wang and Guo (1989); Kunwar (2014); Fang (2012). Let , , . The operator

 LB=x(1−x)∂2+(b1−(a1+a2+1)x)∂−a1a2

is called Gauss hypergeometric differential operator (GHDO). The solution space of in a universal extension Fang (2012) has dimension 2 because the order of is 2. One of the solutions of at is the Gauss hypergeometric function. It is denoted by and defined by the Gauss hypergeometric series

 2F1(a1,a2;b1;x)=∞∑k=0(a1)k(a2)k(b1)kk!xk.

Here denotes the Pochammer symbol. It is defined as and .

has three regular singularities at the points , , and , with exponents , , and respectively. We denote the exponent differences of a GHDO as , , . We may assume , otherwise is reducible (it has exponential solutions).

Let be if , and the denominator of if . We will only consider for which has no Liouvillian solutions. From the Schwarz list Schwarz (1873) one finds that this is equivalent to .

### 2.3 Transformations and Singularities

We summarize properties of transformations in this section. These properties can also be found in Kunwar (2014); Fang (2012); Debeerst (2007); Yuan (2012).

Let be a differential operator of order , and let be a solution of . We consider the following transformations that send solutions of to solutions of another second order differential operator .

1. Change of variables:
, where .
For this means substituting .
Notation: .

2. Exp-product:
, where .
For this means .
Notation: .

3. Gauge transformation:
, where .
For this means computing the least common left multiple of and , and right-dividing it by .
Notation: .

###### Remark 2.1.

Transformations can affect singularities and exponents.

1. [label=()]

2. If a transformation can send a singular point to a regular point , then we call a false singularity.

3. A singularity is a false singularity Debeerst (2007) if and only if is not logarithmic and the exponent difference is 1.

4. If is a singularity of and if transformation can send to an equation for which all solutions of are analytic at , then we call a removable singularity.

5. A point is removable Debeerst (2007) if and only if is not logarithmic and the exponent difference is an integer. Non-removable singularities are called true singularities.

6. A point is a true singularity if and only if the exponent difference is not an integer or is logarithmic.

7. If , then and are called gauge equivalent. If and are the solution spaces of and respectively, then maps to , i.e., .

###### Remark 2.2.

The quotient method (Remark 1.1 in Section 1) can only use true singularities, otherwise, , the quotients of solutions of , would only be known up to a Möbius transformation instead of up to a constant.

###### Remark 2.3.

At the moment, Algorithms

3.1 and 4.5 are only implemented for rational function coefficients. However, if are algebraic, then the three transformations may turn an operator with rational function coefficients into an operator with algebraic function coefficients.

###### Theorem 2.1.

Bostan et al. (2011) Let the GHDO have exponent differences at , at , and at . Let . If , then has the following exponent difference at :

• if has a zero at with multiplicity ,

• if has a zero at with multiplicity ,

• if has a pole at with order .

## 3 Computing Solutions of a Second Order Linear Differential Operator in the form of (1) by using Quotients of Formal Solutions

This section gives our first algorithm, which looks for solutions of the form of (1).

### 3.1 Problem Statement

Given a second order linear differential operator , irreducible and regular singular, we want to find a -type solution of the differential equation of the form of (1). This is equivalent to finding transformations 1 and 2 from a GHDO to . Therefore, we need to find

1. (i.e., find ),

2. parameters and of the change of variables and exp-product transformations such that .

###### Algorithm 3.1.

General Outline of find_2f1.

1. INPUT: and (optional) where

1. a second order regular singular irreducible operator,

2. bound for the algebraic degree (See Example 1.2). If omitted, then which means our implementation tries and .

2. OUTPUT: Solutions of in the form of (1), or an empty list.

3. For each :

4. Use Section 3.4 to compute candidates for and . This is the combinatorial part of the algorithm.

5. For a candidate , compute formal solutions of and at a non-removable singularity (see Remark 2.2 in Section 2.3) up to precision . Take the quotients of formal solutions and compute series expansions for and which will be used to compute

 f=q−1(cQ(x)) (7)

in the next step.

6. Choose a good prime number and try to find mod by looping as in Section 3.5. For each :

1. Compute mod from equation (7) and use it to reconstruct mod (the image of in ). If it fails for every , then proceed with the next candidate GHDO (if any) in Step 2. If no candidates remain, then return an empty list.

2. If rational reconstruction in Step 3.1 succeeds for some values, then apply Hensel lifting (Section 3.6) to find mod a power of . Then try rational number reconstruction. If it does not fail for at least one value, then we have . If no solution is found (see Remark 3.4 in Section 3.6.1), then proceed with the next candidate GHDO (if any) in Step 2. If no candidates remain, then return an empty list.

3. Use Section 2.3 to compute the parameter of the exp-product transformation.

7. Return a basis of -type solutions of .

Step 2 is explained in Sections 3.2 and 3.4. Step 3 is the quotient method, see Section 3.5 for more. Steps 3.2 and 3.3 are explained in Sections 3.6 and 3.7 respectively. A Maple implementation of Algorithm 3.1 and some examples can be found at Imamoglu (2015a).

### 3.2 Degree Bounds for Pullback Functions

###### Theorem 3.1 (Riemann-Hurwitz Formula).

Let and be two algebraic curves with genera and respectively. If is a non-constant morphism, then

 2gX−2=deg(f)(2gY−2)+∑p∈X(ep−1). (8)

Here denotes the ramification order at . See Hartshorne (1977) for more details.

Let be a GHDO, , and assume that

 LBf:P1⟼P1−−−−−−→Cr→ELinp.

Section 3.1 of Imamoglu and van Hoeij (2015) gives an a priori bound for ,

 df≤{6(ntrue−2),logarithmic case,36(ntrue−73),non-logarithmic case. (9)

where is the number of true singularities of . Algorithm 3.1 uses this only as an initial degree bound.

### 3.3 Riemann-Hurwitz Type Formula For Differential Equations

###### Remark 3.1.

Let be any algebraic curve and be its function field. The ring is the ring of differential operators on . Here . An element is a differential operator defined on the algebraic curve .

###### Theorem 3.2.

(Baldassari and Dwork, 1979, Lemma 1.5) Let , be two algebraic curves with genera , , and function fields , respectively. Let be a non-constant morphism. The morphism corresponds to a homomorphism , which in turn corresponds to a homomorphism . If with and is the corresponding element in , then

 Covol(L2,X)=deg(f)⋅Covol(L1,Y) (10)

where

 Covol(L,X):=2gX−2+∑p∈X(1−Δ(L,p))

and where is the absolute value of the exponent difference of at .

###### Proof..

Following Baldassari and Dwork (1979), take finite sets and in such a way that all singularities of are in , all singularities of are in , and all branching points in are in as well.

 #T=∑p∈T1 =∑p∈Tep+∑p∈T(1−ep) (11) =deg(f)⋅#S+∑p∈X(1−ep) (12) =deg(f)⋅#S−(2gX−2−deg(f)(2gY−2)). (13)

From (12) to (13) we used (8). Then,

 ∑p∈X(1−Δ(L2,p)) =∑p∈T(1−Δ(L2,p)) (14) =∑p∈T1 − ∑p∈TΔ(L2,p) (15) =#T − deg(f)∑s∈SΔ(L1,s). (16)

Then, combine (13) and (16) and get

 2gX−2+∑p∈X(1−Δ(L2,p))=deg(f)(2gY−2+∑s∈Y(1−Δ(L1,s))) (17)

which is the same as (10). ∎

###### Corollary 3.1.

Let and suppose that where is a GHDO with exponent differences at . Since an exp-product transformation does not affect exponent differences, Theorem 3.2 gives the following equation for :

 −2+∑p∈P1(1−Δ(Linp,p))=deg(f)⎛⎝−2+∑i∈{0,1,∞}(1−αi)⎞⎠. (18)
###### Corollary 3.2.

Let and be as in Corollary 3.1. Both have rational function coefficients. This time, suppose that in are algebraic functions. Then for an algebraic curve whose function field is an algebraic extension of both and . Let and denote the degrees of these extensions.

Applying (10) to both field extensions gives:

 Covol(Linp,P1)=dfaf⎛⎝−2+∑i∈{0,1,∞}(1−αi)⎞⎠. (19)

### 3.4 Candidate Exponent Differences

This section explains how to obtain exponent differences for candidate GHDOs.

###### Algorithm 3.2.

General Outline of find_expdiffs.

1. INPUT: , , and where

1. [label=()]

2. the list of exponent differences of at its true singularities,

3. the (possibly empty) list of exponent differences of at its removable singularities,

4. = candidate algebraic degree.

2. OUTPUT: A list of all lists of integers or rational numbers where is a list of candidate exponent differences and is a candidate degree for such that:

1. [label=()]

2. For every exponent difference in there exists with such that for some .

3. The multiplicities are consistent with (8), and their sums are compatible with , see the last paragraph in Step .

3. Let . After reordering we may assume that , , and for . For each we use CoverLogs in Imamoglu (2015a) to compute candidates for .

Algorithm CoverLogs computes candidates that meet these requirements:

• Logarithmic singularities are true singularities with integer exponent differences. If has at least one logarithmic singularity with exponent difference , then a candidate must have at least one logarithmic singularity; at least one of the must be an integer that divides , and for every there must be at least one such that divides .

• for some      .

• Theorem 2.1.

If , then algorithm CoverLogs also computes the exact degree of using Theorem 2.1 which shows that must be the sum of the logarithmic exponent differences of . Otherwise, it uses (9) to compute a bound for , and uses it as to compute a candidate degree.

4. We will explain only the case , and only , which is the case , where and .

5. Let . Let be one of the candidates from algorithm CoverLogs. We need to find candidates for and .

6. The logarithmic singularities of come from the point . Non-integer exponent differences of must be multiples of or . Let be the set of non-logarithmic exponent differences of and be the set of exponent differences of at its removable singularities. Consider the set

 Γ1={ΓA={max(SN)b:b=1,…,df}if SN≠∅,ΓB={ab:a∈SR∪{1},b=1,…,df}otherwise% .

(or , but if so, we may interchange them) must be one of the elements of . We loop over all elements of . Assume that a candidate for is chosen. Let . Now consider the set

 Γ∞={ΓA∪ΓBif Ω=∅,{gb:g=gcd(Ω):b=1,…,df}otherwise.
7. Now take all pairs satisfying (19), , , with additional restrictions on , as follows:

8. For every potential non-zero value for one of the ’s we pre-compute a list of integers by dividing all exponent differences of by and then selecting the quotients that are integers. Next, let be the set of all that can be written as the sum of a sublist of . Each time a non-zero value is taken for one of the , it imposes the restriction . This means that we need not run a loop for , instead, we run a (generally much shorter) loop for (taking values in the intersection of the ’s so far) and then for each such compute from (19). We also check if .

9. Return the list of candidate exponent differences with a candidate degree, the list of lists , for candidate GHDOs.

### 3.5 Quotient Method

In this section, we explain a method to recover the pullback function . We will explain our algorithm for rational pullback functions. For algebraic pullback functions, the only difference is the lifting algorithm, which is explained in Section 3.6. Note that we can always compute the formal solutions of a given differential equation up to a finite precision.

#### 3.5.1 Non-logarithmic Case

Let the second order differential equation be given. Let be a GHDO such that Let and . If is a singularity of and is a singularity of , then we say that comes from when .

After a change of variables we can assume that is a singularity of that comes from the singularity of . This means and we can write where , is the multiplicity of , and the dots refer to an element in .

Let and be the formal solutions of at . The following diagram shows the effects of the change of variables and exp-product transformations on the formal solutions of ,

 yi(x) f→Cyi(f)r→EYi(x)=exp(∫rdx)yi(f),i∈{1,2}

where and are solutions of .

Let be a quotient of formal solutions of . The change of variables transformation sends to