Approximation with error term / bounds for double summation

Question

I'm looking for ways to find continuous functions that approximate double summations of the form $S(n)=\sum _{j=1}^n \sum _{i=1}^n f(n-i j)$ for functions $f:\mathbb{R}\to \mathbb{R}$, and for large $n$. Take, for example,

$$f(n):=\sum _{j=1}^n \sum _{i=1}^n \frac{1}{(n-i j)^4+8}$$

Heuristically, the harmonic function $h(n)=H_{\sqrt{n}}-\sqrt{\gamma }$ (and its implied analytic continuation) would seem to be a pretty good approximation to $f$ (although of course, without upper and lower bounds, its value would be limited):

But this only a rough heuristic, and it could be wrong. I am not trying to prove this specifically. Rather:

How would I go about finding an approximation algebraically?

I tried using integrals to approximate $f$ but this led (via Mathematica) to an amazingly hairy expression with around $100$ terms. Presumably this is because the behaviour of the function is somewhat erratic, but the upshot is that it's not very useful.

So, how do I find something usable (ideally with bounds, but all help gratefully received)?

UPDATE:

I was wrong to dismiss the double integral so quickly - the expression may be hairy, but the result is better than I anticipated:

Which leads to an

UPDATED QUESTION:

1. Is approximation by integration the best approach? I have looked into the Euler-MacLaurin summation formula, but the double-integration defeats me.

2. If no, what approach should I use?

3. If integration is the way to go, then given the double-integral approximation

$$\sum _{j=1}^n \sum _{i=1}^n \frac{1}{(n-i j)^4+8}\approx\int _1^n \int _1^n \frac{1}{(n-t u)^4+8} \mathrm dt \mathrm du$$

and the implied error term

$$\sum _{j=1}^n \sum _{i=1}^n \frac{1}{(n-i j)^4+8}=\int _1^n \int _1^n \frac{1}{(n-t u)^4+8} \mathrm dt \mathrm du + \epsilon_n$$

Then how do I go about finding the error term $\epsilon_n$ as a function of $n$?

(For those who want to see the full double integral expression, the Mathematica code is below.)

Integrate[Integrate[1/((n - t*u)^4 + 8), {t, 1, n}, 
   Assumptions -> Element[n, Reals] && n > 1 && Element[t, Reals] && t >= 1 && 
    Element[u, Reals] && u >= 1], {u, 1, n}, 
  Assumptions -> Element[n, Reals] && n > 1]

Answer 1

This is not an answer. It is more heuristics, to explain why $f(n)$ behaves like $\log(n)$. I hope it can be formalized.

First, by changing the order of summation :

$$f(n) = \sum_{k = 1}^{n^2} Card\{(i, j) | ij = k\}\frac1{(n - k)^4 +8}$$

where $(i,j)$ is ordered.

$Card\{(i, j) | ij = k\}$ is of order $d(k)$, the number of divisors of $k$. But : $$\frac1{N} \sum_{k=1}^{N} d(k) = \log(N) + O(1)$$ (https://en.wikipedia.org/wiki/Divisor_summatory_function)

It has to be formalized but, we may consider that, $d(k) \simeq \log(k)$. Thus, maybe we can prove that :

$$f(n) \simeq \sum_{k = 1}^{n^2} \frac{\log(k)}{(n - k)^4 +8}$$

This sum is of order $\log(n)$, by spliting $k < n; k = n; k > n$.