Método más rápido para calcular la solución exacta de la siguiente integral (basado en las ideas de Fedor Petrov: https://mathoverflow.net/users/4312/fedor-petrov, Doble integral con logaritmos, URL (versión: 2019-04-15): https://mathoverflow.net/q/328126):
$$J\equiv \int_{0}^{1}{\int_{0}^{1}{\frac{\ln x-\ln y}{x-y}}}dxdy .$$
Ya que $$f\left( x,y \right)=\frac{\ln x-\ln y}{x-y}=f\left( y,x \right),$$ la superficie $f\left( x,y \right) $ es simétrica respecto al plano bisector $x = y$; entonces,
$$\frac{J}{2}=\int_{0}^{1}{dx\int_{0}^{x}{\frac{\ln x-\ln y}{x-y}}}dy.$$
Con el cambio de variable$$y\equiv tx,\ t\in \left( 0,\ 1 \right),$$ la integral $$\int_{0}^{x}{\frac{\ln x-\ln y}{x-y}}dy,$$ se transforma en la siguiente que no depende de $x$, $$ I\equiv -\int_{0}^{1}{\frac{\ln t}{1-t}\,}dt.$$
La integración sobre el cuadrado unidad $(0, 0), (1, 0), (1, 1), (0,1)$ se reduce a la integración sobre el triángulo $(0, 0), (1, 0), (1, 1).$
Para resolver la integral $I$ realizaremos el nuevo cambio de variable, $$s\equiv 1-t,$$ por el cual $I$ se transforma en la integral que define el dilogaritmo, cuyo valor para $s = 1$ coincide con la Función zeta de Riemann $\zeta \left( 2 \right)$, cuyo valor es bien conocido:
$$I=\int_{1}^{0}{\frac{\ln \left( 1-s \right)}{s}\,}ds=\text{L}{{\text{i}}_{2}}\left( 1 \right)=\zeta \left( 2 \right)=\frac{{{\pi }^{2}}}{6}.$$ Por lo tanto, la solución a la integral propuesta es $$J=\frac{{{\pi }^{2}}}{3}.$$
Nota. Lo intenté con transformaciones polilogarítmicas, pero no pude obtener el resultado $\frac{{{\pi }^{2}}}{3}.$
