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:
Is approximation by integration the best approach? I have looked into the Euler-MacLaurin summation formula, but the double-integration defeats me.
If no, what approach should I use?
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]