A Simple Method for Computing Some Pseudo-Elliptic Integrals in Terms of Elementary Functions
aa r X i v : . [ c s . S C ] A p r A Simple Method for Computing Some Pseudo-EllipticIntegrals in Terms of Elementary Functions
Sam Blake
The University of Melbourne
DRAFT version: 22/04/2020:00:40:20
Abstract
We introduce a method for computing some pseudo-elliptic integrals in terms of ele-mentary functions. The method is simple and fast in comparison to the algebraic case ofthe Risch-Trager-Bronstein algorithm[13][15][1]. This method can quickly solve manypseudo-elliptic integrals, which other well-known computer algebra systems (CAS) ei-ther fail, return an answer in terms of special functions, or require more than 20seconds of computing time. Unlike the symbolic integration algorithms of Risch[13],Davenport[3], Trager[15], Bronstein[1] and Miller[9]; our method is not a decision pro-cess. The implementation of this method is less than 200 lines of Mathematica codeand can be easily ported to other CAS that can solve systems of linear equations.
The problem of finding elementary solutions to integrals of algebraic functions has chal-lenged mathematicians for centuries. In 1905, Hardy conjectured that the problem maybe unsolvable[7]. Sophisticated algorithms have been developed, including the famed Rischalgorithm[13] and its modern variants by Davenport[3], Trager[15], Bronstein[1] and Miller[9].However, it is known that the implementation of these algorithms are highly complex andsometimes fail, are incomplete or hang[5].We will be describing a seemingly new method for computing pseudo-elliptic integrals.
Definition 1.
For our purposes, a pseudo-elliptic integral is of the form Z a ( x ) b ( x ) p ( x ) n/m dx, where a ( x ) , b ( x ) , p ( x ) are polynomials, deg x ( p ( x )) >
2, gcd( n, m ) = 1 . derivative divides method is a substitution method that finds all composite functions, u = g ( x ), of the integrand f ( x ), and tests if f ( x ) divided by the derivative of u is independentof x after the substitution of u = g ( x ). In other words, up to a constant factor the derivativedivides method simplifies integrals of the form R f ( g ( x )) g ′ ( x ) dx to R f ( u ) du . The followingintegral illustrates the derivative divides method Z x √ x dx = Z √ u du, where u = 1 + x and du = 3 x dx . This method was first implemented in Moses sym-bolic integrator, SIN, in 1967[10] and is used in some CAS prior to calling more advancedalgorithms[6, pp. 473-474].A slightly more difficult example where the derivative divides method fails is Z q x − √ x − dx. In this case we make the substitution u = √ x −
1, then 2 udu = dx . Furthermore, we needto express x − u , which is u = x −
1. Then the integral becomes Z u √ u + u du. While this integral was somewhat more difficult than the previous example, we could stillpick a composite function, √ x −
1, of our integrand, p x − √ x − , to use as our u substitution. In the following well-known example, this is not immediately possible Z x − x + 1) √ x dx. A common approach to solve this integral is rearranging the integrand into Z x (1 − x ) x ( x + 1 /x ) p x ( x + 1 / x ) dx = Z − x ( x + 1 /x ) p ( x + 1 /x ) − dx. Then the substitution u = x + 1 /x , du = (1 − x ) dx yields the integral Z duu √ u − , which can be transformed into a rational function using the Euler substitution[4].Moses devised methods for integrating specific classes of integrals, these included Z R ( x ) exp( p ( x )) dx, where R ( x ) is a rational function of x and p ( x ) is a polynomial in x [10, pp. 85]. Shortly after,these methods were phased-out with the work of Risch. However in most CAS the algebraiccase of the Risch algorithm is either partially implemented, not implemented, or containscomputational bottlenecks that result in long computations. So for algebraic functions adomain-specific approach still has merit. We require that the integrand is of the form p ( x ) q ( x ) r ( x ) n/m , where p ( x ) , q ( x ) , r ( x ) ∈ Z [ x ] andgcd( n, m ) = 1. If the integrand is not in this form then we attempt to write it in this formplus a rational function of x . For example, x + 1 + x √ x + 1( x + 1) √ x + 1 = y x + 1 + xx + 1 , where y = x + 1. The rational part, xx +1 , is integrated using known algorithms for ra-tional function integration[2]. An algorithmic way to create this representation is given inTrager[15], however it requires an integral basis and is beyond the scope of this paper. Following on from our example integral R x − x +1) √ x dx , where making the substitution u = x +1 /x resulted in the integral R duu √ u − . We would like to generalise this method, however thedifficulty was in the choice of the algebraic manipulation to the form R − / x ( x +1 /x ) √ ( x +1 /x ) − dx in order to discover a rational substitution, which simplifies the integral. Consequently ourapproach does not directly rely on such an algebraic manipulation of the integrand.Our method attempts to parameterise constants a , a , a , polynomials a ( u ) , b ( u ) and asubstitution of the form u = s ( x ) x k such that Z p ( x ) q ( x ) r ( x ) n/m dx = Z a ( u ) b ( u ) (cid:0) a u + a u + a (cid:1) n/m du, (3.1)where deg x ( r ( x )) > n, m ) = 1. Consequently, our method does not directlycompute the integral, and requires a recursive call to an algebraic integrator . We notethat a reduction to this form does not guarantee an elementary solution (for example, when m > , a = 0 an elementary form is often not possible).Our method is broken into two parts. The first part is computing the radicand part ofthe integral , which is a parameterisation of a , a , a and the u substitution. The second partis computing the rational part of the integral , which is a parameterisation of a ( u ) and b ( u ). The radicand part of the integral.
Clearly, if we cannot parameterise the radicand r ( x ) tothe form a u + a u + a for a given substitution, then we cannot find a parameterisation of(3.1). Thus, we begin by computing the radicand part of the integral, which requires solving r ( x ) = num (cid:0) a u + a u + a (cid:1) , for the constants a , a , a and the substitution u = s ( x ) /x k . We do this by iterating over0 < d ≤ N such that s ( x ) = d P i =0 c i x i , where for each d we iterate over 0 < h ≤ M such that u = s ( x ) /x h . Given a candidate u , we then iterate over the radicands r ( u ) = a u + a , r ( u ) = a u + a , r ( u ) = a u + a u + a , where for each radicand we require m | deg x (den ( r i ( u ))) (otherwise we would have a frac-tional power in the denominator of r i ( u ) n/m ), then we solve r ( x ) = num ( r i ( u )) . (3.2)If deg x ( r ( x )) = deg x (num ( r i ( u ))), then a solution does not exist (as we have already con-sidered lower-degree polynomials). Otherwise we solve (3.2) by equating coefficients of x and solving the system of equations for the unknowns a , a , a , c , c , · · · , c d . If a solution(or multiple solutions) exists we move to computing the rational part of the integral, oth-erwise we move onto the next radicand or candidate substitution. If no solution exists tothe radicand part of the integral for any candidate u -substitutions, then our method fails tocompute the integral. The rational part of the integral.
Given the substitution and solution set of the radicandpart of the integral, we now look to solve the rational part of the integral, which is given by p ( x ) q ( x ) = a ( u ) b ( u ) u ′ ( x )den ( a u + a u + a ) n/m , (3.3)where a , a , a , u ( x ) are known and a ( u ), b ( u ) are unknown. The degree bound estimateof a ( u ) and b ( u ) is given by D = deg x ( u ( x )) + deg x ( u ′ ( x )) + max (deg x ( p ( x )) , deg x ( q ( x ))). If such an integrator is not available then a reasonable implementation could call a lookup table ofalgebraic forms[11] followed by a rational function integrator[2].
We solve (3.3) by increasing the degree, d , of a ( u ) and b ( u ) from 1 to the degree bound, D ,where for each iteration we solve p ( x ) b ( u ( x )) den (cid:0) a u + a u + a (cid:1) n/m − q ( x ) a ( u ( x )) u ′ ( x ) = 0 , where a ( u ) = d P i =0 v i u i , b ( u ) = d P i =0 v d + i +1 u i , and a ( u ( x )), b ( u ( x )) are rational functions in x after replacing u with the candidate substitution. As before, we equate powers of x and solvefor the unknowns v , v , · · · , v d − . If a solution is found, then we have a complete solutionto (3.1) and we stop. Otherwise if we have iterated up to the degree bound, and iteratedthrough all solution sets from the radicand part of the integral and we have not computeda solution, then the candidate substitution is rejected and must return to the radicand partof the integral to try the next substitution. Example 3.1.
We will apply the method detailed above to compute the following integral Z ( x − √ x − x + 1( x + 1) dx. This integral is already in our canonical form.
The radicand part of the integral . We find the substitution u = ( c x + c ) /x yields asolution to the radicand part of the integral, and is parameterised as follows x − x + 1 = num (cid:0) a u + a u + a (cid:1) = a c x + a c x + a x + 2 a c c x + a c x + a c . Then, equating coefficients of powers of x yields the system of equations a c = 1 a c = − a c c = 1 a = 0 a c = 0 a c = 0 , which has no solution, so we use a linear radicand in u as follows x − x + 1 = num ( a u + a ) = a c x + a x + a c . Again, we equate coefficients of powers of x and we have the system of equations a c = 1 a = − a c = 1 , which has the solution a = − , a = 1 , c = 1 , c = 1. Thus, the radicand part of theintegral is u −
1, where u = (1 + x ) /x . The rational part of the integral . Now we see if a solution exists to the rational part of theintegral. The degree bound on the solution to the rational part is 3. When the degree is 1,we have no solution. When the degree is 2, we have a ( u ) b ( u ) = v u + v u + v v u + v u + v . For the rational part, we are solving the following equation( x − x + 1) = (cid:18) v u + v u + v v u + v u + v (cid:19) den( u − − / u ′ ( x ) , where den( u − − / = 1 /x . After replacing u with (1 + x ) /x and u ′ ( x ) with ( x − /x we have x − x + 1) = ( x −
2) ( u v + uv + v ) x ( u v + uv + v ) = ( x −
2) ( x v + x v + x v + v + 2 x v + x v ) x ( x v + x v + x v + v + 2 x v + x v ) , which after clearing denominators is a polynomial equation in x , given by − v x − v x + ( v − v ) x + ( v − v ) x + ( v − v ) x + (2 v − v ) x +(3 v − v ) x +(3 v − v ) x +(8 v − v ) x +5 v x +(2 v − v ) x +7 v x +2 v x +2 v = 0 , which we solve for the undetermined coefficients v , v , v , v , v , v . Then equating coeffi-cients of powers of x yields the system of equations2 v = 02 v = 07 v = 02 v − v = 05 v = 08 v − v = 03 v − v = 03 v − v = 02 v − v = 0 − v + v = 0 − v + v = 0 − v + v = 0 − v = 0 − v = 0 , which has the solution v = v , v = 0 , v = 0 , v = 0 , v = 0. Thus, the rational part ofthe integral is v v u = 1 u and the integral is given by Z ( x − √ x − x + 1( x + 1) dx = Z √ u − u du = − √ u − u + tan − (cid:0) √ u − (cid:1) = − x √ x − x + 1 x + 1 + tan − (cid:18) √ x − x + 1 x (cid:19) . Our implementation in Mathematica took 0.085 seconds to compute this integral.
In this section we will look at two generalisations of this method to extend the classes ofintegrals with elementary solutions. We also consider solutions in terms of special functions.
This method can be generalised to include reductions of the form Z p ( x ) q ( x ) r ( x ) n/m dx = Z a ( u ) b ( u ) (cid:0) a u k + a u k + a (cid:1) n/m du and Z p ( x ) q ( x ) r ( x ) n/m dx = Z a ( u ) b ( u ) (cid:0) a u k + a (cid:1) n/m du. For both these forms we use essentially the same procedure as in the linear or quadraticreduction. However, in terms of efficiency, obviously the more forms we try to solve for theradicand part of the integral, the longer the process takes to iterate through these possibleforms.
If the integration of the reduced form is permitted to contain special functions and therecursive call to the integrator has the capability to integrate algebraic functions in terms ofspecial functions, then we can obtain integrals to a much larger class of algebraic functions.
Example 4.1.
For the following integral, our method makes the substitution u = (2 x − /x ,to obtain the reduction Z x + 1(2 x − √ x + 2 x − x − x − x + 1 dx = Z du u √ u + u − . This integral does not have a solution in terms of elementary functions, however Mathematicacan express this integral in terms of the Appell hypergeometric function, F Z du u √ u + u − − s − √ u + 1 s √ u + 1 √ u + u − × F
12 ; 14 ,
14 ; 32 ; − √ u , − − √ u ! . Then the integral is given by Z x + 1(2 x − √ x + 2 x − x − x − x + 1 dx = − x s x − √ x + x − x − s x + √ x + x − x − √ √ x + 2 x − x − x − x + 1 × F
12 ; 14 ,
14 ; 32 ; (cid:16) − − √ (cid:17) x x − , (cid:16) − + √ (cid:17) x x − . In the following example we recursively call our method to obtain an elementary form to adifficult integral.
Example 4.2.
For the following integral, our method makes the substitution u = √ x − /x ,then for the recursive integration, we make the substitution u = t to obtain Z x + 1(2 x − √ x − x − x + 2 dx = Z du u √ u − Z dt t √ t − − √ − x − √ √ x − x − x + 2 x ! +12 √ − x + √ √ x − x − x + 2 x ! − √ − (cid:18) x + √ x − x − x + 2 √ x √ x − x − x + 2 (cid:19) . In the above example, the integral Z du u √ u − , does not return an elementary form in either Mathematica (12.1) or Maple (2018.1). We will compare our method with the Mathematica (12.1.0), Maple (2018.1), REDUCE(5286, 1-Mar-20), and FriCAS (version 1.2.6) computer algebra systems. We will also includein the comparison a table lookup package, RUBI[12], which has been ported to a number ofcomputer algebra systems and compares favourably with most built-in integrators on a largesuite of problems[14]. We have also included an experimental algebraic integration packagedeveloped in Mathematica by Manuel Kauers[8]. Within the Kauers package we have re-placed the calls to Singular in favour of Mathematica’s built-in Groebner basis routine.We have included results from Maple twice. Once with a call of int(integrand, x) andonce with int(convert(integrand, RootOf),x) . This is because the default behaviourof Maple is to not use the Risch-Trager-Bronstein integration algorithm[15][1] for algebraicfunctions unless the radicals in the integrand are converted to the Maple
RootOf notation.Within REDUCE we have used the algint package by James Davenport[3].Our test suite is 190 integrals that can be found on github[16]. All the integrals in thesuite have a solution in terms of elementary functions.We will show the results from all the systems and packages for one integral from the testsuite. It is intriguing to see the variety of forms for this integral.Our method returns:
Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx = − √ − √ x √ x + 1 ! − √ − √ x √ x + 1 ! Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx =18 √ (cid:18) x + 1 (cid:16) x + 4 x + √ (cid:0) x + 4 x + 1 (cid:1) − p x + 1 (cid:16) / (cid:0) x + 2 x (cid:1) + 4 √ x (cid:17)(cid:17)(cid:19) − √ (cid:18) − x + 1 (cid:16) x + 4 x + √ (cid:0) x + 4 x + 1 (cid:1) + p x + 1 (cid:16) / (cid:0) x + 2 x (cid:1) + 4 √ x (cid:17)(cid:17)(cid:19) +12 √ − − x − x + √ (cid:0) x + 4 x + 1 (cid:1) √ − x −
1) + √ x + 1 (cid:16) / (2 x + 2 x ) − √ x (cid:17) Kauer’s algorithm returns:
Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx = X α − α log (cid:16) α p x + 1 − x (cid:17) Maple (default) returns:
Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx = 12 √ − √ x + 1 √ x ! − √ √ x +1 √ x + √ √ x +1 √ x − √ Maple (RootOf) returns:
Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx = 14 √ × / x − √ x + 1 x + 4 √ x + 2 × / − x + 2 √ x − ! + i √ × / ix − √ x + 1 x − i √ x + 2 × / i x + 2 √ x + 2 ! Mathematica returns:
Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx = 12 √− (cid:0) − F (cid:0) i sinh − (cid:0) √− x (cid:1)(cid:12)(cid:12) − (cid:1) +Π (cid:0) − √− i sinh − (cid:0) √− x (cid:1)(cid:12)(cid:12) − (cid:1) + Π (cid:0) √− i sinh − (cid:0) √− x (cid:1)(cid:12)(cid:12) − (cid:1) +Π (cid:16) − ( − / ; i sinh − (cid:0) √− x (cid:1)(cid:12)(cid:12) − (cid:17) + Π (cid:16) ( − / ; i sinh − (cid:0) √− x (cid:1)(cid:12)(cid:12) − (cid:17)(cid:17) where F is the elliptic integral of the first kind, and Π is the incomplete elliptic integral.REDUCE (using algint package) returns: Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx = Z x √ x + 1 x + 1 dx − Z √ x + 1 x + 1 dx : Z (cid:0) x − (cid:1) √ x + 1 x + 1 dx = − tan − (cid:18) √ x √ x +1 (cid:19) √ − tanh − (cid:18) √ x √ x +1 (cid:19) √ (cid:0) x + 1 (cid:1) q x +1( x +1) F (cid:0) − ( x ) | (cid:1) √ x + 1 + (cid:0) ( − − i ) − i √ (cid:1) (cid:0) x + 1 (cid:1) q x +1( x +1) F (cid:0) − ( x ) | (cid:1) √ x + 1 + (cid:0) √ − i ) (cid:1) i (cid:0) x + 1 (cid:1) q x +1( x +1) F (cid:0) − ( x ) | (cid:1) √ x + 1 + (cid:0) √ i ) (cid:1) i (cid:0) x + 1 (cid:1) q x +1( x +1) F (cid:0) − ( x ) | (cid:1) √ x + 1 − (cid:0) − i (cid:1) (cid:0) − / (cid:1) (cid:0) x + 1 (cid:1) q x +1( x +1) F (cid:0) − ( x ) | (cid:1) √ x + 1 where F is the elliptic integral of the first kind, and Π is the incomplete elliptic integral.The table below summarises the comparison between all systems on the test suite ofintegrals. The expression size is the leaf count( LeafCount ) in Mathematica.Table 1: A comparison of our method with major CAS and algebraic integration packages.
CAS/package Elementaryforms [%] Contains R dx [%] Containsspecialfunctions [%] Timed-out( > RootOf ) 91.6 2.0 0.0 6.3 4.32 136–238FriCAS 70.0 4.7 0.0 25.3 0.26 . –129Kauers 62.6 7.9 0.0 29.4 0.40 88–116RUBI 13.7 60.0 16.3 10.0 0.34 244–405Maple 11.6 55.3 33.2 0.0 0.34 436–943Mathematica 9.5 44.7 43.2 2.63 1.21 574–1019REDUCE (algint) 6.3 65.3 0.0 28.4 1.16 .
In regards to this table: The median time excludes integrals that timed-out. For FriCAS After posting a preprint of this comparison on the sci.math.symbolic newsgroup, Albert Rich (the creatorof RUBI) devised a general rule for integrals of this form which will be included in the next release of RUBI:
Z (cid:0) f + gx (cid:1) √ d + ex a + bx + cx = e f cd r de − be c tan − x r de − be c √ d + ex + e f cd r de − be c tanh − x r de − be c √ d + ex , when ef + dg = 0 and cd − ae = 0. RootOf notation in Maple, the time requiredincreases significantly, however the results improve significantly. We are then left wonderingwhy Maple does not make this conversion internally if the initial algebraic integration rou-tines fail?FriCAS is a fork of AXIOM. FriCAS inherits the most complete implementation of theRisch-Trager-Bronstein algorithm for integration in finite (elementary) terms. The resultssuggest that the FriCAS (version 1.2.6) integrator is either fast, hangs, or returns an error.We have reported these issues to the FriCAS developers and most issues will be fixed in thenext release[18]. Thus, the performance of FriCAS on our test suite of integrals will be muchbetter, with the exception of incomplete parts of the Risch-Bronstein-Trager algorithm foralgebraic functions, for example after 16 seconds of computing time FriCAS (version 1.2.6)returns (1) -> integrate(((4+x^6)*(-4+x^4+2*x^6)*(32-14*x^4-32*x^6-4*x^8+7*x^10+8*x^12)^(1/2))/(x^9*(-2+x^6)),x)>> Error detected within library code:integrate: implementation incomplete (residue poly has multiple non-linear factors)
When Kauer’s heuristic does not time out, then its fast, with only 11 problems tak-ing longer than 2 seconds. It also often returns a concise form. Miller[9] has extendedKauer’s heuristic to a complete algorithm for the logarithmic part of the integral of a mixedalgebraic-transcendental function. Unfortunately we did not have access to an implementa-tion for comparison.
We have shown that we can efficiently solve some pseudo-elliptic integrals. Our methodcompares favourably with major CAS and algebraic integration packages.The computational burden of our method is low, as the core computational routine re-quires solving multiple systems of linear equations, and consequently our method shouldbe tried before the more computationally expensive algorithms of Trager[15], Bronstein[1],3Kauers[8], or Miller[9].Our method, relative to the algebraic case of the Risch-Bronstein-Trager algorithm, isvery simple to implement. Our exemplar implementation in Mathematica is only a coupleof hundred lines of code and relies heavily on
SolveAlways for solving polynomial equationswith undetermined coefficients. The implementation is available on github[17].
References [1] Bronstein, M. (1990). “Integration of Elementary Functions”.
Journal of Symbolic Com-putation . 9(2), pp. 117-173.[2] Bronstein, M. (1997). “Symbolic Integration 1: Transcendental Functions”. Springer-Verlag.[3] Davenport, J. (1979). “Integration of algebraic functions”,
EUROSAM ‘79: Proceedingsof the International Symposium on Symbolic and Algebraic Computation. pp. 415-425.[4] https://en.wikipedia.org/wiki/Euler_substitution [5] http://fricas-wiki.math.uni.wroc.pl/RischImplementationStatus [6] Geddes, K. Czapor, S. Labahn, G. (1992). “Algorithms for Computer Algebra”, SpringerUS, ISBN .[7] Hardy, G. (1916).
The Integration of Functions of a Single Variable . Cambridge Uni-versity Press. Cambridge, England.[8] Kauers, M. (2008). “Integration of algebraic functions: a simple heuristic for finding thelogarithmic part”.
ISSAC ‘08 . pp. 133-140.[9] Miller, B. (2012). “On the Integration of Elementary Functions: Computing the Loga-rithmic Part”. Thesis (Ph.D.) Texas Tech University, Dept. of Mathematics and Statis-tics.[10] Moses, J. (1967). “Symbolic Integration”. MAC-TR-47, MIT, Cambridge, MA.[11] Prudnikov, A.P. Brychkov, A. Marichev, O.I. “Integrals and Series, Volume 1, Elemen-tary Functions”. CRC Press. edition 1. ISBN-10: 2881240976.[12] Rich, A. Scheibe, P. Abbasi, N. (2018). “Rule-based integration: An extensive systemof symbolic integration rules”. Journal of Open Source Software. 3.32, p1073 1-3.[13] Risch, R. (1969). “The Problems of Integration in Finite Terms”,
Trans. AMS . 139(1).pp. 167-189.4[14] https://rulebasedintegration.org/testResults [15] Trager, B. (1984). “Integration of algebraic functions”. Thesis (Ph.D.) MassachusettsInstitute of Technology, Dept. of Electrical Engineering and Computer Science.[16] https://github.com/stblake/algebraic_integration/blob/master/comparisonIntegrands.m [17] https://github.com/stblake/algebraic_integration/blob/master/AlgebraicIntegrateHeuristic.m [18][18]