Уравнение Пуассона, решение дискретным преобразованием Фурье

Основные авторы описания: В.М.Степаненко, Е.В.Мортиков

1 Описание свойств и структуры алгоритма

1.1 Общее описание задачи

Уравнение Пуассона для многомерного пространства имеет вид:

[math] \begin{equation} \label{eq:poisson} \sum_{i=1}^{N}\frac{\partial^2 \phi}{\partial x_i^2}=f,~\mathbf{x}\in D, \end{equation} [/math]

где [math]D \in \mathbb{R}^N[/math] - область определения решения [math]\phi(\mathbf{x})[/math], [math]\mathbf{x}=(x_1,...,x_N)^T[/math] - вектор независимых переменных. Уравнение Пуассона дополняется граничными условиями:

[math] \begin{equation} B(\phi)=F,~ \text{на}~ \mathbf{x}\in \Gamma(D), \end{equation} [/math]

где [math]\Gamma(D)[/math] - граница области [math]D[/math], a [math]B(\phi)[/math] - оператор, определяющий граничные условия. [math]B(\phi)=\phi[/math] соответствует граничным условиям Дирихле, [math]B(\phi)=\partial\phi/\partial n[/math] ([math]\mathbf{n}[/math] - внешняя нормаль к границе [math]\Gamma(D)[/math]) - условиям Неймана. Иногда задают также смешанное граничное условие: [math]B(\phi)=С\phi+\partial\phi/\partial n[/math] ([math]C[/math] - константа). Встречаются также так называемые "периодические граничные условия", при которых задача ставится для бесконечной области, но предполагается периодичность решения по подмножеству переменных из [math]\mathbf{x}[/math].

Уравнение Пуассона возникает во многих задачах математической физики, например, в электростатике (в этом случае [math]\phi[/math] - потенциал электрической силы) и гидродинамике ([math]\phi[/math] - давление жидкости или газа); при этом [math]N=2,3[/math] для плоской и объемной задач, соответственно.

Решение уравнения Пуассона для произвольной правой части и неоднородных граничных условий в аналитическом виде не может быть получено, поэтому в большинстве приложений оно находится численно. Наиболее распространенная дискретизация уравнения имеет вид

[math] \begin{equation} \sum_{i=1}^{N}\frac{\phi_{k_1,...,k_i+1,...,k_N}-2\phi_{k_1,...,k_i,...,k_N}+\phi_{k_1,...,k_i-1,...,k_N}}{\Delta x_i^2}=f_{k_1,...,k_N},~(k_1,...,k_N) \in D_N. \end{equation} [/math]

Здесь вторые производные заменены центральными разностями второго порядка точности, а решение ищется на дискретном множестве точек [math]N[/math]-мерного пространства, [math]D_N[/math]. Конечными разностями аппроксимируются при этом также граничные условия.

1.2 Математическое описание

В настоящей статье мы рассмотрим конечно-разностную схему для наиболее часто встречающейся задачи решения уравнения Пуассона в трехмерном пространстве:

[math] \begin{equation} \label{eq:poisdiscr} \frac{\phi_{i+1,j,k}-2\phi_{i,j,k}+\phi_{i-1,j,k}}{ \Delta x^2} + \frac{\phi_{i,j+1,k}-2\phi_{i,j,k}+\phi_{i,j-1,k}}{ \Delta y^2} + \frac{\phi_{i,j,k+1}-2\phi_{i,j,k}+\phi_{i,j,k-1}}{ \Delta z^2} = f_{i,j,k},~(i,j,k) \in D_h, \end{equation} [/math]

где [math]D_h=\{0:N_x-1\}\times\{0:N_y-1\}\times\{0:N_z-1\} [/math] - параллелепипед в сеточной области. В качестве граничных условий для простоты примем так называемые троякопериодические граничные условия

[math] \begin{align} \phi_{0,j,k}=\phi_{N_x,j,k},\\ \phi_{i,0,k}=\phi_{i,N_y,k},\\ \phi_{i,j,0}=\phi_{i,j,N_z}. \end{align} [/math]

Периодические граничные условия удовлетворяются "автоматически", если представить решение в виде стандартного дискретного обратного преобразования Фурье:

[math] \begin{equation} \label{eq:fourier} \phi_{i,j,k}=\frac{1}{N_x N_y N_z}\sum_{l=0}^{N_x-1}\sum_{m=0}^{N_y-1}\sum_{n=0}^{N_z-1}\Phi_{l,m,n}e^{2\pi \overline{i} \left(\frac{il}{N_x}+\frac{jm}{N_y}+\frac{kn}{N_z}\right)}. \end{equation} [/math]

Здесь [math]\overline{i}=\sqrt{-1}[/math]. Аналогичное разложение применяется также к правой части уравнения, [math]f_{i,j,k}[/math]. Решение дисретного уравнения Пуассона преобразованием Фурье удобно тем, что базисные функции разложения Фурье являются собственными функциями дискретного оператора Лапласа. А именно, подставляя разложение Фурье неизвестной функции [math]\phi_{i,j,k}[/math] и правой части в исходное уравнение, получаем:

[math] \begin{equation} \Phi_{l,m,n}=-\frac{F_{l,m,n}}{4\left[\sin^2\left(\frac{\pi l}{N_x}\right) + \sin^2\left(\frac{\pi m}{N_y}\right) + \sin^2\left(\frac{\pi n}{N_z}\right) \right]}. \end{equation} [/math]

где [math]F_{l,m,n}[/math] - преобразование Фурье правой части. Отсюда очевиден алгоритм решения уравнения: сначала правая часть уравнения разлагается в ряд Фурье, затем по приведенной выше формуле вычисляются коэффициенты Фурье решения и, наконец, решение восстанавливается обратным преобразованием Фурье.

1.3 Вычислительное ядро алгоритма

Вычислительным ядром алгоритма является одномерное преобразование Фурье. В самом деле, обратное дискретное преобразование Фурье может быть записано как

[math] \begin{equation} \label{eq:fourier2} \phi_{i,j,k}=\frac{1}{N_x} \sum_{l=0}^{N_x-1} \left[ e^{2\pi \overline{i}\left(\frac{il}{N_x}\right)} \frac{1}{N_y} \sum_{m=0}^{N_y-1} \left[ e^{2\pi \overline{i}\left(\frac{jm}{N_y}\right)} \frac{1}{N_z} \sum_{n=0}^{N_z-1} \Phi_{l,m,n}e^{2\pi \overline{i} \left(\frac{kn}{N_z}\right)}\right]\right], \end{equation} [/math]

откуда видно, что трехмерное преобразование Фурье сводится в последовательности трех одномерных. Для вычисления одномерного преобразования широко используется эффективный алгоритм быстрого преобразования Фурье (FFT, Fast Fourier Transform).

1.4 Макроструктура алгоритма

Из приведенного выше ясно, что основной макрооперацией алгоритма решения уравнения Пуассона является одномерное быстрое преобразование Фурье. В дальнейшем будем обозначать ее [math]\text{FFT}_i,~i=x,y,z[/math] по каждому из трех направлений, соответственно, а обратное преобразование Фурье - [math]\text{IFFT}_i,~i=x,y,z[/math].

1.5 Описание схемы реализации последовательного алгоритма

Введем для краткости следующие обозначения сеточных функций [math]\phi_h=\{\phi_{i,j,k},~(i,j,k) \in D_h\},~F_h=\{F_{l,m,n},~l=0,...,N_x-1,~m=0,...,N_y-1,~n=0,...,N_z-1\}[/math] и аналогично [math]f_h[/math] и [math]\Phi_h[/math]. Тогда алгоритм запишется следующим образом:

1. Вычисление [math]\text{FFT}_x(f_h)[/math]
2. Вычисление [math]\text{FFT}_y(\text{FFT}_x(f_h))[/math]
3. Вычисление [math]\text{FFT}_z(\text{FFT}_y(\text{FFT}_x(f_h)))=F_h[/math]
4. Расчет [math]\Phi_h[/math] по [math]F_h[/math]
5. Вычисление [math]\text{IFFT}_x(\Phi_h)[/math]
6. Вычисление [math]\text{IFFT}_y(\text{IFFT}_x(\Phi_h))[/math]
7. Вычисление [math]\text{IFFT}_z(\text{IFFT}_y(\text{IFFT}_x(\Phi_h)))=\phi_h[/math]

Примечательно, что весь описанный алгоритм реализуется с помощью одного трехмерного массива [math]N_x \times N_y \times N_z[/math], поскольку результат каждого одномерного преобразования Фурье можно записать во входной массив и затем использовать его на вход для следующего преобразования, а операция (4) - это поэлементное изменение массива.

1.6 Последовательная сложность алгоритма

Рассмотрим случай куба, [math]N_x=N_y=N_z=N[/math]. Тогда сложность одномерного БПФ по каждому из направлений составляет [math]N(\text{log}_2 N)[/math] операций. Поскольку одномерное БПФ ка каждом этапе алгоритма осуществляется [math]N^2[/math] раз, таких этапов 6, а количество операций на этапе 4 составляет [math]N^3[/math], то общее количество действий составляет [math]N^3(2\text{log}_2 N+1)[/math].

1.7 Информационный граф

Информационный граф алгоритма выглядит слеующим образом

Красным кругом обозначена операция одномерного БПФ, например [math]FFT_{x11}[/math] означает БПФ вдоль оси x при [math]j=1,~k=1[/math]. Красными кругами представлены поэлементные преобразования массива для получения коэффициентов Фурье решения (этап 4 алгоритма), например [math]L_{2N,1,...,N}[/math] означает преобразование элементов [math]i=2,~j=N,~k=1,...,N[/math]. Целочисленные индексы для краткости записи перенумерованы от [math]0,...,N-1[/math] до [math]1,...,N[/math]

Видно, что