Taxa de convergência do FFT Poisson solver

Qual é a taxa de convergência teórica para um solucionador de FFT Poison?

Estou resolvendo uma equação de Poisson: com n ( x , y , z ) =

$$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$$ no domínio[0,2]×[0,2]×[0,2]com condição de contorno periódica. Essa densidade de carga é neutra em termos líquidos. A solução é dada por: VH(x)=∫n $$n(x, y, z) = {3\over\pi} ((x-1)^2 + (y-1)^2 + (z-1)^2 - 1)$$ $[0, 2] \times [0, 2] \times [0, 2]$ ondex=(x,y,z). No espaço recíproco VH(G)=4πn(G $$V_H({\bf x}) = \int {n({\bf y})\over |{\bf x}-{\bf y}|} d^3 y$$ ${\bf x}=(x, y, z)$ ondeGsão os vetores espaciais recíprocos. Estou interessado na energia Hartree: EH= $$V_H({\bf G}) = 4\pi {n({\bf G})\over G^2}$$ ${\bf G}$ No espaço recíproco, isso se torna (após discretização): EH=2π∑G≠0 | n( G ) | $$E_H = {1\over 2} \int {n({\bf x}) n({\bf y})\over |{\bf x}-{\bf y}|} d^3 x d^3 y = {1\over 2} \int V_H({\bf x}) n({\bf x}) d^3 x$$ OtermoG=0é omitido, o que efetivamente neutraliza a rede de densidade de carga (e, como já é neutra, tudo é consistente). $$E_H = 2\pi \sum_{{\bf G}\ne 0} {|n({\bf G})|^2\over G^2}$$ ${\bf G}=0$

$$E_H = {128\over 35\pi} = 1.16410...$$

Aqui está um programa usando o NumPy que faz o cálculo.

from numpy import empty, pi, meshgrid, linspace, sum
from numpy.fft import fftn, fftfreq
E_exact = 128/(35*pi)
print "Hartree Energy (exact):      %.15f" % E_exact
f = open("conv.txt", "w")
for N in range(3, 384, 10):
    print "N =", N
    L = 2.
    x1d = linspace(0, L, N)
    x, y, z = meshgrid(x1d, x1d, x1d)

    nr = 3 * ((x-1)**2 + (y-1)**2 + (z-1)**2 - 1) / pi
    ng = fftn(nr) / N**3

    G1d = N * fftfreq(N) * 2*pi/L
    kx, ky, kz = meshgrid(G1d, G1d, G1d)
    G2 = kx**2+ky**2+kz**2
    G2[0, 0, 0] = 1  # omit the G=0 term

    tmp = 2*pi*abs(ng)**2 / G2
    tmp[0, 0, 0] = 0  # omit the G=0 term
    E = sum(tmp) * L**3
    print "Hartree Energy (calculated): %.15f" % E
    f.write("%d %.15f\n" % (N, E))
f.close()

E aqui está um gráfico de convergência (apenas plotando o conv.txtscript acima, aqui está um caderno que faz isso se você quiser brincar com isso):

Como você pode ver, a convergência é linear, o que foi uma surpresa para mim, pensei que a FFT converge muito mais rápido que isso.

Atualização :

A solução tem uma cúspide na fronteira (eu não percebi isso antes). Para que a FFT converja rapidamente, a solução deve ter todos os derivados suaves. Então, eu também tentei o seguinte lado direito:

nr = 3*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)/4

$V_H=\sin(\pi x)\sin(\pi y)\sin(\pi z)$$E_H = \frac{3\pi}{8}$

Alguém conhece alguma referência em 3D para que eu possa ver uma convergência mais rápida do que linear?

Deixe-me primeiro responder a todas as perguntas:

Qual é a taxa de convergência teórica para um solucionador de FFT Poison?

A convergência teórica é exponencial desde que a solução seja suficientemente suave.

Quão rápido essa energia deve convergir?

$E_H$

Alguém conhece alguma referência em 3D para que eu possa ver uma convergência mais rápida do que linear?

Qualquer lado direito que produza uma solução periódica e infinitamente diferenciável (inclusive além da fronteira periódica) deve convergir exponencialmente.

---

No código acima, ocorre um erro, que faz com que a convergência seja mais lenta que a exponencial. Usando o código de densidade suave ( https://gist.github.com/certik/5549650/ ), o seguinte patch corrige o erro:

@@ -6,7 +6,7 @@ f = open("conv.txt", "w")
 for N in range(3, 180, 10):
     print "N =", N
     L = 2.
-    x1d = linspace(0, L, N)
+    x1d = linspace(0, L, N+1)[:-1]
     x, y, z = meshgrid(x1d, x1d, x1d)

     nr = 3*pi*sin(pi*x)*sin(pi*y)*sin(pi*z)/4

O problema era que a amostragem no espaço real não pode repetir o primeiro e o último ponto (que tem o mesmo valor devido à condição de contorno periódica). Em outras palavras, o problema estava em configurar a amostragem inicial.

Após essa correção, a densidade converge em uma iteração, como Matt disse acima. Então, eu nem plotei os gráficos de convergência.

No entanto, pode-se tentar uma densidade mais difícil, por exemplo:

     nr = 3*pi*exp(sin(pi*x)*sin(pi*y)*sin(pi*z))/4

então a convergência é exponencial, como esperado. Aqui estão os gráficos de convergência para essa densidade: