The function to be integrated (and any related function that will ultimately solve the problem) hates and wants to avoid the points $0,1$, so the pink contour would be a first guess to apply complex analysis (C.A. for short below):
Let $a$ be $a=2\pi i$.
We quickly find a polynomial $P=P_a$, so that $P(X+a)-P(X)=X^2$, this is
$$
P(X) = \frac 1{3a}\left(x^3-\frac 32ax^2+\frac 12a^2x\right)
\ .
$$
Then
$$
\int_{\color{magenta}{\text{Pink contour}}}
\frac {P(\ln z)}{(z-1)^2}
=2\pi\; i\sum_{r\text{ Residue}}\operatorname{Res}_{z=r}
\frac {P(\ln z)}{(z-1)^2}
=0
\ .
$$
Now we let $R$ go to $+\infty$.
The upper blue segment produces an integral of the shape $$\displaystyle
\int_{0+\varepsilon}^\infty
\frac{P(\ln x)}{(x-1)^2}\; dx\ ,
$$
then the semicircle does not contribute, since for $z=R e^{it}$ the numerator contributes with $\ln^3 R+O(1)$, the denominator is $O(R^{-2})$, the whole fraction is $\sim R^{-2}\ln R$, and the contour is in $O(R)$,
and finally the lower pink segment produces an integral of the shape
$$\displaystyle
\int_\infty^{0+\varepsilon}
\frac{P(\ln x+2\pi i)}{(x-1)^2}\; dx\ ,
$$
if we succeed to push it beyond the holes.
Here i have to appologize for using $\ln$ for the real function $\log$, and also for the used branch of the logarithm, which is the reason for getting $\ln(x+2\pi i)$ by the monodromy of the logarithm when the contour comes back.
By the choice of $P$, the two integrals come together to build
(up to sign)
$$
\int_{0+\varepsilon}^\infty
\frac{P(\ln x+2\pi i)-P(\ln x)}{(x-1)^2}\; dx
=
\int_{0+\varepsilon}^\infty
\frac{x^2}{(x-1)^2}\; dx
\ ,
$$
and this value comes then during "pushing" exactly from the holes in $0,1$. Let us compute the two contributions. Let $C(z_0,\varepsilon)$ the contour / the circle centered in $z_0\in\Bbb C$ with radius $\varepsilon>0$. Then
using the substitution $z=z_0+\varepsilon e^{it}$, $t\in[0,2\pi]$,
$$
\left|
\int_{C(0,\varepsilon)}
\frac
{\ln^k z}
{(z-1)^2}\; dz
\right|
\sim 2\pi\epsilon\cdot |\ln\varepsilon+O(1)|^k
\in O(\varepsilon^{1/2})
$$
produces no contribution, and
$$
\begin{aligned}
\int_{C(1,\varepsilon)}
\frac
{\ln^k z}
{(z-1)^2}\; dz
&=
\int_0^{2\pi}
\frac
{\displaystyle\ln^k(1+\varepsilon e^{it})}
{\displaystyle(1+\varepsilon e^{it}-1)^2}
\; i\varepsilon e^{it}\; dt
\\
&=
\int_0^{2\pi}
\frac
{\displaystyle
\left(
\frac 11\varepsilon e^{it}
-\frac 12\varepsilon^2 e^{2it}
+\frac 13\varepsilon^3 e^{3it}
\pm\dots
\right)^k}
{\displaystyle \varepsilon^2 e^{2it}}
\; i\varepsilon e^{it}\; dt
\\
&\to
\begin{cases}
2\pi i = a&\text{ for }k=1\ ,\\
0&\text{ for }k\ge 2\ ,
\end{cases}
\qquad
\text{ when }
\varepsilon\to0\ .
\end{aligned}
$$
The coefficient of $x$ in $P$ was $\frac a6$, so we get the contribution from $1$ in the form (well, computations were done up to sign)
$$
\pm \frac 16a^2=\pm\frac 16(2\pi i)^2=\mp \color{blue}{\frac23\pi^2}\ .
$$
(We are choosing the plus sign of course, the integral is built w.r.t. a positive function.)
$\color{red}\square$
Job done, but...
This method applies mot-a-mot to the integrals
$$
\int_0^\infty
\frac{\ln^4 x}{(x-1)^2}\; dx\ ,\qquad
\int_0^\infty
\frac{\ln^6 x}{(x-1)^2}\; dx\ ,\qquad
\int_0^\infty
\frac{\ln^8 x}{(x-1)^2}\; dx\ ,\qquad
\dots
$$
in the following way:
- we build the corresponding Bernoulli polynomial $Q$ so that $Q(x+1)-Q(x)=(2n+1)x^{2n}$. An $a$-homogenized, $1/(2n+1)$-normed version of it delivers $P$ with $P(0)=0$, $P(x+a)-P(x)=x^{2n}$, and store the coefficient $c_1$ of $P(x)$ in $x^1$,
- we use the same contour of integration, the same argument applies,
- there is no contribution from the singularity in zero,
- the singularity in one contributes with $a=2\pi i$ again.
Numerical experiments.
Numerically, i prefer to compute the half value,
$$
\begin{aligned}
J(2n)
&= \int_1^\infty \frac{\ln^{2n} x}{(x-1)^2}\; dx
\qquad\text{ Substitution: }
y=\frac{1}{x},\ x=\frac{1}{y},\ dx=-\frac 1{y^2}\; dy
\\
&= \int_0^1 \frac{\ln^{2n} y}{(1-y)^2}\; dy
=\frac 12\int_0^\infty \frac{\ln^{2n} x}{(x-1)^2}\; dx
\ .
\end{aligned}
$$
Then the expected value of the integral is numerically validated:
sage: for n in [2, 4, 6, 8, 10]:
....: Q = bernoulli_polynomial(x, n+1)/(n+1)
....: c1 = Q.coefficient( x^1 )
....: print "n = %s" % n
....: print "B(x, %s) / %s is %s" % (n+1, n+1, Q)
....: print "Its coefficient in x^1 is %s" % c1
....: f(x) = log(x)^n / (x-1)^2
....: myintegral = numerical_integral( lambda x: f(x), (0,1) )[0]
....: print "f(x) is %s" % (f(x))
....: print "Integral of f on (0, 1) is ~ %f" % myintegral
....: print "| c1 * (2 pi)^%2s / 2 | is ~ %f" % (n, (abs(c1)*(2*pi)^n/2).n())
....: print
....:
n = 2
B(x, 3) / 3 is 1/3*x^3 - 1/2*x^2 + 1/6*x
Its coefficient in x^1 is 1/6
f(x) is log(x)^2/(x - 1)^2
Integral of f on (0, 1) is ~ 3.289868
| c1 * (2 pi)^ 2 / 2 | is ~ 3.289868
n = 4
B(x, 5) / 5 is 1/5*x^5 - 1/2*x^4 + 1/3*x^3 - 1/30*x
Its coefficient in x^1 is -1/30
f(x) is log(x)^4/(x - 1)^2
Integral of f on (0, 1) is ~ 25.975758
| c1 * (2 pi)^ 4 / 2 | is ~ 25.975758
n = 6
B(x, 7) / 7 is 1/7*x^7 - 1/2*x^6 + 1/2*x^5 - 1/6*x^3 + 1/42*x
Its coefficient in x^1 is 1/42
f(x) is log(x)^6/(x - 1)^2
Integral of f on (0, 1) is ~ 732.487004
| c1 * (2 pi)^ 6 / 2 | is ~ 732.487005
n = 8
B(x, 9) / 9 is 1/9*x^9 - 1/2*x^8 + 2/3*x^7 - 7/15*x^5 + 2/9*x^3 - 1/30*x
Its coefficient in x^1 is -1/30
f(x) is log(x)^8/(x - 1)^2
Integral of f on (0, 1) is ~ 40484.398950
| c1 * (2 pi)^ 8 / 2 | is ~ 40484.399002
n = 10
B(x, 11) / 11 is 1/11*x^11 - 1/2*x^10 + 5/6*x^9 - x^7 + x^5 - 1/2*x^3 + 5/66*x
Its coefficient in x^1 is 5/66
f(x) is log(x)^10/(x - 1)^2
Integral of f on (0, 1) is ~ 3632409.110395
| c1 * (2 pi)^10 / 2 | is ~ 3632409.114224
sage:
Also for some two bigger values:
sage: for n in [20, 30 ]:
....: Q = bernoulli_polynomial(x, n+1)/(n+1)
....: c1 = Q.coefficient( x^1 )
....: print "n = %s" % n
....: print "c1 is %s" % QQ(c1).factor()
....: f(x) = log(x)^n / (x-1)^2
....: myintegral = numerical_integral( lambda x: f(x), (0, 1), eps_abs=1e-4 )[0]
....: print "f(x) is %s" % (f(x))
....: print "Integral of f on (0, 1) is ~ %r" % myintegral
....: print "| c1 * (2 pi)^%2s / 2 | is ~ %r" % (n, (abs(c1)*(2*pi)^n/2).n())
....: print
....:
n = 20
c1 is -1 * 2^-1 * 3^-1 * 5^-1 * 11^-1 * 283 * 617
f(x) is log(x)^20/(x - 1)^2
Integral of f on (0, 1) is ~ 2.432904323980626e+18
| c1 * (2 pi)^20 / 2 | is ~ 2.43290432907279e18
n = 30
c1 is 2^-1 * 3^-1 * 5 * 7^-1 * 11^-1 * 31^-1 * 1721 * 1001259881
f(x) is log(x)^30/(x - 1)^2
Integral of f on (0, 1) is ~ 2.6525284549431265e+32
| c1 * (2 pi)^30 / 2 | is ~ 2.65252860059228e32
sage:
Best Answer
If you choose to evaluate this using a contour it will not just be similar, but identical to the contour for $\cot(\pi z)$ with a variable change.
It is literally identical, with only a variable change, and you can also choose to do this variable change before or after the integration itself. What you are asking for in the question is doing the variable change before the contour integration, what follows below is doing the variable change after:
If you want to do the variable change before the contour integration, just look at $z\rightarrow \frac{iz}{2}$. Then we are integrating the function $\pi z \coth \left( \frac{\pi z}{2}\right)$, which really the same function as $\pi z\cot\left(\pi z\right)$.