Tasa de convergencia del solucionador FFT de Poisson

¿Cuál es la tasa de convergencia teórica para un solucionador de veneno FFT?

Estoy resolviendo una ecuación de Poisson: con en el dominio con periódico condición de frontera. Esta densidad de carga es neta neutra. La solución viene dada por: donde . En el espacio recíproco donde son los vectores espaciales recíprocos. Estoy interesado en la energía Hartree: n ( x , y , z ) =

$$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$$ [0,2]×[0,2]×[0,2]VH(x)=∫n( y $$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]$x=(x,y,z)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)$ GEH= $$V_H({\bf G}) = 4\pi {n({\bf G})\over G^2}$$ ${\bf 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$$ En el espacio recíproco esto se convierte (después de la discretización): $$E_H = 2\pi \sum_{{\bf G}\ne 0} {|n({\bf G})|^2\over G^2}$$ El término ${\bf G}=0$ se omite, lo que hace que la densidad de carga sea neutra neta (y dado que es ya neutral, entonces todo es consistente).

Para el problema de prueba anterior, esto puede evaluarse analíticamente y uno obtiene:

$$E_H = {128\over 35\pi} = 1.16410...$$ ¿Qué tan rápido debería converger esta energía?

Aquí hay un programa que usa NumPy que hace el 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()

Y aquí hay un gráfico de convergencia (solo trazando el conv.txtscript anterior, aquí hay un cuaderno que lo hace si quieres jugar con esto tú mismo):

Como puede ver, la convergencia es lineal, lo que fue una sorpresa para mí, pensé que FFT converge mucho más rápido que eso.

Actualización :

La solución tiene una cúspide en el límite (no me di cuenta de esto antes). Para que FFT converja rápidamente, la solución debe tener todas las derivadas suaves. Así que también probé el siguiente lado derecho:

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

Que puedes poner en el script anterior ( script actualizado ). La solución exacta es , que debería ser infinitamente diferenciable. La integral exacta en este caso es . Sin embargo, el solucionador de FFT todavía converge linealmente hacia esta solución exacta, como se puede verificar ejecutando el script anterior y trazando la convergencia ( cuaderno actualizado con gráficos).$V_H=\sin(\pi x)\sin(\pi y)\sin(\pi z)$$E_H = \frac{3\pi}{8}$

¿Alguien conoce algún punto de referencia en 3D para que pueda ver una convergencia más rápida que lineal?

Permítanme primero responder todas las preguntas:

¿Cuál es la tasa de convergencia teórica para un solucionador de veneno FFT?

La convergencia teórica es exponencial siempre que la solución sea lo suficientemente suave.

¿Qué tan rápido debería converger esta energía?

La energía de Hartree debería converger exponencialmente para una solución suficientemente suave. Si la solución es menos fluida, entonces la convergencia es más lenta.$E_H$

¿Alguien conoce algún punto de referencia en 3D para que pueda ver una convergencia más rápida que lineal?

Cualquier lado derecho que produzca una solución que sea periódica e infinitamente diferenciable (incluso a través del límite periódico) debe converger exponencialmente.

---

En el código anterior, hay un error que hace que la convergencia sea más lenta que exponencial. Usando el código de densidad suave ( https://gist.github.com/certik/5549650/ ), el siguiente parche corrige el error:

@@ -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

El problema era que el muestreo en el espacio real no puede repetir el primer y el último punto (que tiene el mismo valor debido a la condición de límite periódica). En otras palabras, el problema estaba en configurar el muestreo inicial.

Después de esta solución, la densidad converge en una iteración, como Matt dijo anteriormente. Así que ni siquiera tracé los gráficos de convergencia.

Sin embargo, uno puede probar una densidad más difícil, por ejemplo:

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

entonces la convergencia es exponencial, como se esperaba. Aquí están los gráficos de convergencia para esta densidad: