In a flat toric universe (up connects down, right connects left and front connects back), every points repeats at $size_x$, $size_y$ and $size_z$ intervals.
In such case the Newtonian gravitational acceleration undergone in one point x, y, from a single mass M placed in one corner of a cube, would be undergone from every «image» of this mass M.
eg: lets take one cube of space with a mass (the red dot) in a corner :
Then let's repeat it ad-infinitum in every direction, we see that the mass appears in all corners of every cubes:
In such a situation, the acceleration undergone in the center of any cube would be null and it's general formula would be: $a(x, y, z) = - M G \lim\limits_{N\to\infty}\sum\limits_{i=-N}^{N}\sum\limits_{j=-N}^{N}\sum\limits_{k=-N}^{N} {(x-size_x.i, y-size_y.i, z-size_z.i) \over |{(x-size_x.i, y-size_y.i, z-size_z.i)}|^3}$
To make computing simpler, it would be nice to have a simple equivalent of this formula. which is obviously a periodic function.
In a first approach, I should be able to bring back this problem to a simpler one dimensional function : $\begin{align}f(x) &= \lim\limits_{N\to\infty}\sum\limits_{n=-N}^{N}{x-n\pi \over ((x-n\pi)^2+a^2)^{3/2}} \\\\ &={x \over (x^2+a^2)^{3/2}} + \lim\limits_{N\to\infty}[\sum\limits_{n=-N}^{-1}{x-n\pi \over ((x-n\pi)^2+a^2)^{3/2}} + \sum\limits_{n=1}^{N}{x-n\pi \over ((x-n\pi)^2+a^2)^{3/2}}] \\\\ &={x \over (x^2+a^2)^{3/2}} + \sum\limits_{n=1}^{\infty}[{x+n\pi \over ((x+n\pi)^2+a^2)^{3/2}} + {x-n\pi \over ((x-n\pi)^2+a^2)^{3/2}}]\end{align}$
But then, I wonder how I could compute f(x)...
As told above, it will obviously be a periodic function, but it doesn't look like an obvious one. On this simulation I made, one can see various graphs for this function (red for a = 0.6, green for a = 0.3 and blue for a = 0). In addition, cyan is a graph of $-\cot x$, magenta is a graph of $-\cot x/f(x)$ and yellow is $|\sin x|$ (The «V» shape on the magenta graph is caused by the 0/0 situation).
A possible solution might be to take the Fourier transform and multiply it by a dirac comb (Ш function) and take the reverse Fourier transform - or said otherwise, to take the value corresponding with ω being an integer multiple of $2 \pi$. I tried it online with the WolframAlpha website but the result are quite obfuscated :
$\begin{align}ℱ_x[{x \over (x^2 + a^2)^{3/2}}](ω) = -ⅈ [&\sqrt{2 \over π} \sqrt{1 \over a^2} K_1(\sqrt{a^2} ω sgn(ω)) Abs'(ω)\\\\ &+ {1 \over \sqrt{2π}}\sqrt{1 \over a^2} |ω| (\sqrt{a^2} ω sgn'(ω) + \sqrt{a^2} sgn(ω)) (-K_0(\sqrt{a^2} ω sgn(ω)) - K_2(\sqrt{a^2} ω sgn(ω)))]\end{align}$
I don't know if it is due to my own misunderstanding, but it seems that many of those notations might be quite simplified. What are $Abs'(⬚)$ and $Sgn'(⬚)$ ? Aren't they just the same as respectively $Sgn(⬚)$ and $2\delta(⬚)$ ? Why do they mix «$|⬚|$», «$⬚ sgn(⬚)$» and «$\sqrt{⬚^2}$» ? Isn't it all the same thing ? If I follow this, I should rewrite this formula as:
$\begin{align}ℱ_x[{x \over (x^2 + a^2)^{3/2}}](ω) &=-ⅈ \left[\sqrt{2 \over π} {1 \over |a|} K_1(|a ω|) sgn(ω) + {1 \over \sqrt{2π}}{1 \over |a|} |ω| (|a| ω 2 \delta(ω) + |a| sgn(ω)) (-K_0(|a ω|) - K_2(|a ω|))\right]\\\\ &= {ⅈ \over |a|}\sqrt{2 \over π}[|aω| (ω \delta(ω) + sgn(ω)/2) (K_0(|aω|) + K_2(|aω|)) - sgn(ω) K_1(|a ω|) ]\\\\ &\text{(for ω=0, ωδ(ω) is undetermined but since I'm looking for a symmetric} \\ &\ \ \text{solution, I will reject a continuous component by setting ωδ(ω)=0)} \\\\ &= {ⅈ \over |a|}\sqrt{2 \over π}[{|a|ω \over 2} (K_0(|aω|) + K_2(|aω|)) - sgn(ω) K_1(|a ω|) ]\end{align}$
Then I might take the values of my transform for all ω multiple of 2pi and use them in a Fourier serie.