Скорость сходимости FFT-пуассоновского решателя

16

Какова теоретическая скорость сходимости для решения FFT Poison?

Я решаю уравнение Пуассона: с n ( x , y , z ) = 3

2ВЧАС(Икс,Y,Z)знак равно-4πN(Икс,Y,Z)
в области[0,2]×[0,2]×[0,2]с периодическим граничным условием. Эта плотность заряда является нейтральной. Решение имеет вид: VH(x)=n(
N(Икс,Y,Z)знак равно3π((Икс-1)2+(Y-1)2+(Z-1)2-1)
[0,2]×[0,2]×[0,2] гдеx=(x,y,z). В обратном пространстве VH(G)=4πn(G)
ВЧАС(Икс)знак равноN(Y)|Икс-Y|d3Y
Иксзнак равно(Икс,Y,Z) гдеG- обратные пространственные векторы. Меня интересует энергия Хартри: EH=1
ВЧАС(грамм)знак равно4πN(грамм)грамм2
грамм В обратном пространстве это становится (после дискретизации): EH=2πG0 | n( G ) | 2
ЕЧАСзнак равно12N(Икс)N(Y)|Икс-Y|d3Иксd3Yзнак равно12ВЧАС(Икс)N(Икс)d3Икс
ЧленG=0опущен, что фактически делает плотность заряда чистой нейтральной (а поскольку она уже нейтральная, то все согласованно).
ЕЧАСзнак равно2πΣграмм0|N(грамм)|2грамм2
граммзнак равно0

ЕЧАСзнак равно12835πзнак равно1,16410 ...

Вот программа, использующая NumPy, которая выполняет вычисления.

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

А вот график сходимости (просто построение приведенного conv.txtвыше сценария, вот блокнот, который делает это, если вы хотите поиграть с этим самостоятельно):

БПФ граф сходимости

Как видите, конвергенция линейная, что было для меня неожиданностью, я думал, что FFT сходится гораздо быстрее, чем это.

Обновление :

У решения есть острый край на границе (я не осознавал этого раньше). Чтобы БПФ быстро сходилось, решение должно иметь все гладкие производные. Поэтому я также попробовал следующую правую часть:

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

ВЧАСзнак равногрех(πИкс)грех(πY)грех(πZ)ЕЧАСзнак равно3π8

Кто-нибудь знает какой-нибудь эталон в 3D, чтобы я мог видеть более быструю конвергенцию, чем линейную?

Ондржей Чертик
источник
Ондрей, не является ли преобразование Фурье вашей гладкой плотности дельта-функцией? Я признаюсь, что слишком ленив, чтобы запустить его, но он должен получить точный ответ с первой попытки.
Мэтт Кнепли
Я думаю, что это. Но это не сходится за одну итерацию, как видно из графиков тетради. Я понятия не имею, что происходит.
Ондржей Чертик
Ондрей, ты уверен, что твоя реализация верна? Я помню, как пытался использовать спектральные решатели для домашнего задания в аспирантуре и полностью отбрасывал константы. Я заметил, что вы измеряете погрешность, глядя на абсолютное расстояние между вычисленной и точной энергией. Как выглядит ваше сближение с фактическим решением проблемы? Это должно быть легко вычислить и даже построить на двухмерной части задачи.
Арон Ахмадиа
Арон --- Я проверил свою реализацию на предмет какого-то другого кода, но я проверил его на предмет неправильной начальной выборки, поэтому у меня была одинаковая ошибка в обоих кодах. Мэтт был прав, теперь он сходится с первой попытки. Смотрите мой ответ ниже.
Ондржей Чертик

Ответы:

10

Позвольте мне сначала ответить на все вопросы:

Какова теоретическая скорость сходимости для решения FFT Poison?

Теоретическая сходимость является экспоненциальной, пока решение является достаточно гладким.

Как быстро эта энергия должна сходиться?

ЕЧАС

Кто-нибудь знает какой-нибудь эталон в 3D, чтобы я мог видеть более быструю конвергенцию, чем линейную?

Любая правая часть, которая дает периодическое и бесконечно дифференцируемое решение (в том числе через периодическую границу), должна сходиться экспоненциально.


В приведенном выше коде есть ошибка, которая вызывает сходимость медленнее, чем экспоненциальная. Используя код с гладкой плотностью ( https://gist.github.com/certik/5549650/ ), следующий патч исправляет ошибку:

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

Проблема заключалась в том, что реальная космическая выборка не может повторить первую и последнюю точку (которая имеет одинаковое значение из-за периодических граничных условий). Другими словами, проблема заключалась в настройке начальной выборки.

После этого исправления плотность сходится за одну итерацию, как сказал Мэтт выше. Поэтому я даже не строил графики сходимости.

Однако можно попробовать более сложную плотность, например:

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

тогда сходимость является экспоненциальной, как и ожидалось. Вот графики сходимости для этой плотности: введите описание изображения здесь введите описание изображения здесь

Ондржей Чертик
источник
Потрясающие. Извините, я больше не помогал!
Арон Ахмадиа