Solution (source code)

= Solution

Let $r=\Delta\tau/h^2$ and impose periodic indices modulo $J$. With $\boldsymbol\Phi^n=(\Phi_0^n,\ldots,\Phi_{J-1}^n)^T$, the <linear system> is
$$
\boxed{A\boldsymbol\Phi^{n+1}=\boldsymbol\Phi^n-4\pi G\Delta\tau\boldsymbol\rho,\quad
A=(1+2r)I-r(S+S^{-1}),}
$$
where $S$ shifts one grid cell cyclically. For $J\ge3$, $A$ has diagonal $1+2r$, adjacent and corner entries $-r$, and all other entries zero. It is a cyclic tridiagonal <matrix>, rather than a nonperiodic tridiagonal one.

The periodic Laplacian annihilates the constant <vector>, so its steady equation requires $\boxed{\sum_j\rho_j=0}$. A strictly positive gravitational <density> violates this condition: summing the update gives $\bar\Phi^{n+1}=\bar\Phi^n-4\pi G\Delta\tau\bar\rho$, an indefinitely drifting mean. Periodic gravitational calculations normally solve for $\rho-\bar\rho$ and fix the arbitrary mean <gravitational potential>. The per-step <matrix> $A$ itself remains invertible; the issue is existence of a periodic steady solution of <Poisson's equation>.