How to calculate $\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx$?

Question

Using the derivative of beta function, find

$$I=\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx$$

setting $x^2=y$ gives

$$I=\frac18\int_0^1\frac{\ln^2y\ln^2(1-y)}{\sqrt{y}(1-y)}\ dy=\frac18\left.\frac{\partial^4}{\partial a^2\partial b^2}\text{B}(a,b)\right|_{a\mapsto 1/2\\b\mapsto0^{+}}$$

Any good software that can find the 4th derivative and also gives the final result? Wolfram fails to calculate it (or maybe I do not know how to use it well) and when I tried to do it manually, some terms involve $\psi(b)$ and if we take the limit, then $\psi(0)$ is undefined and even if I take the limit of $\psi(b)$ together with other terms, still undefined. I do not know how to avoid this problem as I am not experienced with the beta function.

Thank you.

---

Note: Solution should be done without using harmonic series.

Answer 1

All the following Mathematica commands calculate your limit, in decreasing order of time (the more naive one uses more time):

Limit[D[Gamma[a]*Gamma[b]/Gamma[a + b], {a, 2}, {b, 2}] /. {a -> 
  1/2 + x, b -> x}, x -> 0] // FunctionExpand // Expand

the above command directly calculates the limit, by choosing a path approaching $(a,b)=(1/2,0)$. It takes $32$ seconds on my machine.

D[Normal[Series[
        Gamma[a]*Gamma[b]/Gamma[a + b], {a, 1/2, 4}, {b, 0, 4}]], {a, 
       2}, {b, 2}] /. a -> 1/2 /. b -> 0 // FullSimplify // 
  Expand 

rather than calculating the limit, this one uses series expansion up to constant term. It takes $12$ seconds.

Normal[Series[
     D[Gamma[a]*Gamma[b]/Gamma[a + b], {a, 2}, {b, 2}], {a, 1/2, 
      0}, {b, 0, 0}]] // FullSimplify // Expand

this one does not even calculate derivatives, instead uses series expansion up to 4th order. It takes $3.5$ seconds.

D[Exp[Series[
        LogGamma[a] + LogGamma[b] - LogGamma[a + b], {a, 1/2, 4}, {b, 
         0, 4}]], {a, 2}, {b, 2}] /. a -> 1/2 /. b -> 0 // 
   FullSimplify // Expand

This use the well-known simple series of log gamma function, it takes only $0.5$ seconds.

---

It's easy to guess why the fourth one is most efficient. To see how much is used for each computation, execute ClearSystemCache[];(your command)//Timing.

Such beta limit arising from logarithm integrals is well-known, it's also not difficult to write down an recursion for it.

Answer 2

A magnificent solution by Cornel without using the derivative of beta function.

---

We have the identity

$$\frac{\ln^2(1-x)}{1-x}=\sum_{n=1}^\infty (H_n^2-H_n^{(2)})x^n$$

replace $x$ with $x^2$, then multiply both sides by $\ln^2x$ and integrate from $x=0$ to $1$ we get

$$\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx=2\sum_{n=1}^\infty\frac{H_n^2-H_n^{(2)}}{(2n+1)^3}\tag1$$

---

By the master theorem, we have

$$3n\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(k+1)(k+n+1)}=H_n^3+3H_nH_n^{(2)}+2H_n^{(3)}$$

differentiate both sides with respect to $n$ to get

$$3\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(k+n+1)^2}=\frac{d}{dn}\left(H_n^3+3H_nH_n^{(2)}+2H_n^{(3)}\right)$$

differentiate again and let $n\mapsto -1/2$ we get

$$-6\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(k+1/2)^3}=\frac{d^2}{dn^2}\left(H_n^3+3H_nH_n^{(2)}+2H_n^{(3)}\right)_{n\mapsto-1/2}$$

or

$$\sum_{k=1}^\infty \frac{H_k^2-H_k^{(2)}}{(2k+1)^3}=\frac{31}{2}\zeta(5)-\frac{45}{8}\ln2\zeta(4)+\frac72\ln^22\zeta(3)-7\zeta(2)\zeta(3)\tag2$$

---

From (1) and (2) we get

$$\int_0^1\frac{\ln^2x\ln^2(1-x^2)}{1-x^2}\ dx=31\zeta(5)-\frac{45}{4}\ln2\zeta(4)+7\ln^22\zeta(3)-14\zeta(2)\zeta(3)$$

Answer 3

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}[1]{\mathcal{#1}} \newcommand{\mrm}[1]{\mathrm{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} I & \equiv {1 \over 8}\int_{0}^{1}{\ln^{2}\pars{y}\ln^{2}\pars{1 - y} \over \root{y}\pars{1 - y}}\,\dd y \\[5mm] & = \left.{1 \over 8}{\partial^{2} \over \partial \mu^{2}\,\partial \nu^{2}}\int_{0}^{1} \pars{y^{\,\mu - 1/2} - 1}\pars{1 - y}^{\,\nu - 1} \,\dd y\, \right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}} \\[5mm] & = \left.{1 \over 8}{\partial^{4} \over \partial \mu^{2}\,\partial \nu^{2}}\bracks{{\Gamma\pars{\mu +1/2} \Gamma\pars{\nu} \over \Gamma\pars{\mu + \nu + 1/2}} - {1 \over \nu}}\, \right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}} \\[5mm] & = \left.{1 \over 8}{\partial^{4} \over \partial \mu^{2}\,\partial \nu^{2}}\braces{{1 \over \nu}\bracks{{\Gamma\pars{\mu +1/2} \Gamma\pars{\nu + 1} \over \Gamma\pars{\mu + \nu + 1/2}} - 1}}\, \right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}} \\[5mm] & = \left.{1 \over 24}{\partial^{5} \over \partial \mu^{2}\,\partial \nu^{3}}{\Gamma\pars{\mu +1/2} \Gamma\pars{\nu + 1} \over \Gamma\pars{\mu + \nu + 1/2}}\, \right\vert_{\ {\Large\mu\ =\ 0} \atop {\Large\nu\ =\ 0}} \\[5mm] & = \bbx{-\,{7 \pi ^2 \over 3}\,\zeta\pars{3} + 31\zeta\pars{5} + 7\zeta\pars{3}\ln^{2}\pars{2} - {1 \over 8}\,\pi^{4}\ln\pars{2}} \\[5mm] & \approx 0.0654 \end{align}