A nasty integral of a rational function

Question

I'm having a hard time proving the following $$\int_0^{\infty} \frac{x^8 - 4x^6 + 9x^4 - 5x^2 + 1}{x^{12} - 10 x^{10} + 37x^8 - 42x^6 + 26x^4 - 8x^2 + 1} \, dx = \frac{\pi}{2}.$$

Mathematica has no problem evaluating it while I haven't the slightest idea how to approach it. Of course, I would like to prove it without the use of a computer. Is this an explicit form of a special function I fail to recognize?

Answer 1

The integral over $\mathbb{R}$ of a meromorphic function $f(z)$, $O(|z|^{-2})$ at infinity, non-vanishing over $\mathbb{R}$, is equal to $2\pi i$ times the sum of residues in the poles located in the complex upper-half plane. Since:

$$p(y) = y^6-2y^5-2y^4+4y^3+3y^2-4y+1 = p_{+}(y)\cdot p_{-}(y),$$ $$p_{+}(y)= y^3-(i+1)y^2+(i-2)x+1,\qquad p_{-}(y)=y^3+(i-1)y^2-(2+i)y+1,$$

(I got this through a numerical calculation of the roots of $p(y)$, followed by a separation of the roots with positive and negative imaginary part, say $\zeta_1,\zeta_2,\zeta_3$ and $\bar{\zeta_1},\bar{\zeta_2},\bar{\zeta_3}$ - so $p_{+}(z)$ is just $\prod_{j=1}^3 (z-\zeta_j)$) we have:

$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = 2\pi i\sum_{j=1}^{3}\operatorname{Res}_{z=\zeta_j}\left(\frac{z^2}{p_{+}(z)\cdot p_{-}(z)}\right),$$

but $p_{-}(x)-p_{+}(x)=2i(x^2-x)$, so:

$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = 2\pi i\sum_{j=1}^{3}\operatorname{Res}_{z=\zeta_j}\left(\frac{z^2}{p_{+}^2(z)+2i(z^2-z)p_{+}(z)}\right).$$

By De l'Hopital theorem, and since $\zeta_j$ is a double zero of $p_{+}^2(x)$:

$$\lim_{z\to\zeta_j}\frac{z^2(z-\zeta_j)}{p_{+}^2(z)+2i(z^2-z)p_{+}(z)}=\frac{\zeta_j^2}{2i(\zeta_j^2-\zeta_j)p_{+}'(\zeta_j)},$$

so:

$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = \pi\sum_{j=1}^{3}\frac{\zeta_j}{(\zeta_j-1)p_{+}'(\zeta_j)}.$$

Now we compute the remainder between $(z-1)p_{+}'(z)$ and $p_{+}(z)$, in order to have:

$$ I = \int_{\mathbb{R}}\frac{y^2 dy}{p(y)} = \pi\sum_{j=1}^{3}\frac{\zeta_j}{-(1+i)+6\zeta_j-(2-i)\zeta_j^2}.$$

If now we take $\alpha=\frac{3+\sqrt{6-i}}{2-i}$ and $\beta=\frac{3-\sqrt{6-i}}{2-i}$ we can re-write the last line as:

$$ I = \frac{\pi}{(i-2)(\alpha-\beta)}\left(\sum_{j=1}^{3}\frac{\alpha}{\zeta_j-\alpha}-\sum_{j=1}^{3}\frac{\beta}{\zeta_j-\beta}\right)=-\frac{\pi}{2\sqrt{6-i}}\left(\Sigma_1-\Sigma_2\right).$$

Now $\Sigma_1$ is the sum of the reciprocal of the roots of the polynomial $p_{+}(\alpha(z+1))$, and $\Sigma_2$ is the sum of the reciprocal of the roots of the polynomial $p_{+}(\beta(z+1))$. This quantities can be computed through the coefficients of $p_{+}$, since the sum of the reciprocal of the roots of a polynomial $q(z)$ is just $-\frac{q'(0)}{q(0)}$. This gives:

$$\Sigma_1 = -\alpha\frac{p_{+}'(\alpha)}{p_{+}(\alpha)},\qquad \Sigma_2 = -\beta\frac{p_{+}'(\beta)}{p_{+}(\beta)}.$$

Up to a massive amount of long but straightforward computations, we get:

$$\Sigma_1 = (i-2)-\sqrt{6-i},\qquad \Sigma_2 = (i-2)+\sqrt{6-i}, $$

from which $\color{red}{I=\pi}$ finally follows.

I am really grateful to Jon Haussmann for the proof that

$$\int_0^{\infty} \frac{x^8 - 4x^6 + 9x^4 - 5x^2 + 1}{x^{12} - 10 x^{10} + 37x^8 - 42x^6 + 26x^4 - 8x^2 + 1}dx = \frac{1}{2}\int_{\mathbb{R}}\frac{y^2 dy}{p(y)},$$

where only the second integral is treated here.

IMPORTANT UPDATE: In fact, there is no need to compute the coefficients of $p_{+}(x)$ and $p_{-}(x)$ (we only need the identity $p_{-}(x)-p_{+}(x)=2i(x^2-x)$), or introduce $\alpha$ and $\beta$. Since $p_{+}(x)$ is a third-degree polynomial with roots in the upper half-plane, $$0=\int_{\mathbb{R}}\frac{dz}{p_{+}(z)}=\sum_{j=1}^{3}\frac{1}{p_{+}'(\zeta_j)}.$$ This gives: $$ I = \pi\sum_{j=1}^{3}\frac{\zeta_j}{(\zeta_j-1)p_{+}'(z)} = \pi\sum_{j=1}^{3}\frac{1}{(\zeta_j-1)p_{+}'(\zeta_j)}, $$ but if we decompose $\frac{1}{p_{+}(z)}$ in simple fractions, we get: $$\frac{1}{p_{+}(z)}=\sum_{j=1}^{3}\frac{1}{p_{+}'(\zeta_j)(z-\zeta_j)}, $$ so the magic gives: $$ I = -\frac{\pi}{p_{+}(1)}.$$ Since $p(x)=p_{+}(x)\cdot p_{-}(x)$, $p(1)=1$, $I\in\mathbb{R}^+$ and $p_{-}(1)$ is the conjugate of $p_{+}(1)$, $p_{+}(1)$ can be only $+1$ or $-1$, so $I=\pi$.

Answer 2

Some progess: The integrand actually decomposes as $$\frac{1}{2} \left( \frac{x^2 + 2x + 1}{x^6 + 4x^5 + 3x^4 - 4x^3 - 2x^2 + 2x + 1} + \frac{x^2 - 2x + 1}{x^6 - 4x^5 + 3x^4 + 4x^3 - 2x^2 - 2x + 1} \right).$$ Note that the second term is the same as the first term, except with $-x$ instead of $x$. Thus, with some substitutions, the integral becomes $$\frac{1}{2} \int_{-\infty}^\infty \frac{y^2}{y^6 - 2y^5 - 2y^4 + 4y^3 + 3y^2 - 4y + 1} \; dy.$$

Answer 3

Let,

$$I=\int_0^{\infty} \frac{x^8 - 4x^6 + 9x^4 - 5x^2 + 1}{x^{12} - 10 x^{10} + 37x^8 - 42x^6 + 26x^4 - 8x^2 + 1} \, dx$$

As noted by Jon Haussmann,

$$2I=\int_0^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx + \int_0^{\infty}\frac{(x+1)^2}{x^6 +4x^5 + 3x^4 - 4x^3 - 2x^2 + 2x + 1}dx$$

Perform the change of variable $x=\dfrac{1}{u-1}$ in the second integral,

$\begin{align}2I&=\int_0^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\int_1^{\infty} \frac{x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\\ &=\left(\int_0^1 \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\\ \int_1^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\right)+\\&\int_1^{\infty} \frac{x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\end{align}$

Perform the change of variable $x=\dfrac{u-1}{u}$ in the first integral of the latter equality,

$\begin{align}2I&=\int_1^{\infty} \frac{x^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\int_1^{\infty} \frac{(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx+\int_1^{\infty} \frac{x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\\ &=\int_1^{\infty} \frac{x^2+(x-1)^2+x^2(x-1)^2}{x^6-4x^5 + 3x^4 +4x^3 - 2x^2 -2x + 1}dx\\ &=\Big[\arctan\left(\dfrac{x^3-2x^2-x+1}{x(x-1)}\right)\Big]_0^{+\infty}\\ &=\pi \end{align}$

(Problem found in American Mathematical Monthly, vol. 112, april 2005.

Solution found in vol. 114)

Answer 4

Among integrals like this one, although some are doable through Glasser's Master Theorem, most are related to theorems about polynomials as follows.

---

Let $p(x)$ be a monic polynomial of degree $n$, whose zeros $\{x_k\}$ all lie in the half upper plane $\mathcal H$. Note that it must not be a real polynomial. Since $p(x)=\prod_{k}(x-x_k)$ , the series expansion holds $$ \frac1{p(1/x)}=\prod_k\frac x{1-x_kx}=x^n\sum_{l\ge0}h_lx^l\\ $$

Here $$ h_l(\{x_k\})=\sum_{i_1+\cdots+i_n=l,i_j\ge0}x_1^{i_1}x_2^{i_2}\cdots x_n^{i_l} $$ is the complete the complete symmetric polynomial.

Given $k\le n=\deg p$ an integer, apply residue theorem to the integral $$ \begin{align} &\text{P.V.} \int_{-\infty}^\infty\frac{x^{k-1}}{p(x)}dx \\ =&\lim_{R\to\infty}\int_{-R}^R\frac{x^{k-1}}{p(x)}dx \\ =&\oint_{\mathcal H}\frac{z^{k-1}}{p(z)}dz-\int_C \frac{z^{k-1}}{p(z)}dz \\ =&-2\pi i\text{Res}\left[\frac{z^{k-1}}{p(z)},\infty\right]-\pi\delta_{k,n} \\ =&2\pi i[x^{k}]\frac1{p(1/x)}-\pi\delta_{k,n} \\ =&\pi\Big(2h_{k-n}-\delta_{k,n}\Big) \\ =&\pi~\delta_{k,n} \end{align} $$ The contour is the upper large semi circle. The integral on the large arc vanishes except the case $k=n$ since the leading term $\dfrac{z^{k-1}}{p(z)}\sim \dfrac{z^{k-1}}{z^n}$ . Meanwhile, $h_{k-n}$ takes nonzero value only when $k\ge n$, and obviously $h_0=1$.

Hence, we can generally conclude $$ \text{P.V.}\int_{-\infty}^\infty\frac{x^{k-1}}{p(x)}d x=\begin{cases} 0 & 0<k<n\\[5pt] \pi i & k=n\\[5pt] \infty &k>n \end{cases} $$ When $k=n$, take the imaginary part to obtain a integral that converges to $\pi$, and we can remove the $\text{P.V.}$ sign.

Given an integral like this, all we have to do is the following:

1. Factor the denominator in the form $|p(x)|^2=p(x)p^*(x)$.
2. Verify all the roots of $p(z)$ lies on the same side of the real line, else the method is not valid. Here we suppose they are all in $\mathcal H$.
3. Rewrite the integrand as the imaginary part of a rational function $\dfrac{q(x)}{p(x)}$,
4. The result is $\pi[x^{\deg p-1}]q(x)$.

---

As for the OP's integral, first extend the range to $\mathbb R$ since the integrand is even, then notice $$ \frac{x^8-4 x^6+9 x^4-5 x^2+1}{x^{12}-10 x^{10}+37 x^8-42 x^6+26 x^4-8 x^2+1}=\Im\frac{{\color{red}1}x^5-2 i x^4-6 x^3+5 i x^2+3 x-i}{{\color{red}1}x^6-2 i x^5-7 x^4+6 i x^3+6 x^2-2 i x-1} $$ therefore the integral on $\mathbb R$ is $\pi$ and the desired one is half of it.

Another example $$ \frac{x^3}{x^8-3 x^6+2 x^5+7 x^4+12 x+13}=\Im\frac{\color{red}{-497}x^3+\cdots}{\color{red}{4148} \left(x^4-3 i x^3-(6+i) x^2-(2-6 i ) x+3+2i\right)} $$ so that $$ \int_{-\infty}^\infty\frac{x^3}{x^8-3 x^6+2 x^5+7 x^4+12 x+13}d x=-\frac{497}{1418}\pi $$ More $$ \int_{-\infty}^\infty\frac{x^4}{x^{10}+6 x^9+65 x^8+80 x^7-83 x^6-146 x^5+56 x^3+21 x^2+4 x+1}dx=\frac{5992 \pi }{53533} $$

$$ \int_{-\infty}^\infty\frac{x^7}{x^{10}+4 x^9-x^8-12 x^7+5 x^6+14 x^5-20 x^4-52 x^3+61 x^2+44 x+20}dx=-\frac{183167 \pi }{28240} $$

$$ \int_{-\infty}^\infty\frac{x^{13}}{p(x)}dx=\frac{13882987895674304464657174937266 \pi }{103271122611655554792443174968701227}\\ p(x)=100 x^{16}-180 x^{15}+990 x^{14}-818 x^{13}+3507 x^{12}-2376 x^{11}+7767 x^{10}-1952 x^9\\ +9354 x^8+1684 x^7+10338 x^6+8556 x^5+9421 x^4+8520 x^3+7483 x^2-3138 x+1066 $$

As for an arbitrary polynomial, one can find the lower bound of the imaginary part of all its roots, and shift it to $\mathcal H$ to obtain a suitable polynomial, so it is very easy to come up with integrals of this type.

Remark 1

The only part that lacks rigor is the verification of the locations of the roots. I am more of a physicist so the numerical result satisfies me. To be more rigorous, one can estimate the maximum magnitude of the roots $R_\max$, then do the numerical integral $$ \frac1{2\pi i}\oint_{C}\frac{p'(z)}{p(z)}dz $$ along the upper semi circle with radius exceeding $R_\max$. The result is always a integer, and should be $\deg p$. For a pure theoretical proof, I suppose it requires tedious estimation or advanced polynomial theorem which is unfamiliar to me.

Remark 2

Actually, we can apply the procedure to all real polynomials without no real roots since their roots are symmetrical to the real axis. However, only when the factorization requires algebraics less complicated than the roots of the polynomial simplifies the problem. Indeed, in most cases we just need $i$ (so the result is just a rational multiple of $\pi$).

---

I wrote a Mathematica code to evaluate similar integrals

rint::fef = "Method not suitable";
Options[rint] = {Extension -> {I}};
rint[p_, x_, OptionsPattern[]] := Module[{num, den, prfl, apl},
  If[Limit[x p, x -> Infinity] != 0, 
   Message[Integrate::"idiv", p, {-\[Infinity], \[Infinity]}],
   {num, den} = NumeratorDenominator[Simplify@p];
   If[PolynomialQ[num, x] && PolynomialQ[den, x],
    prfl = 
     Select[Rest@FactorList[den, Extension -> OptionValue[Extension]],
       And @@ (Solve[#[[1]] == 0, x] // Values // Flatten // Im // 
          Positive) &];
    If[Length[prfl] == 0, Message[rint::fef],
     apl = 
      List @@ (Apart@Factor[p, Extension -> OptionValue[Extension]]);
     2 Pi I Sum[
         Cases[t[[1]]/Coefficient[t[[1]], x, Exponent[t[[1]], x]] apl,
           X_ /; PolynomialQ[X, x] :> 
           Coefficient[X, x, Exponent[X, x]]], {t, prfl}][[1]] // 
      Simplify
     ]
    , Message[rint::fef]]
   ]]

You can enter rint[(1 - 5 x^2 + 9 x^4 - 4 x^6 + x^8)/( 1 - 8 x^2 + 26 x^4 - 42 x^6 + 37 x^8 - 10 x^10 + x^12), x] to get OP's integral, or play around with similar ones on the site.

The algorithm will output the result or state that the method fails. If the latter occurs and the integral indeed converges, you may analyze the roots and possibly find some simple algebraics to do the factorization, then simply add it to the Extension option. For example, rint[x^2/(1 + 6 x^2 + x^4)^2, x] fails, while rint[x^2/(1 + 6 x^2 + x^4)^2, x, Extension -> {I, Sqrt[2]}] does the job.

I am not good at Mathematica, so if you find any flaws or have any suggestions please inform me.