Tingkat konvergensi pemecah FFT Poisson

Berapa tingkat konvergensi teoretis untuk pemecah FFT Poison?

Saya memecahkan persamaan Poisson: dengan n ( x , y , z ) =

$$\nabla^2 V_H(x, y, z) = -4\pi n(x, y, z)$$ pada domain[0,2]×[0,2]×[0,2]dengan kondisi batas periodik. Kerapatan muatan ini netral bersih. Solusi diberikan oleh: 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]$ manax=(x,y,z). Dalam ruang timbal balik 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)$ di manaGadalah vektor ruang timbal balik. Saya tertarik pada energi Hartree: EH= $$V_H({\bf G}) = 4\pi {n({\bf G})\over G^2}$$ ${\bf G}$ Dalam ruang timbal balik, ini menjadi (setelah diskritisasi): 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$$ IstilahG=0dihilangkan, yang secara efektif membuat kerapatan muatan menjadi netral (dan karena sudah netral, maka semuanya konsisten). $$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...$$

Berikut adalah program menggunakan NumPy yang melakukan perhitungan.

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

Dan di sini adalah grafik konvergensi (hanya memplot conv.txtdari script di atas, di sini adalah notebook yang melakukan itu jika Anda ingin bermain dengan ini sendiri):

Seperti yang Anda lihat, konvergensi linier, yang merupakan kejutan bagi saya, saya pikir FFT bertemu lebih cepat dari itu.

Perbarui :

Solusinya memiliki titik puncak pada batas (saya tidak menyadari ini sebelumnya). Agar FFT dapat konvergen dengan cepat, solusinya harus memiliki semua turunannya lancar. Jadi saya juga mencoba sisi kanan berikut:

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

Adakah yang tahu benchmark apa pun dalam 3D sehingga saya bisa melihat konvergensi yang lebih cepat daripada linear?

Biarkan saya menjawab semua pertanyaan:

Berapa tingkat konvergensi teoretis untuk pemecah FFT Poison?

Konvergensi teoretis adalah eksponensial selama solusinya cukup lancar.

Seberapa cepat energi ini bertemu?

Energi Hartree $E_H$harus menyatu secara eksponensial untuk solusi yang cukup halus. Jika solusinya kurang mulus, maka konvergensi lebih lambat.

Adakah yang tahu benchmark apa pun dalam 3D sehingga saya bisa melihat konvergensi yang lebih cepat daripada linear?

Sisi kanan mana pun yang menghasilkan solusi yang periodik dan dapat dibedakan tanpa batas (termasuk melintasi batas periodik) harus bertemu secara eksponensial.

---

Dalam kode di atas ada bug, yang menyebabkan konvergensi lebih lambat daripada eksponensial. Menggunakan kode kepadatan halus ( https://gist.github.com/certik/5549650/ ), tambalan berikut memperbaiki bug:

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

Masalahnya adalah bahwa pengambilan sampel ruang nyata tidak dapat mengulangi titik pertama dan terakhir (yang memiliki nilai yang sama karena kondisi batas periodik). Dengan kata lain, masalahnya adalah menyiapkan sampel awal.

Setelah perbaikan ini, kepadatan berkumpul dalam satu iterasi, seperti yang Matt katakan di atas. Jadi saya bahkan tidak memetakan grafik konvergensi.

Namun, seseorang dapat mencoba kepadatan yang lebih sulit, misalnya:

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

maka konvergensi bersifat eksponensial, seperti yang diharapkan. Berikut adalah grafik konvergensi untuk kerapatan ini: