How to find exact value of integral $\int_{0}^{\infty} \frac{1}{\left(x^{4}-x^{2}+1\right)^{n}}dx$?

Question

When I first encountered the integral $\displaystyle \int_{0}^{\infty} \frac{1}{x^{4}-x^{2}+1} d x$, I am very reluctant to solve it by partial fractions and search for any easier methods. Then I learnt a very useful trick to evaluate the integral.

Noting that $$I(1):= \int_{0}^{\infty} \frac{d x}{x^{4}-x^{2}+1} \stackrel{x \mapsto \frac{1}{x}}{=} \int_{0}^{\infty} \frac{x^{2}}{x^{4}-x^{2}+1} $$

Combining them yields

\begin{aligned} I(1)&=\frac{1}{2} \int_{0}^{\infty} \frac{x^{2}+1}{x^{4}-x^{2}+1} d x\\&= \frac{1}{2}\int_{0}^{\infty} \frac{1+\frac{1}{x^{2}}}{x^{2}+\frac{1}{x^{2}}-1} d x \\ &= \frac{1}{2}\int_{0}^{\infty} \frac{d\left(x-\frac{1}{x}\right)}{\left(x-\frac{1}{x}\right)^{2}+1} \\ &= \frac{1}{2}\tan ^{-1}\left(x-\frac{1}{x}\right)_{0}^{\infty} \\ &= \frac{\pi}{2} \end{aligned}

Later, I started to investigate the integrands with higher powers.

Similarly, $$ I(2):= \int_{0}^{\infty} \frac{d x}{\left(x^{4}-x^{2}+1\right)^{2}} \stackrel{x \mapsto \frac{1}{x}}{=}\int_{0}^{\infty} \frac{x^{6}}{\left(x^{4}-x^{2}+1\right)^{2}} d x $$

By division, we decomposed $x^6$ and obtain $$ \frac{x^{6}}{\left(x^{4}-x^{2}+1\right)^{2}}=\frac{x^{2}+1}{x^{4}-x^{2}+1}-\frac{1}{\left(x^{4}-x^{2}+1\right)^{2}} $$$$ I(2)=\int_{0}^{\infty} \frac{x^{2}+1}{x^{4}-x^{2}+1} d x-\int_{0}^{\infty} \frac{1}{\left(x^{4}-x^{2}+1\right)^{2}}dx $$ We can now conclude that $$I(2)=I(1)=\frac{\pi}{2} $$

My Question:

How about the integral $$\displaystyle I_{n}=\int_{0}^{\infty} \frac{d x}{\left(x^{4}-x^{2}+1\right)^{n}}$$ for any integer $n\geq 3$?

Answer 1

In order to tackle the general integral, I define a more general integral

$$ I_{n}(a):=\int_{0}^{\infty} \frac{1}{\left(x^{4}-a x^{2}+1\right)^{n}} d x. $$ where $a<2.$

I am attempting to create a “Reduction Formula” for $I_{n}(a)$.

Differentiating $I_{n}(a)$ twice w.r.t. $a$ yields $$ \int_{0}^{\infty} \frac{x^{4}}{\left(x^{4}-a x^{2}+1\right)^{n+2}} d x=\frac{1}{n(n+1)} I_{n}^{\prime \prime}(a) \tag*{(1)} $$

Differentiating $I_{n+1}(a)$ w.r.t. $a$ yields $$ \int_{0}^{\infty} \frac{x^{2}}{\left(x^{4}-a x^{2}+1\right)^{n+2}} d x=\frac{1}{n+1} I_{n+1}^{\prime}(a) \tag*{(2)} $$

(1)-$a \times $(2) yields

$$ \begin{aligned} \\ \int_{0}^{\infty} \frac{x^{4}-a x^{2}+1-1}{\left(x^{4}-a x^{2}+1\right)^{n+2}} d x=\frac{I_{n}^{\prime \prime}(a)}{n(n+1)}-\frac{a I_{n+1}^{\prime}(a)}{n+1} \\ I_{n+1}(a)-I_{n+2}(a)=\frac{I_{n}^{\prime \prime}(a)}{n(n+1)}-\frac{a I_{n+1}^{\prime}(a)}{n+1} \\ \boxed{I_{n+2}(a)=I_{n+1}(a)-\frac{I_{n}^{\prime \prime}(a)}{n(n+1)}+\frac{a I_{n+1}^{\prime}(a)}{n+1}}, \\ \end{aligned} $$ which is the “Reduction Formula” for $I_n(a).$ Now we need two “initial” integrals, $I_1(a)$ and $I_2(a)$ which can be found in my answer in Quora,

$\displaystyle \boxed{I_{1}(a)=\frac{\pi}{2 \sqrt{2-a}} \text { and } I_{1}^{\prime}(a)=\frac{\pi}{4(2-a)^{\frac{3}{2}}}} \tag*{} $

$\boxed{ \displaystyle I_{2}(a)=\frac{\pi(3-a)}{4(2-a)^{\frac{3}{2}}} \text { and } I_{2}^{\prime}(a)=\frac{\pi(5-a)}{8(2-a)^{\frac{5}{2}}}} \tag*{} $

Now we can apply the formula. $$ \begin{aligned} I_{3}(a) &=I_{2}(a)-\frac{I_{1}^{\prime \prime}(a)}{2}+\frac{a I_{2}^{\prime}(a)}{2} \\ &=\frac{\pi(3-a)}{4(2-a)^{\frac{3}{2}}}-\frac{3 \pi}{16(2-a)^{\frac{5}{2}}}+\frac{a \pi(5-a)}{16(2-a)^{\frac{5}{2}}}\\&= \frac{3 \pi\left(a^{2}-5 a+7\right)}{16(2-a)^{\frac{5}{2}}} \end{aligned} $$ $$\boxed{\int_{0}^{\infty} \frac{1}{\left(x^{4}-a x^{2}+1\right)^{3}} d x= \frac{3 \pi\left(a^{2}-5 a+7\right)}{16(2-a)^{\frac{5}{2}}}}\tag*{}$$

In particular, when $a=1$,

$$ \boxed{I_{n+2}(1)=I_{n+1}(1)-\frac{I_{n}^{\prime \prime}(1)}{n(n+1)}+\frac{I_{m+1}^{\prime}(1)}{n+1}} $$

$$ \int_{0}^{\infty} \frac{d x}{\left(x^{4}-x^{2}+1\right)^{3}}=\frac{9 \pi}{16}, $$ $$ \begin{aligned} \int_{0}^{\infty} \frac{d x}{\left(x^{4}-x^{2}+1\right)^{4}} &=I_{3}(1)-\frac{I_{2}^{\prime \prime}(1)}{6}+\frac{I_{3}^{\prime}(1)}{3} =\frac{9 \pi}{16}-\frac{9 \pi}{48}+\frac{9 \pi}{32} =\frac{21 \pi}{32} \end{aligned} $$

By the way, I had also found 8 other integrals. $ \begin{aligned}&\int_{0}^{\infty} \frac{d x}{x^{4}+x^{2}+1}=I_{1}(-1)=\frac{\pi}{2 \sqrt{3}} \\&\int_{0}^{\infty} \frac{d x}{x^{4}-x^{2}+1}=I_{1}(1)=\frac{\pi}{2} \\&\int_{0}^{\infty} \frac{d x}{\left(x^{4}+x^{2}+1\right)^{2}}=I_{2}(-1)=\frac{\pi}{3 \sqrt{3}} \\&\int_{0}^{\infty} \frac{d x}{\left(x^{4}-x^{2}+1\right)^{2}}=I_{2}(1)=\frac{\pi}{2} \\&\int_{0}^{\infty} \frac{d x}{\left(x^{4}+x^{2}+1\right)^{3}}=I_{3}(-1)=\frac{13}{48 \sqrt{3}} \\&\int_{0}^{\infty} \frac{d x}{x^{4}+1}=I_{1}(0)=\frac{\pi}{2 \sqrt{2}} \\&\int_{0}^{\infty} \frac{d x}{\left(x^{4}+1\right)^{2}}=I_{2}(0)=\frac{3 \pi}{8 \sqrt{2}} \\&\int_{0}^{\infty} \frac{d x}{\left(x^{4}+1\right)^{3}}=I_{3}(0)=\frac{21\pi}{64 \sqrt{2}}\end{aligned} \tag*{} $

Theoretically, I can find the exact values of the integrals one by one. But is there any hope to find the closed forms of $I(a)$ using the formula? Your help is highly appreciated.

Answer 2

Too long for a comment

You could be interested by $$\displaystyle I_{n}=\int_{0}^{t} \frac{d x}{\left(x^{4}-x^{2}+1\right)^{n}}=t \, F_1\left(\frac{1}{2};n,n;\frac{3}{2};e^{+\frac{i \pi }{3}} t^2,e^{-\frac{i \pi }{3}} t^2\right)$$ where appears the Appell hypergeometric function of two variables and the roots of the quadratic in $x^2$.

I did not find any interesting reduction formula or any asymptotics for $t \to \infty$.

Nevertheless, $$\displaystyle J_{n}=\int_{0}^{\infty} \frac{d x}{\left(x^{4}-x^{2}+1\right)^{n}}=\pi \,\frac {a_n}{2^{b_n}}$$ where the $a_n$ correspond to the sequence $$\{1,1,9,21,25,483,2359,2907,115533,288915,363363,3674619,37326471,\cdots\}$$ and the $b_n$ the sequence $$\{1,1,4,5,5,9,11,11,16,17,17,20,23,\cdots\}$$ None of these has been found in $OEIS$.

Edit

For $$I_{m,n}=\int_0^\infty\frac{dx}{(x^m+1)^n}=\frac{\Gamma \left(1+\frac{1}{m}\right) \Gamma \left(n-\frac{1}{m}\right)}{\Gamma(n)} \qquad \text{if} \qquad \Re(m n)>1\land \Re(m)>0$$

gives the simple $$I_{m,n+1}=\left(1-\frac{1}{m n}\right)I_{m,n}$$

Answer 3

Note

$$ \int_0^\infty \frac1{(x^4-x^2+1)^{n+1}}dx=\frac{(-1)^{n}}{n!} \frac{d^{n} J(a)}{da^{n}}\bigg|_{a=1} $$ where $$J(a)=\int_0^\infty \frac{dx}{x^4-x^2+a}=\frac1{2\sqrt a} \int_0^\infty \frac{d(x-\frac{\sqrt a}x)}{(x-\frac{\sqrt a}x)^2+2\sqrt a-1} = \frac\pi{2\sqrt{a(2\sqrt a-1)} } $$

Then \begin{align} & \int_0^\infty \frac1{x^4-x^2+1}dx= J(1)=\frac\pi2\\ &\int_0^\infty \frac1{(x^4-x^2+1)^2}dx= -\frac{dJ(a)}{da}\bigg|_{a=1}=\frac\pi2\\ &\int_0^\infty \frac1{(x^4-x^2+1)^3}dx=\frac12\frac{d^2J(a)}{da^2}\bigg|_{a=1}=\frac{9\pi}{16}\\ &\int_0^\infty \frac1{(x^4-x^2+1)^4}dx=-\frac16\frac{d^3 J(a)}{da^3}\bigg|_{a=1}=\frac{21\pi}{32}\\ &\int_0^\infty \frac1{(x^4-x^2+1)^5}dx=\frac1{24} \frac{d^4 J(a)}{da^4}\bigg|_{a=1}=\frac{25\pi}{32}\\ &\int_0^\infty \frac1{(x^4-x^2+1)^6}dx=-\frac1{120} \frac{d^5 J(a)}{da^5}\bigg|_{a=1}=\frac{483\pi}{512}\\ &\hspace{5mm}\cdots\hspace{2mm}\cdots \end{align}

Answer 4

There are many approaches to solving this integral that I found in papers. We first generalize the integral to $$ N_{0,4}(a ; m) = \int_0^\infty \frac{dx}{(x^4+2ax^2+1)^{m+1}} $$

The formula for this general expression is $$ N_{0,4}(a ; m) = \frac{\pi}{2} \frac{P_m(a)}{[2(a+1)]^{m+\frac 12}} $$

where $$ P_m(a) = \sum_{l=0}^m \left[2^{-2m} \sum_{k=l}^m 2^k \binom{2m-2k}{m-k}\binom{m+k}{m}\binom{k}{l}\right]a^l $$

I will be providing brief details on these computations , and not much more than that , in order to highlight merely the key computations. Some sections overlap with each other, but each one of them is, in their own right, capable of solving and generating further discussion.

---

1

The first approach is to set $x = \tan \theta$. This gives $$ \int_{0}^{\frac \pi 2} \frac{d\theta}{\cos^2 \theta}\left(\frac{\cos^4 \theta}{\sin^4 \theta + 2a \sin^2 \theta \cos^2 \theta + \cos^4 \theta}\right)^{m+1} $$

Observe that $$ \sin^4 \theta + 2a \sin^2 \theta \cos^2 \theta + \cos^4 \theta = (1+a)+(1-a)\cos^2 (2 \theta) $$

We set $u = 2 \theta$ now: $$ 2^{-(m+1)} \int_0^\pi \left(\frac{(1+\cos u)^2}{(1+a) + (1-a)\cos^2 u}\right)^{m+1} \frac{du}{1+\cos u} $$

which simplifies to $$ 2^{-(m+1)} \int_0^{\pi} \left(\frac{(1+\cos u)^{2m+1}}{[(1+a) + (1-a)\cos^2 u]^{m+1}}\right) du $$

Now we expand $(1+\cos u)^{2m+1}$ using the binomial theorem. Since this is a linear combination of powers of $\cos u$, we'll merely show how to tackle the more general $$ \int_0^{\pi} [(1+a) + (1-a)\cos^2 u]^{-(m+1)} \cos^j u \ du $$

The first observation is that for $j$ odd, the integrand is an odd function around $\frac \pi 2$. Therefore, the corresponding integral is zero. For $j$ odd, one applies the double angle substitution again since we have a polynomial in $\cos^2$ as the integrand, and notes that the resulting function following substitution is even, so a factor of $2$ comes out and the integral limits change. The end of this should lead to $$ N_{0,4}(a ; m) = \sum_{j=0}^m 2^{-j} \binom{2m+1}{2j} \int_0^{\pi} [(3+a)+(1-a)\cos v]^{-(m+1)} (1+\cos v)^{j} dv $$ Finally, we can simplify this using the typical $\tan(v/2)$ substitution, and the result is already specified above.

---

2

We can find a recursion as follows. Let $V(x) = x^4+2ax^2+1$. The idea is to recognize that we can generate the recursion by integrating a very simple polynomial identity. We start by noticing that $V,V'$ are co-prime, therefore we can find polynomials $B,C$ such that $CV+BV' = \frac{-1}{m}$. A quick use of the inverse Euclid algorithm will tell you $$ B(x) = \frac{ax^3-x(1-2a^2)}{4(a^2-1)m} \quad ; \quad C(x) = \frac{\frac{ax^2}{1-a^2} - 1}{m} $$

Once we divide both sides of the equation by $V^{m+1}$ we get $$ CV^{-m} + BV'V^{-(m+1)} = \frac{-V^{-(m+1)}}{m} $$ and integrating this from $0$ to $\infty$ gives $$ \int_0^\infty CV^{-m} dx + \int_0^\infty BV'V^{-(m+1)} dx = \frac{-1}{m} \int_0^\infty V^{-(m+1)} dx $$

Which leads to $$ -\int_0^\infty CmV^{-m} dx - \int_0^\infty BmV'V^{-(m+1)} dx = N_{0,4}(a ; m) $$

In the first term, one can use the definition of $C$ to write the integral as a sum of two integrals, one involving $x^2V^{-m}$ and the other involving just $V^{-m}$. For the other, one calculates $BmV'$ ,performs the polynomial division $\frac{BmV'}{V}$ into quotient and remainder and then splits the integral based on that calculation, leading once again to a sum of two integrals : one involving $x^2V^{-m}$ and the other involving $V^{-m}$. Once these calculations are performed, we are eventually led to $$ N_{0,4}(a ; m) = \left(1+\frac{1-2a^2}{4m(a^2-1)}\right) N_{0,4}(a ; m-1) + \frac{(4m-3)a}{4m(a^2-1)} \int_{0}^\infty x^2V^{-m} dx $$

However, we now make the observation that $$ \int_{0}^\infty x^2V^{-m} dx = \frac{1}{2(m-1)}\frac{d}{da} \int_0^\infty V^{-(m-1)} dx $$ Since the derivative and the integral can be interchanged. This leads to the following recursion $$ N_{0,4}(a ; m) = \left(1+\frac{1-2a^2}{4m(a^2-1)}\right) N_{0,4}(a ; m-1) -\frac{(4m-3)a}{8m(m-1)(a^2-1)}\frac{d}{da} N_{0,4}(a ; m-2) $$

Now, one can check that the polynomial given as an answer, satisfies these conditions and the initial conditions, completing the proof. Having said that, the intuition behind why such a polynomial should come into play is hard to explain without knowing hypergeometric functions. Extremely roughly speaking , however : provided we latch on a multiplicative term involving $a$, we can find a cousin of the solution to the recurrence above, that satisfies a second-order recurrence in $m$, which is also of the orthogonal type. That leads us to the Jacobi polynomials and the resulting expansion.

---

3

We are now going to use reduction formulas that express the given integrand in terms of simpler, easier to integrate functions. The theory applies to the integration of even rational polynomials with positive coefficients, whose denominator has symmetric coefficients i.e. reads the same from (Like $(x^4+2ax^2+1)^m$)

Theorem : Define the "general symmetric polynomial" $$ D_p(d_1,\ldots,d_p ; z)=\sum_{k=0}^p d_{p+1-k}(z^{2k}+z^{4p-2k}) $$ Define the auxiliary polynomials $$ E_p(d_1,\ldots,d_p;z) = \left(\sum_{j=1}^{p} d_j\right) z^{2p} + \sum_{i=1}^p 2^{2i-1} z^{2(p-i)}\sum_{k=1}^{p-i+1}\frac{j+i-1}{i}\binom{j+2i-2}{j-1}d_{j+i} $$ For $d_i \in \mathbb R^+$. We have the following integral identity $$ \int_0^{\infty} \frac{z^{2n}dz}{D_p^{m+1}} = 2^{-m}\sum_{j=0}^{(m+1)p-n-1} 4^j\binom{(m+1)p-n-1+j}{2j}\int_0^{\infty} \frac{z^{2[(m+1)p-1-j]}dz}{E_p^{m+1}} $$

In fact, the proof isn't really very difficult at all! It's actually just a repetition of what was done in section $1$ : a couple of double-angle substitutions topped off by a $\tan(v/2)$ type substitution. Nevertheless, it is important to see that the generalization occurs in the first place.

Using this theorem on $D_1 = z^4+2d_1z^2+1$ and the computed $E_1 = (1+d_1)z^2+2$ gives the following identity with $n=0$ : $$ \int_0^{\infty} \frac{dz}{(z^4+2d_1z^2+1)^{m+1}} = 2^{-m}\sum_{j=0}^m 4^j\binom{m+j}{2j} \int_0^\infty \frac{z^{2(m-j)}dz}{((1+d_1)z^2+2)^{m+1}} $$

A change of variable involving $z^2$ on the latter integral will tell you that you only need to be able to evaluate $\displaystyle\int_0^\infty \frac{u^{r-\frac 12}du}{(1+u)^s}$. This can be done in a few ways, perhaps the best being to let $t = (1+u)^{-s}$, then recognizing a Beta integral which simplifies to a ratio of Gamma functions, followed by using Euler's reflection formula for the $\Gamma$ function. All in all, it's a simpler integrand to take care of.

---

4

We now introduce yet another transform. This transform is called the Landen transform, and it's origin in this situation is observed when one uses a substitution relating to $\cot \theta$. You can try , first, to prove the following identity

Theorem : Let $f$ be a rational function which is integrable over $\mathbb R$. Then $$ \int_{-\infty}^{\infty} f(x)dx = \int_{-\infty}^{\infty} \left[f(y+\sqrt{y^2+1}) + f(y - \sqrt{y^2+1})\right]\left(1+\frac{y}{\sqrt{y^2+1}}\right)dy $$

If $f$ is even, then the above integral limits can also be $0,+\infty$. However, the most important observation is that this transform retains rationality : on the RHS, you still end up getting a rational function to integrate. The hope, of course, is that it will be simple enough to integrate.

In our case, if $Q(x) = \frac 1{(x^4+2ax^2+1)^{m+1}}$ then we use the result once to get the integrand on the right hand side $$ Q_1 = \left[Q(y+\sqrt{y^2+1}) + Q(y - \sqrt{y^2+1})\right]\left(1+\frac{y}{\sqrt{y^2+1}}\right) $$

Remarkably enough, due to the nature of the $Q$ present here, we can simplify this and get the following expression $$ Q_1 = \sum_{k=0}^m \binom{m+k}{m-k} \frac{(2y)^{2k}}{2^m(1+a+2y^2)^{m+1}} $$

Without going too much into that particular proof, the idea is that if $\phi = y + \sqrt {y^2+1}$ then $-\phi^{-1} = y - \sqrt{y^2+1}$. One can find a relation between $Q_1$ and the polynomial $\frac{\phi^{2m+1}+\phi^{-(2m+1)}}{\phi+\phi^{-1}}$. The latter polynomial is reminiscent of polynomials involving roots of quadratic equations : it is seen to satisfy a typical three-term recurrence. This three-term recurrence, interestingly enough, is once again of the "orthogonal form" as mentioned before, and the RHS comes in as a solution of the recurrence.

However, once $Q_1$ is identified, it's clear that the integral of $Q_1$ again reduces to the integration of $\frac{y^mdy}{(1+y)^n}$, and we are done.

---

5

We use Ramanujan's master theorem now, to provide yet another way of evaluating the integral. However, this one is more surprising, because it yields the answer quite straightforwardly and without too much fuss. It is also capable of harbouring better intuition. While the theorem itself is more difficult to state, it is the following version that is of use, and easier to state.

Theorem : Let $f(x) = \sum_{k=0}^\infty \frac{\lambda(k)(-x)^k}{k!}$ be expandable in a series. Then $$ \int_0^\infty x^{s-1} f(x)dx = \Gamma(s)\lambda(-s) $$

That is, the ability to write $f$ as a power series in $-x$ allows you to then find certain integrals relating to $f$ (These are called the Mellin transforms). We want to use this to get our result. How do we do that?

The idea is to find a function , whose $\lambda$ are the $N_{0,4}(a;m)$. However, the integrands of $N_{0,4}(a ; m)$ form a geometric progression. It is then natural to expect that $f$ is likely to be of the type $\frac 1{1-\ldots}$ or something i.e. a sum, whose expansion is a geometric progression of coefficients.

With this in mind, we observe the following identity by solving the base case $$ \frac 1{\pi \sqrt 2} \int_0^\infty \frac{dx}{x^4+2ax^2+1+c} = \frac{d}{dc}\sqrt{a+\sqrt{1+c}} $$

Successive differentiations with respect to $c$ lead to an important geometric pattern being established. The result of such an operation is that we obtain the Taylor expansion $$ \sqrt{a+\sqrt{1+c}} = \sqrt{a+1} + \frac{1}{\pi \sqrt 2} \sum_{k=1}^\infty \frac{(-1)^{k-1}}{k} N_{0,4}(a ; k-1)c^k $$

and we are in the territory of the theorem. Applying that, tells us $$ \Gamma(s) \lambda(-s) = \int_0^\infty c^{s-1}\sqrt{a+\sqrt{1+c}} dc $$ However, did you notice that $\lambda(-s) = -\frac{(-s-1)!}{\pi \sqrt 2} N_{0,4}(a ; -s-1)$ seems to be some negative indexed analogue of $N_{0,4}$? This might be somewhat disconcerting. However, we are in the realm of complex analysis , where analytic continuation allows us to make meaning of such terms. One may notice the following derivative identity to achieve this $$ \lambda(m+1) = \left(\frac{-1}{2} \frac{d}{da}\right)^{2m+1} \lambda(-m) $$ Following the differentiation and some routine substitutions, we arrive at $$ N_{0,4}(a ; m) = \frac{m \pi \sqrt 2}{2^{6m+1}}\binom{4m}{2m}\binom{2m}{m}\int_1^\infty u(u^2-1)^{m-1}(a+u)^{-(2m+1/2)}du $$

The latter may be integrated by-parts : indeed, we take $u(u^2-1)^m$ to be the differentiated term, and integrate the other term repeatedly. The derivatives of this function eventually run out to $0$ . Furthermore, we can find the values of its derivatives at the point $1$ until they run out, hence we can derive every boundary term occurring in the expansion. Performing these routine calculations then lead to the proof.

---

6

This is just Gradshteyn-Ryzhik and their formula for this case, using the hypergeometric function. Formula 3.252.11 of GS reads $$ \int_0^\infty \frac{x^{\nu-1}dx}{(x^2+2ax+1)^{\mu + 1/2}} = \frac{2^\mu \Gamma(1+\mu)B(-\nu+2\mu+1,\nu) P^{-\mu}_{\mu - \nu}(a)}{(1-a^2)^{\mu/2}} $$

where $\Gamma,B$ are the Gamma and Beta functions respectively, and $$ P^{\beta}_{\alpha}(x) = \frac 1{\Gamma(1-\beta)} \left(\frac{a+1}{a-1}\right)^{\beta/2} {}_2F_1(-\alpha,\alpha+1 ; 1-\beta ; \frac{1-a}{2}) $$

is the Legendre function with parameters $\alpha,\beta$ expressed in parametric form. It is possible to use this representation to connect it to Jacobi polynomials, which leads to the formula in the starting.

I do not know how to derive this identity or how the hypergeometric functions link up with the Jacobi functions, however.

---

References

• Gradshteyn, Ryzhik : Table of integrals, series and products

• Boros, George; Moll, Victor H., An integral hidden in Gradshteyn and Ryzhik, J. Comput. Appl. Math. 106, No. 2, 361-368 (1999). ZBL0939.33007

• Boros, George; Moll, Victor H., The double square root, Jacobi polynomials and Ramanujan’s master theorem, J. Comput. Appl. Math. 130, No. 1-2, 337-344 (2001). ZBL1011.33005

• Amdeberhan, Tewodros; Espinosa, Olivier; Gonzalez, Ivan; Harrison, Marshall; Moll, Victor H.; Straub, Armin, Ramanujan’s master theorem, Ramanujan J. 29, No. 1-3, 103-120 (2012). ZBL1258.33001

• Amdeberhan, Tewodros; Moll, Victor H., A formula for a quartic integral: a survey of old proofs and some new ones, Ramanujan J. 18, No. 1, 91-102 (2009). ZBL1178.33002

• Boros, George; Moll, Victor H., Landen transformations and the integration of rational functions, Math. Comput. 71, No. 238, 649-668 (2002). ZBL0988.33009

Answer 5

Glad to see there are several nice answers in various perspectives. I am now going to give one more by inverse substitution followed by integration by parts.

Using $x \mapsto \frac{1}{x}$ yields $$ I(m, n, r):=\int_{0}^{\infty} \frac{x^{r} d x}{\left(x^{m}+1\right)^{n}} \stackrel{x \mapsto \frac{1}{x}}{=} \int_{0}^{\infty} \frac{\frac{1}{x^{r}}}{\left(\frac{1}{x^{m}}+1\right)^{n}}\left(\frac{d x}{x^{2}}\right) $$ Simplifying and then performing integration by parts, $$ \begin{aligned} I(m, n, r) &=\int_{0}^{\infty} \frac{x^{m n-r-2}}{\left(1+x^{m}\right)^{n}} d x \\ &=-\frac{1}{m(n-1)} \int_{0}^{\infty} x^{(n-1)m-r-1} d\left(\frac{1}{\left(1+x^{m}\right)^{n-1}}\right) \\ &=\frac{m(n-1)-r-1}{m(n-1)} \int_{0}^{\infty} \frac{x^{m(n-1)-r-2}}{\left(1+x^{m}\right)^{n-1}} d x \\ &=\left(1-\frac{r+1}{m(n-1)}\right) I(m, n-1, r) \\ &\qquad\qquad \vdots\\ &=\left(1-\frac{r+1}{m(n-1)}\right)\left(1-\frac{r+1}{m(n-2)}\right) \cdots\left(1-\frac{r+1}{m}\right) I\left(m,1,r\right) \\ &=\prod_{j=1}^{n-1}\left(1-\frac{r+1}{j m}\right)\int_{0}^{\infty} \frac{x^{r}}{x^{m}+1} d x \end{aligned} $$ Using the formula proved in my post $$ \int_{0}^{\infty} \frac{x^{r}}{x^{m}+1} d x=\frac{\pi}{m} \csc \frac{(r+1) \pi}{m}, $$

by putting $r=2k$ and $m=6$, we can conclude $$ \begin{aligned} \int_{0}^{\infty} \frac{1}{\left(x^{4}-x^{2}+1\right)^{n}} d x =& \int_{0}^{\infty}\left(\frac{x^{2}+1}{x^{6}+1}\right)^{n} d x \\ =& \sum_{k=0}^{n}\left(\begin{array}{c} n \\ k \end{array}\right) \int_{0}^{\infty} \frac{x^{2 k}}{\left(1+x^{6}\right)^{n}} d x \\=& \boxed{\sum_{k=0}^{n}\left[\left(\begin{array}{c} n \\ k \end{array}\right) \frac{\pi}{6} \csc \frac{(2 k+1) \pi}{6} \prod_{j=1}^{n-1}\left(1-\frac{2 k+1}{6 j}\right)\right]} \end{aligned} $$

Answer 6

Differentiation sometimes is easier than integration!

Inverse substitution rewrites the integral

$\displaystyle \int_{0}^{\infty} \frac{x^{r}}{x^{m}+a} d x \stackrel{x\mapsto\frac{1}{\sqrt[m]{a}}}{=}\frac{1}{a^{1-\frac{r+1}{m}}} \int_{0}^{\infty} \frac{x^{r}}{x^{m}+1} d x \tag*{} $ where $0<a<2.$

Using the formula proved in my post $ \displaystyle \int_{0}^{\infty} \frac{x^{r}}{x^{m}+1} d x=\frac{\pi}{m } \csc \frac{(r+1) \pi}{m},\tag*{} $

we get $ \displaystyle \int_{0}^{\infty} \frac{x^r}{x^{m}+a} d x=a^{-\left(1-\frac{r+1}{m}\right)} \frac{\pi}{m} \csc \frac{(r+1) \pi}{m}\tag*{} $

Differentiating the integral w.r.t. $ a$ by$n-1 $ times yields

$ \displaystyle \int_{0}^{\infty} \frac{(-1)^{r-1}(n-1) ! x^{r}}{\left(x^{m}+a\right)^{n}}=\frac{\pi}{m} \csc \frac{(r+1) \pi}{m}(-1)^{n-1}\left(1-\frac{r+1}{m}\right)\left(2-\frac{r+1}{m}\right)\cdots\left(n-1-\frac{r+1}{m}\right)a^{-\left(n-\frac{r+1}{m}\right)} \tag*{}$

Now we can conclude that

$$\boxed{ \displaystyle \int_{0}^{\infty} \frac{x^{r} d x}{\left(x^{m}+a\right)^{n}}=\frac{\pi}{m(n-1) !} \csc \frac{(r+1) \pi}{m} \prod_{j=1}^{n-1}\left(j-\frac{r+1}{m}\right) a^{-\left(n-\frac{r+1}{m}\right)}}$$

Come back to our integral and use Binomial Theorem.

\begin{aligned}\int_{0}^{\infty} \frac{1}{\left(x^{4}-x^{2}+1\right)^{n}} d x =& \int_{0}^{\infty}\left(\frac{x^{2}+1}{x^{6}+1}\right)^{n} d x =\sum_{k=0}^{n}\left(\begin{array}{l}n \\k\end{array}\right) \int_{0}^{\infty} \frac{x^{2 k}}{\left(x^{6}+1\right)^{n}} d x\end{aligned}

Putting $ m=6, r=2k $ and $ a=1 $ yields

$\displaystyle \boxed{\begin{aligned}\int_{0}^{\infty} \frac{1}{\left(x^{4}-x^{2}+1\right)^{n}} d x &=\sum_{k=0}^{n}\left(\begin{array}{c}n \\k\end{array}\right) \frac{\pi}{6(n-1) !} \csc \frac{(2k+1) \pi}{6} \prod_{j=1}^{n-1}\left(j-\frac{2k+1}{6}\right)\\&=\sum_{k=0}^{n}\left[\left(\begin{array}{l}n \\k\end{array}\right) \frac{\pi}{6} \csc \frac{(2 k+1) \pi}{6} \prod_{j=1}^{n-1}\left(1-\frac{2 k+1}{6 j}\right)\right] \end{aligned} }\tag*{} $

Answer 7

With CAS and help of MellinTransfrom:

$$\int_0^{\infty } \frac{1}{\left(x^4-x^2+1\right)^n} \, dx=\\\mathcal{M}_s^{-1}\left[\int_0^{\infty } \mathcal{M}_a\left[\frac{1}{\left(x^4-a x^2+1\right)^n}\right](s) \, dx\right](1)=\\\mathcal{M}_s^{-1}\left[\int_0^{\infty } \frac{(-1)^{-s} x^{-2 s} \left(1+x^4\right)^{-n+s} \Gamma (n-s) \Gamma (s)}{\Gamma (n)} \, dx\right](1)=\\\mathcal{M}_s^{-1}\left[\frac{(-1)^{-s} \Gamma \left(\frac{1}{4}-\frac{s}{2}\right) \Gamma \left(-\frac{1}{4}+n-\frac{s}{2}\right) \Gamma (s)}{4 \Gamma (n)}\right](1)=\\\frac{\Gamma \left(\frac{1}{4}\right) \Gamma \left(-\frac{1}{4}+n\right) \, _2F_1\left(\frac{1}{4},-\frac{1}{4}+n;\frac{1}{2};\frac{1}{4}\right)}{4 \Gamma (n)}+\frac{\Gamma \left(\frac{3}{4}\right) \Gamma \left(\frac{1}{4}+n\right) \, _2F_1\left(\frac{3}{4},\frac{1}{4}+n;\frac{3}{2};\frac{1}{4}\right)}{4 \Gamma (n)}$$