So, I've been attempting to design a simple solver for a problem of finding the gravitational potential of a system using Poisson's equation (let's call the potential phi, $\phi$). The goal is that I apply this potential to the system, move it forward one time step, and then repeat.

Right now, I am working with a periodic system (so assuming the system operates while being bound within a circle), and thus need to use periodic boundary conditions. Up to this point, I have been using this post to develop my method with Finite Differences.

Here's the issue. Let's say that the system is represented by a wave in the direct "center" of the circle, then the potential should by pseudo-linear from zero until the center where it drops off and returns pseudo-linearly to zero. enter image description here

Now, if I was to shift this wave halfway across the circle and do this again, we should see the same potential also shifted halfway across the circle. In this case, the "center" would be 0, and pseudo-linearly towards the edges before tapering off into a peak. Instead, what I get is nearly the exact opposite, a non-zero flatline everywhere except the edges. enter image description here

So, now the question is, how do I solve this so that it produces appropriate potential for a wave that sits on the "edges" of the circle? I have been coding this in Julia, so I will include the solution code below.

#M is the number of points, g is the waveform, which we get from |psi|^2
#minx = 0, maxx = 2pi
#dx is the step size, C is the FDM matrix
using SparseArrays, LinearAlgebra, LinearSolve, Plots, Printf, LaTeXStrings

#Define the FDM matrix
ctop = zeros(Float64, M-1)
cmid = zeros(Float64, M)
cbot = zeros(Float64, M-1)
@inbounds for i in 1:M-1
    cmid[i] = -2/dx^2
    ctop[i] = 1 /dx^2
    cbot[i] = 1/dx^2
cmid[M] = -2/dx^2
C = spdiagm(-1=>cbot, 0=>cmid, 1=>ctop)
C[M,1] = 1/dx^2

#Set up g, the waveform, it should be noted that init_func produces a gaussian on the circle centered at pi of width sigma=1
psi= init_func(minx, maxx, M, [pi,0,1])
abspsi = abs.(psi).^2
g = zeros(Float64, M)

@inbounds for i in 1:M
    g[i] = abs(psi[i])^2
    #g[i] = -4 *pi * abspsi[i]

#Now solve then plot
phi = zeros(Float64, M)
phi = C\g
  • $\begingroup$ Firstly shouldn't you set C[ 1, M ]? Secondly for periodic boundaries the matrix is singular, are you aware of the implications? $\endgroup$
    – Ian Bush
    Commented Jul 13, 2023 at 19:16
  • $\begingroup$ Are you solving $\partial_t u(t,x) = \partial_{xx} u(t,x) + f(t,x)$ with periodic boundaries, or $\partial_{t}u(t,x) = \partial_{xx} u(t,x) + f(x)$, or $\partial_t u(t,x) = \partial_{xx} u(t,x)$? (I can't read Julia sorry.) What do you want to happen at time infinity, what do you want the steady state to be? $\endgroup$
    – lightxbulb
    Commented Jul 13, 2023 at 19:22
  • 1
    $\begingroup$ @lightxbulb From the Julia there's no time stepping here, the problem is Poisson's equation in 1D with PBC if I understand things correctly. $\endgroup$
    – Ian Bush
    Commented Jul 13, 2023 at 19:44
  • $\begingroup$ @IanBush Thank you for the clarification. TheAkashain actually said Poisson equation, idk why I thought of diffusion equations, my brain was fried probably. $\endgroup$
    – lightxbulb
    Commented Jul 13, 2023 at 21:06
  • $\begingroup$ What are $\phi,\psi$? What is the equation being solved? $\endgroup$ Commented Jul 14, 2023 at 4:24

1 Answer 1


Here's what I think the problem is. You have the equation: $$\Delta u(x) = f(x), \, x\in \Omega.$$ where $\Omega$ is a circle/torus. First you have some compatibility constraints you have to satisfy. Using the divergence theorem you get: $$\int_{\Omega} f(x) \,dx = \int_{\Omega} div(\nabla u(x)) \,dx = \int_{\partial\Omega} \partial_n u(x) \,dx = 0.$$ So if your potential does not integrate to zero you must make it integrate to zero. You can achieve this by replacing it with $g_i = f_i - f_{avg}$ where $f_{avg} = (1^Tf)/n$ is the average of $f$, then $1^T g = 0$.

Next, note that you have infinitely many solutions: you can add a constant to any of those and it will still be a solution. So we also need a constraint on the integral of the solution, e.g. $\int_{\Omega}u(x)\,dx = |\Omega|\mu$, if you want to have the mean value of $f$ then $\mu = (\int_{\Omega}f(x)\,dx)/|\Omega|$. In the discrete setting this reads $1^T u = n\mu = 1^Tf$. Let $M$ discretise $\Delta$, you can now use $M' = M - 11^T/n$ where $11^T$ is the matrix of ones everywhere, and your new system is: $$M' u = g - \mu 1.$$ I picked a minus because $M$ is negative semi-definite. That's why I generally prefer to write: $$-\Delta u(x) = f(x).$$ In that case with $M\approx -\Delta$, $M$ is positive semi-definite, and $M'= M+11^T/n$ is positive definite, and the system is $M' u = g +\mu 1$. However if $\mu = f_{avg}$ you can skip the $g$ and just solve $M' u = f$.

I should also mention here that the matrices $M$ and $M'$ are diagonalized on a torus with the discrete Fourier transform (since they are circulant), e.g. $M = V\Lambda V^*$, where $V^*$ is the discrete Fourier transform, and $V$ is the backtransform. So you can solve this as follows: $u = M^{+} f = V \Lambda^{+} V^* f + \mu 1$ or if $\mu = f_{avg}$ then just $u = (M')^{-1} = V (\Lambda')^{-1} V^* f$ (it's the same thing).

Also a note if you use a (reasonable) iterative solver - you don't really need to keep around a matrix with ones (in fact you don't need to keep any matrix around if you implement the mat-vec multiply yourself), you can just set the initial guess to have mean $\mu$ and then use the matrix $M$ in the iterations with right-hand side $g$, the residual for the constant eigenvector would then be zero, and you would converge to the correct thing. For example this works with the conjugate gradient method.

  • $\begingroup$ So, let me see if I can apply my problem to this! Let's say that the potential is a gravitational self potential defined by a normalized wave, $\psi$. In this case, the potential is $f=|\psi|^2$. So, a valid statement of the poisson equation should be $\partial_x \phi = |\psi|^2$. If we use a crank nicolson setup, then would what you wrote basically be "the triangular matrix from the comments, subtract a matrix of the same sized filled with 1/n in every entry, multipled by u, equals $g$ - the integral of f divided by the length of the circle"? Assuming I define g as you did, using $f-f_{avg}$. $\endgroup$ Commented Jul 13, 2023 at 22:16
  • $\begingroup$ It should read $d_{xx} \phi = |\psi|^2$ (or $-d_{xx}\phi = |\psi|^2$, I prefer this one since the operator is positive semi-definite). The first step is yo subtract the mean from the right hand side, i.e. $ f= |\psi|^2$, compute $g = f - f_{avg}1$ where $f_{avg} = 1^Tf/n$. Now if you use $d_{xx}$, i.e. $h^{-2}\begin{bmatrix} 1 & -2 & 1 \end{bmatrix}$ then add $-1/n$ to every element of your matrix resulting in $M'$. Let $\mu$ be the mean you want your solution to have, then solve $M'\phi =g - \mu 1$. If you use $-d_{xx}$, i.e. $\begin{bmatrix} -1 & 2 & -1\end{bmatrix}$ then add $1/n$ and... $\endgroup$
    – lightxbulb
    Commented Jul 13, 2023 at 22:34
  • $\begingroup$ ... solve $M'\phi = g + \mu 1$. If $\mu = 1^T f / n$ then you do not need to compute $g$ in this second case and you can just solve $M'\phi = f$ instead. I don't see how you can use Crank-Nicolson considering this is not an evolution equation. I personally would just compute the spectral solution using DFT/FFT. Another alternative is to use the conjugate gradient method since $M$ is symmetric positive semi-definite. Richardson/Jacobi would also do, though it is not the fastest (you need a step size of $\tau\in (0,1/2)$ for those to be stable). Gauss-Seidel and SOR are ok, but CG beats them. $\endgroup$
    – lightxbulb
    Commented Jul 13, 2023 at 22:36
  • $\begingroup$ thank you so much! Apologies for using the CN term, I forgot the term for the [1 -2 1]/h^2 method! I'll try this out tonight, and I will mark your answer as the solution! $\endgroup$ Commented Jul 13, 2023 at 22:44
  • $\begingroup$ There are some typos in my two comments above, I fix them here: $-d_{xx} \approx h^{-2}\begin{bmatrix} -1 & 2 & -1 \end{bmatrix}$, I had missed the $h^{-2}$ (this is not CN, just finite differences). The step size should read $\tau \in (0,h^2/2)$. Note that in 2-D this is $\tau \in (0,h^2/4)$, and I think in $m$-D you have $\tau \in(0,h^2/(2m))$. Then again the conjugate gradient solver (CG) doesn't require you to specify a step size, and you can implement the mat-vec mult matrix-free. I still recommend dft/fft. Full multigrid may be comparable in terms of speed to FFT. $\endgroup$
    – lightxbulb
    Commented Jul 13, 2023 at 22:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.