Integration of a Generalized Ratio of Polynomials

This paper provides a closed-form solution to the indefinite integral of a ratio of generalized polynomials where the denominator polynomial is raised to the general order r∈Z+. Such an integral arises in physics and engineering, the solution of which allows for closed-form analysis.


Introduction
This paper provides a proof for the following theorem: Theorem 1. Given two real valued polynomials ∑ N n a n x n and ∑ M m b m x m , where the latter polynomial has roots x m which may be real or complex conjugate pairs, then a 0 + a 1 x + a 2 x 2 + · · · + a N x N (−1) k A k,m a n x n+1 x k m (n + 1) 2 F 1 k, n + 1; n + 2; where r ∈ Z + ; N and M denote the order of the numerator and denominator polynomials, respectively; x m , m ∈ [1, M] are again the roots of the denominator polynomial; 2 F 1 (. . . ) is the Gauss Hypergeometric function; and A k,m are the coefficients of the partial fraction decomposition of the denominator, which require an algorithmic approach to determine (to be explained soon).
Such an integral is useful if one needs to integrate a complicated ratio of functions, and can expand the numerator and denominator into their Taylor series representations. Integrals of this form also appear in physics and engineering, some examples being the calculation of the Petzval curvature in gradient index optics [1,2], and inverse square law theories. Furthermore, the result of this integration condenses perhaps hundreds of solutions showcased in integral tables involving ratios of various types of polynomials, into one generalized result. Such a result does not appear to be found in integral table resources [3][4][5][6][7]. References [6,7] come close, following a similar approach initially, but fall short by (i) not writing the results in a generalized closed form, and (ii) not allowing the denominator polynomial be raised to arbitrary integer exponent.
The development of the proof will be pedagogical, starting with r = 1, then r = 2, and finally general r ∈ Z + . This is to very clearly show the logical progression of starting with the simplest case, and then generalizing. Note that for simplicity, the proceeding analysis does not include denominator polynomials with degenerate roots. Such an analysis will be left for future work.

Denominator of Linear Total Exponentiation
We begin by writing down the generalized ratio of polynomials with r = 1, as an integrand: where we have written I 1 to be equal to the integral in question. The subscript 1 refers to the total exponentiation of the denominator. This notation will continue throughout this article as we analyze different values of r.
By the fundamental theorem of algebra, the denominator can be factored into a product of binomials consisting of the roots x m , m ∈ {1, 2, . . . M}.
where the integration and summation operations have been swapped due to their linearity.

Remark 1.
Note how this indicates that the roots can be complex in nature, meaning the approach at this point is valid for all polynomials. Isolating, for a moment, the term involving the denominator, one realizes that it can be decomposed via partial fraction decomposition: To find any particular A j , one simply multiplies both sides by (x − x j ) and then substitutes in x = x j -namely: where the substitution of x = x j kills all of the terms except A j . Knowing that the A coefficients are now known, we can proceed with the integration: To continue with this and the following proofs, we make use of the following Lemma: Proof. The Taylor expansion of the integrand can be found as: where we have utilized the rising factorial (sometimes referred to as the Pochhammer symbol) a (i) = (a+i−1)! (a−1)! . Integrating term-wise gives: Multiplying by 1 in a judicious manner yields: where we have used the definition of the Gauss Hypergeometric function, namely: Returning back to Equation (6), using Lemma 1 for a = 1 gives: where we have absorbed the m and n summation into the constant C.
An example of the agreement between the analytical solution and the numerical calculation is shown in Figure 1 for a denominator polynomial with one real root and two complex roots. In all numerical examples, the integration was performed via the MATLAB function cumtrapz, which utilizes the cumulative trapezoidal technique. The cumtrapz function simply has inputs of (i) the independent variable x defined by a vector and (ii) the integrand function, also defined by a vector, evaluated at the points x. After the integration, an arbitrary constant is generated and offsets the curve. Therefore, in all of the presented results the integration constants from both the numerical calculation and the closed-form solution are both set to zero.

Denominator of Squared Total Exponentiation
Next, we extend this analysis to the condition such that the denominator is squared. Much of the procedure is the same as shown above, so these steps will not be shown in such detail. We begin by rewriting the integral: Again, by the fundamental theorem of algebra, the denominator can be factored into a product of binomials consisting of the roots x m , m ∈ {1, 2, . . . M}.
Again isolating, for a moment, the term involving the denominator, one again decomposes via partial fraction decomposition: Note the extra term on the right-hand side due to the squared binomials in the denominator. We initially begin much the same way as shown in the previous section, the difference this time being that we are solving for the coefficients on the squared binomial term first by multiplying on both sides by (x − x j ) 2 and then substituting in x = x j . This yields: The same trick cannot be done to the linear binomial term; namely, multiplying on both sides by (x − x j ) and substituting x = x j . Doing so would cause a divergent term on the lefthhand side of the equation. Thus, we choose instead to multiply both sides by the entire denominator term, yielding: where we have changed the dummy index on the product to not be ambiguous. Next we isolate the term involving A 1,m : From this equation, all of the constants A 1,m can be found by realizing this expression is valid for all values of x (where x = x m so as not to cause singularities). Thus, one can choose arbitrary values of x and create the necessary number of equations to form a well-defined system. This system can be solved by conventional matrix methods. Indeed, the right-hand side when x = c m is simply a constant, and the left-hand side forms the equation with all A 1,m variables. In other words: where where again the c m are arbitrary values that are not equal to the zeros of the denominator polynomial. We can now return to Equation (14) knowing now that the A constants are all known.
The first term is the identical integral to the previous section. The second integral is again obtained by Lemma 1, but this time for a = 2. Doing this gives the final result: Figure 2 shows one particular numerical example validating the partial fraction coefficient calculation algorithm, by showing agreement between the numerical integration and the overall closed-form solution. The particular choices of numerator and denominator polynomials are the same as those of Figure 1, though in this case, the denominator is squared. The figure shows only the range of x that is between two singularities. The reason for this will be explained more thoroughly in Section 5. In short, for r > 1, integrating through multiple poles causes the resultant area under the curve to diverge.

Denominator of Arbitrary Total Exponentiation r
The above process, along with Lemma 1, can be easily extended to any denominator order r. Namely, extending the analysis to solve the integral: The critical step is the partial fraction decomposition. For r = 1, this amounts to one term for each root, shown in Equation (4). For r = 2 it amounts to two terms for each root, one with a linear denominator, and one squared, shown in Equation (15). This trend continues for arbitrary order r. Namely, for denominator order r, each binomial will yield r terms of increasing denominator order: The integration can proceed in an identical fashion as the previous section, with continued use of Lemma 1 for each value of k. Thus, the total generalized result is: Finding the coefficients A k,m can be done in an identical manner as the order 2 case. The final coefficient A r,m can be determined by first multiplying both sides of Equation (26) by (x − x j ) r : 1 and then substituting x = x j . Doing so causes the right-most term to vanish and allows A r,j to be easily found: To obtain all of the other coefficients A k,m , k = r, we again multiply both sides of Equation (26) by the entire denominator on the left-hand side (changing the dummy index again on the product): Inserting arbitrary values x = c m,k , c m,k = x m gives a system of equations of A k,m . Using the same notation as the previous section: where:

Spanning Multiple Poles in the Integration Range
When integrating over a pole, the underlying hypergeometric functions reach one of their singular points (argument equal to 1) at the pole itself and must be analytically continued for values of x/x m > 1. These analytic continuations occur in the complex plane, and a question arises of how these complex outputs contribute to the overall integration result, as it must be purely real since it represents an area under the curve of the integrand. We investigate this behavior by expanding the hypergeometric function: Performing an index shift, such that k = k + n + 1, and then renaming k to k gives: We now investigate the behavior of this expression for different values of r, starting with r = 1: By realizing that the Taylor expansion for the natural logarithm is given by one can then add and subtract the quantity ∑ n k=1 (z k /k) to Equation (38) and obtain: Recall that z = x/x m , and as such, when x approaches the value of the root x m , then z becomes equal to 1. As x continues (either positively or negatively, depending on the location of x m and the domain of interest), the value of this ratio becomes larger than 1. In this circumstance, the second term in Equation (39) must converge since the sum is finite. The argument for the first term then becomes negative. To allow for computation, one must extend the input values into the complex plane. For purely real inputs, the phase is simply equal to odd integer multiples of π, where we restrict the domain of values to be just π, which accounts for the negative sign. In this case, one can simply define y = 1 − z ∈ R | y < 0, and then write it in polar form with phase equal to π, namely, y = Re iπ . Substituting this into Equation (39) gives: We can now substitute it into Equation (12) with z = x/x m , and see that the prefactor term (n + 1)/x n+1 cancels. Upon summing, there will be a constant imaginary component. The consequence of this result is that the constant imaginary part can be subtracted away by the integration constant C if one imposes that C ∈ C, and the result becomes purely real.
If there exist complex roots, then Equation (39) must be treated differently. The complex roots will cause z = x/x m to be complex, and therefore yield some non-trivial phase value θ. Thus, we can write the argument of the logarithm as: The last step allows one to use the same procedure as the purely real zero case in regards to the complex logarithm. We must determine how the phase changes as z changes value. This is most easily done by first equating imaginary parts: and then solving for θ : Realizing that since R = Re(z) 2 + Im(z) 2 and R = (1 − Re(z)) 2 + Im(z) 2 , where Re(z) and Im(z) denote the real and imaginary parts of z, respectively, one can solve for Re(z) in the expression for R to get: This expression can then be substituted into the R expression and simplified to yield: Thus, we see from Equation (45), the imaginary part (dictated by the phase) will not remain constant, as was the previous case. However, because complex roots in real polynomials always come in conjugate pairs, the term corresponding to the conjugate root will have a phase of −θ. Substituting −θ into Equation (45) and using the fact that both sine and arcsin are odd functions results in an equal and opposite phase. Therefore, when summing the total contributions from both roots, the imaginary portions will cancel out completely and the resulting integration result will be purely real. Note that this analytic continuation choice does not include z = 0, and the divergence at this point (corresponding to x/x m = 1) cannot be removed, and thus, integration ranges should exclude these points.
In summary, these results indicate that for the r = 1 case, as long as one does not include the exact singular points, there is no need to worry about divergence. Furthermore, there will be no complex output, even though the complex nature of the analytic continuation is a necessity. Next we examine the r > 1 case: Referring back to Equation (36), if r > 1, the expression becomes: The computational implementation of the hypergeometric function for this article applies some but not all of these formulas automatically, and restricts the divergent nature to just the usual three points when possible.

Closing Remarks
This paper showcased how to integrate a generalized ratio of polynomials, where the denominator as a whole is raised to the general integer order r ∈ Z + . The results were then studied to investigate the possible complex output due to the analytic continuation of the hypergeometric function. Examples for r = 1 and r = 2 were considered and validated using numerical integration. For all cases, poles will cause divergence, but only at the value of the pole itself with the proper choice of analytic continuation. This issue can be circumvented if one omits the exact location of the poles in the integration region.
The methods described here are valid for any numerator or denominator polynomials, except in cases where the denominator has degenerate roots. In such a case, however, the same methodology presented here can be used, the only difference being the particular form of the partial fraction decomposition. The generalized approached to include denominator polynomials with degenerate roots is left as future work.
In regard to the exponent r, as mentioned previously, the analysis is limited to r ∈ Z + . However, it may be possible to extend this work to non-integer values of r. Whether or not the above techniques, namely, partial fraction decomposition, can be generalized to include this scenario, remains an open question left for future work.
Finally, it is important to mention the sensitivity of the closed-form solution with respect to calculating the roots of the denominator polynomial. Figure 3 shows one particular choice of numerator and denominator polynomial for differing amounts of error with respect to the location of the zeros. The image shows a domain between two roots x = 1 and x = 2. As can be seen, the error is minimal as one gets further from the poles. However, the error is substantial close to the poles. Considerations such as these should be taken into account over the required integration domain of the problem at hand.