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≡∫10∫10lnx−lnyx−ydxdy.
Ya que f(x,y)=lnx−lnyx−y=f(y,x), la superficie f(x,y) es simétrica respecto al plano bisector x=y; entonces,
J2=∫10dx∫x0lnx−lnyx−ydy.
Con el cambio de variabley≡tx, t∈(0, 1), la integral ∫x0lnx−lnyx−ydy, se transforma en la siguiente que no depende de x, I≡−∫10lnt1−tdt.
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≡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 ζ(2), cuyo valor es bien conocido:
I=∫01ln(1−s)sds=Li2(1)=ζ(2)=π26. Por lo tanto, la solución a la integral propuesta es J=π23.
Nota. Lo intenté con transformaciones polilogarítmicas, pero no pude obtener el resultado π23.