How prove this limit $\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{i+j}{i^2+j^2}=\frac{\pi}{2}+\ln{2}$

show that: this limit $$I=\lim_{n\to\infty}\dfrac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n}\dfrac{i+j}{i^2+j^2}=\dfrac{\pi}{2}+\ln{2}$$

My try: $$I=\lim_{n\to\infty}\dfrac{1}{n^2}\sum_{i=1}^{n}\sum_{j=1}^{n}\dfrac{\dfrac{i}{n}+\dfrac{j}{n}}{\left(\dfrac{i}{n}\right)^2+\left(\dfrac{j}{n}\right)^2}=\int_{0}^{1}\int_{0}^{1}\dfrac{x+y}{x^2+y^2}dxdy?$$

This idea is true? and have other methods?

Evaluation of the integral(after Thomas Andrews's calculation)

$$ \int_0^1 \log(1+y^2) dy= y\log(1+y^2)|_0^1 -\int_0^1 \frac{2y^2}{1+y^2} dy$$ $$=\log 2-\int_0^1 \frac{2y^2+2-2}{1+y^2} dy$$ $$=\log 2 - 2 + \frac{\pi}{2} $$

$$-\int_0^1 \log (y^2) dy = -\int_0^1 2\log y dy = 2\int_0^{\infty} e^{-x} dx = 2$$

Hence, the integral all together is $$\log 2 + \frac{\pi}{2}$$

Application of Stolz-Cesaro:

If we can prove that the limit exists $$ \lim_{n\rightarrow\infty}\left( \sum_{i=1}^n \frac{i+n}{i^2+n^2} +\sum_{j=1}^{n-1} \frac{n+j}{n^2+j^2}\right) $$ then the answer is the same as the limit.

Evaluation of this limit requires integral, but this time not a 'double integral'.

$$\lim_{n\rightarrow\infty} \sum_{i=1}^n \frac{i+n}{i^2+n^2}$$ $$=\lim_{n\rightarrow\infty} \sum_{i=1}^n \frac{ \frac{i}{n}+1}{\frac{i^2}{n^2}+1} \frac{1}{n}$$ $$=\int_0^1 \frac{x+1}{x^2+1} dx$$ $$=\frac{\pi}{4}+\frac{1}{2}\log 2.$$ The other part can be dealt similarly, it has the same value as above. Thus, all together the answer becomes $$\frac{\pi}{2}+\log 2.$$

[Incomplete answer]

Probably a good place to start is to note: $$\int_0^1\int_0^1\frac{x}{x^2+y^2}\,dx\,dy = \int_0^1\int_0^1\frac{y}{x^2+y^2}\,dx\,dy$$

by symmetry, assuming the integral exists.

So you really need to compute $$\int_0^1\int_0^1 \frac{2x}{x^2+y^2}\,dx \,dy$$

That looks easier to compute.

$$\int_0^1\frac{2x}{x^2+y^2} dx = \log (1+y^2)-\log(y^2) = \log(1+y^{-2})$$

That leads to $$\int_0^1 \log(1+y^{-2})\,dy$$

Not sure if that's going to do any good.

(1) Integral is finite.

The function in the integral is $\leq \frac{2}{x+y}$. The integral of $1/(x+y)$ on the square is bounded, because $\frac{1}{x+y} \in [N,N+1)$ for points in the square with $x+y \in [1/(N+1), 1/N]$ is (for large $N$) a trapezoid with height approximately $1/N^2$ and width of approximately $1/N$, for a total contribution to the integral of about $N \times \frac{1}{N^3} (1 + o(1)) \sim \frac{1}{N^2} $. This converges when summed on integer $N \geq n_0$.

(2) Sum converges to the integral.

The sum is a "minimum" Riemann sum for the integral on a partition of the square into $n^2$ equal small. The "maximum" Riemann sum on the same squares except the one containing $(0,0)$ (whose contribution to the sum and to the integral both converge to $0$, and can be ignored) is the same sum for $i,j < n$, plus a small modification along the boundary of the square ($ij=0$ or $i=n$ or $j=n$). These changes to the sum are essentially line integrals of a the function along the sides of the square, weighted by $1/n$, and the modifications go to $0$.

(3) Evaluation of integral.

This can be done using polar coordinates, and maybe by easier methods. I did not work it out, since parts 1+2 already identified what the limit "is" and showed convergence.