Konvergenzrate des FFT-Poisson-Lösers

Was ist die theoretische Konvergenzrate für einen FFT-Giftlöser?

Ich löse eine Poisson-Gleichung:

$$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$$ mit in der Domäne[0,2]×[0,2]×[0,2]mit periodischer Randbedingung. Diese Ladungsdichte ist netto neutral. Die Lösung ist gegeben durch: 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]$ mitx=(x,y,z). Im reziproken Raum ist 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)$ wobeiGdie reziproken Raumvektoren sind. Ich interessiere mich für die Hartree-Energie: EH= $$V_H({\bf G}) = 4\pi {n({\bf G})\over G^2}$$ ${\bf G}$ Im reziproken Raum wird dies (nach Diskretisierung): 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$$ DerTermG=0wird weggelassen, wodurch die Ladungsdichte effektiv netto neutral wird (und da sie bereits neutral ist, ist alles konsistent). $$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...$$

Hier ist ein Programm mit NumPy, das die Berechnung durchführt.

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()

Und hier ist ein Konvergenzdiagramm (nur conv.txtdas obige Skript, hier ist ein Notizbuch, das es macht, wenn Sie selbst damit spielen möchten):

Wie Sie sehen, ist die Konvergenz linear, was für mich eine Überraschung war. Ich dachte, dass FFT viel schneller konvergiert.

Update :

Die Lösung hat eine Spitze an der Grenze (das habe ich vorher nicht bemerkt). Damit FFT schnell konvergiert, müssen alle Ableitungen der Lösung glatt sein. Also habe ich auch die folgende rechte Seite ausprobiert:

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}$

Kennt jemand einen Benchmark in 3D, so dass ich eine schnellere Konvergenz als linear sehen kann?

Lassen Sie mich zunächst alle Fragen beantworten:

Was ist die theoretische Konvergenzrate für einen FFT-Giftlöser?

Die theoretische Konvergenz ist exponentiell, solange die Lösung ausreichend glatt ist.

Wie schnell soll diese Energie konvergieren?

$E_H$

Kennt jemand einen Benchmark in 3D, so dass ich eine schnellere Konvergenz als linear sehen kann?

Jede rechte Seite, die eine Lösung erzeugt, die periodisch und unendlich differenzierbar ist (auch über die periodische Grenze hinweg), sollte exponentiell konvergieren.

---

Im obigen Code ist zufällig ein Fehler aufgetreten, der dazu führt, dass die Konvergenz langsamer als exponentiell ist. Mit dem Code für die glatte Dichte ( https://gist.github.com/certik/5549650/ ) behebt der folgende Patch den Fehler:

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

Das Problem bestand darin, dass die Realraumabtastung den ersten und den letzten Punkt nicht wiederholen kann (die aufgrund der periodischen Randbedingung denselben Wert haben). Das Problem bestand also darin, die Erstbemusterung einzurichten.

Nach dieser Korrektur konvergiert die Dichte in einer Iteration, wie Matt oben sagte. Also habe ich nicht einmal die Konvergenzgraphen gezeichnet.

Man kann jedoch eine schwierigere Dichte ausprobieren, zum Beispiel:

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

dann ist die Konvergenz erwartungsgemäß exponentiell. Hier sind die Konvergenzdiagramme für diese Dichte: