Tasa de convergencia del solucionador FFT de Poisson

16

¿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 ) = 3

2VH(X,y,z)=-4 4πnorte(X,y,z)
[0,2]×[0,2]×[0,2]VH(x)=n( y )
norte(X,y,z)=3π((X-1)2+(y-1)2+(z-1)2-1)
[0 0,2]×[0 0,2]×[0 0,2]x=(x,y,z)VH(G)=4πn(G)
VH(X)=norte(y)El |X-yEl |re3y
X=(X,y,z) GEH=1
VH(sol)=4 4πnorte(sol)sol2
sol
miH=12norte(X)norte(y)El |X-yEl |re3Xre3y=12VH(X)norte(X)re3X
En el espacio recíproco esto se convierte (después de la discretización):
miH=2πsol0 0El |norte(sol)El |2sol2
El término sol=0 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:

miH=12835π=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):

Gráfico de convergencia FFT

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).VH=pecado(πX)pecado(πy)pecado(πz)miH=3π8

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

Ondřej Čertík
fuente
Ondrej, ¿no es la transformación de Fourier de su densidad suave una función delta? Admito que soy demasiado vago para ejecutarlo, pero debería obtener la respuesta exacta en el primer intento.
Matt Knepley
Creo que es. Pero no converge en una iteración, como se puede ver en los diagramas del cuaderno. No tengo idea de lo que está pasando.
Ondřej Čertík
Ondrej, ¿estás seguro de que tu implementación es correcta? Recuerdo haber intentado implementar solucionadores espectrales para una tarea en la escuela de posgrado y sacudir totalmente las constantes. Noto que está midiendo el error al observar la distancia absoluta entre la energía calculada y la exacta. ¿Cómo se ve su convergencia con la solución real del problema? Esto debería ser fácil de calcular e incluso trazar en una porción 2D del problema.
Aron Ahmadia
Aron --- Verifiqué mi implementación con algún otro código, pero estaba comprobando que mi muestreo inicial era incorrecto, por lo que tuve el mismo error en ambos códigos. Matt tenía razón, ahora converge en el primer intento. Vea mi respuesta a continuación.
Ondřej Čertík

Respuestas:

10

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.miH

¿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: ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

Ondřej Čertík
fuente
Increíble. Lo siento, no fui más ayuda!
Aron Ahmadia