(Actualización: Algunas correcciones se han agregado ya he resuelto el problema con $a_{ij}=\frac{2x_ix_j}{x_i+x_j}$ en lugar de $a_{ij}=\frac{2x_ix_j}{x_i^2+x_j^2}$. Gracias a Paata Ivanisvili por su comentario.)
El truco es escribir la expresión de la $\frac{1}{x^2+y^2}$ como una integral.
Para cada vector distinto de cero $(u_1,\ldots,u_n)$ de reales,
$$
(u_1,\ldots,u_n) \begin{pmatrix}
a_{11} & \dots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \dots & a_{nn} \\
\end{pmatrix}
\begin{pmatrix} u_1\\ \vdots \\ u_n\end{pmatrix} =
\sum_{i=1}^n \sum_{j=1}^n a_{ij} u_i u_j =
$$
$$
=2 \sum_{i=1}^n \sum_{j=1}^n u_i u_j x_i x_j \frac1{x_i^2+x_j^2} =
2\sum_{i=1}^n \sum_{j=1}^n u_i u_j x_i x_j \int_{t=0}^\infty
\exp \Big(-(x_i^2+x_j^2)t\Big) \mathrm{d}t =
$$
$$
=2\int_{t=0}^\infty
\left(\sum_{i=1}^n u_i x_i \exp\big(-{x_i^2}t\big)\right)
\left(\sum_{j=1}^n u_j x_j \exp\big(-{x_j^2}t\big)\right)
\mathrm{d}t =
$$
$$
=2\int_{t=0}^\infty
\left(\sum_{i=1}^n u_i x_i \exp\big(-{x_i^2}\big)\right)^2
\mathrm{d}t \ge 0.
$$
Si $|x_1|,\ldots,|x_n|$ son distintas y distinto de cero, entonces el último integrando no puede ser constante cero: para un gran $t$ el mínimo de $|x_i|$ $u_i\ne0$ determina el orden de magnitud. De modo que la integral es estrictamente positivo.
Si hay valores iguales entre $|x_1|,\ldots,|x_n|$, entonces la integral puede desaparecer. En consecuencia, las filas correspondientes de la matriz son iguales o el negativo de la otra, por lo que la matriz es sólo positivo semidefinite.
Usted puede encontrar muchas variantes de esta desigualdad. Ver
problema A. 477. de la KöMaL de la revista.
La entrada $a_{ij}$ puede ser sustituido por $\left(\frac{2x_ix_j}{x_i+x_j}\right)^c$ con un fijo arbitrarios positivo $c$; véase la
problema A. 493. de KöMaL.
