Участник:Jhuighuy/Уравнение Пуассона, решение методом Галеркина

Уравнение Пуассона, решение методом Галеркина, Бутаков О.Б. 404 группа.

Свойства и структура алгоритмов

Общее описание алгоритма

Метод Галеркина (Галеркина-Бубнова) метод решения краевых задач для дифференциальных уравнений. В данной статье рассмотрим алгоритм метода Галеркина для решения уравнения Пуассона с краевыми условиями первого рода:

[math]\displaystyle{ \newcommand{\laplace}{\mathop{}\!\mathbin\bigtriangleup} \begin{cases} \laplace u = -f, u = \phi(x,y,z), (x,y,z) \in G, G = (0,L_x)\times(0,L_y)\times(0,L_z), \\ u(x,y,z) = 0, (x,y,z) \in \partial G. \end{cases} }[/math]

Алгоритм предполагает выбор полного набора базисных функций, удовлетворяющих краевым условиям. Решение строится в виде линейной комбинации базисных функций, коэффициенты находятся исходя из ортогональности невязки базисной фукнции.

Математическое описание алгоритма

Выберем тригонометрическую систему базисных функций:

[math]\displaystyle{ \phi_{m,n,l} = sin(\frac{\pi m x}{L_x})sin(\frac{\pi n y}{L_y})sin(\frac{\pi l z}{L_z}), m=1..N_x, n=1..N_y, l=1..N_z. }[/math]

Очевидно, эта система удовлетворяет краевым условиям. Проверим ортогональность функций, введя скалярное произведение [math]\displaystyle{ (u,v) = \int_{0}^{L_x}\int_{0}^{L_y}\int_{0}^{L_z} u v dx dy dz }[/math]:

[math]\displaystyle{ (\phi_{m,n,l} \phi_{i,j,k}) = \int_{0}^{L_x}\int_{0}^{L_y}\int_{0}^{L_z} \phi_{m,n,l} \phi_{i,j,k} dx dy dz = \int_{0}^{L_x} sin(\frac{\pi m x}{L_x}) sin(\frac{\pi i x}{L_x}) dx \int_{0}^{L_y} sin(\frac{\pi n x}{L_y}) sin(\frac{\pi j x}{L_y}) dy \int_{0}^{L_z} sin(\frac{\pi l x}{L_z}) sin(\frac{\pi k x}{L_z}) dz = \delta_{m,i}\frac{L_x}{2} \delta_{n,j}\frac{L_y}{2} \delta_{l,k}\frac{L_z}{2}. }[/math]

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

[math]\displaystyle{ \hat{u} = \Sigma_{m,n,l=1}^{L_x,L_y,L_z} C_{m,n,l} \phi_{m,n,l} }[/math].

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

[math]\displaystyle{ R(x,y,z) = \laplace\hat{u} + f = -\Sigma_{m,n,l=1}^{N_x,N_y,N_z} C_{m,n,l} \pi^2(\frac{m^2}{L_x^2}+\frac{n^2}{L_y^2}+\frac{l^2}{L_z^2}) \phi_{m,n,l} + f }[/math].

Разложив выражение для невязки, потребуем ортогональность невязки базисным функциям:

[math]\displaystyle{ \begin{split} &R_{i,j,k} = (R, \phi_{i,j,k}) = -C_{i,j,k} \pi^2(\frac{m^2}{L_x^2}+\frac{n^2}{L_y^2}+\frac{l^2}{L_z^2}) \frac{L_x}{2}\frac{L_y}{2}\frac{L_z}{2} + (f, \phi_{i,j,k}) = 0 \Rightarrow \\ &C_{i,j,k} = \frac{8 L_x L_y L_z (f, \phi_{i,j,k})}{\pi^2 ((i L_y L_z)^2 + (j L_x L_z)^2 + (k L_x L_y)^2)}, i=1..N_x, j=1..N_y. k=1..N_z. \end{split} }[/math].

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

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

[math]\displaystyle{ \omega_h = [(x_i,y_j,z_k); x_i=i\Delta x, y_j=i\Delta y, z_k=k\Delta k; L_x=M_x\Delta x, L_y=M_y\Delta y, L_z=M_z\Delta z]. }[/math]

Тогда трехмерный интеграл будет приближенно описываться выражением:

[math]\displaystyle{ \int_{0}^{L_x}\int_{0}^{L_y}\int_{0}^{L_z} f dx dy dz \approx \Sigma_{m,n,l=1}^{M_x,M_y,M_z} f(x_i,y_j,z_k)\Delta x\Delta y\Delta z. }[/math]

Вычисление функции на сетке тривиально:

[math]\displaystyle{ \hat{u}_{h,i,j,k} = \Sigma_{m,n,l=1}^{N_x,N_y,N_z} C_{m,n,l} sin(\frac{\pi m x_i}{L_x})sin(\frac{\pi n y_j}{L_y})sin(\frac{\pi l z_k}{L_z}). }[/math]

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

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

Первый вариант:

1. Для всех [math]\displaystyle{ i=1..N_x, j=1..N_y, k=1..N_z }[/math] вычислить [math]\displaystyle{ C_{m,n,l} }[/math].
2. Для всех [math]\displaystyle{ i=1..M_x, j=1..M_y, k=1..M_z }[/math] вычислить [math]\displaystyle{ \hat{u}_{h,i,j,k} }[/math].

Второй вариант:

1. Определим разбиение, например, [math]\displaystyle{ N_z }[/math] на [math]\displaystyle{ S }[/math] не пересекающихся подмножеств [math]\displaystyle{ [k_{0,s}, k_{1,s}]: \sqcup_{s=1}^{S}[k_{0,s}, k_{1,s}]=[1,N_z] }[/math]. Для всех [math]\displaystyle{ s=1..S: }[/math]
   1. Для всех [math]\displaystyle{ i=1..N_x, j=1..N_y, k=k_{0,s}..k_{1,s} }[/math] вычислить [math]\displaystyle{ C_{m,n,l} }[/math].
   2. Для всех [math]\displaystyle{ i=1..M_x, j=1..M_y, k=1..M_z }[/math] вычислить [math]\displaystyle{ \hat{u}_{h,i,j,k}^{s} = \Sigma_{m,n=1,l=k_{0,s}}^{N_x,N_y,k_{1,s}} C_{m,n,l} sin(\frac{\pi m x_i}{L_x})sin(\frac{\pi n y_j}{L_y})sin(\frac{\pi l z_k}{L_z}) }[/math].
2. Для всех [math]\displaystyle{ i=1..M_x, j=1..M_y, k=1..M_z }[/math] вычислить [math]\displaystyle{ \hat{u}_{h,i,j,k} = \Sigma_{s=1}^{S}\hat{u}_{h,i,j,k}^{s} }[/math]

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

Опишем лишь первый вариант макроструктуры, так как второй имеет смысл рассматривать лишь в случае последовательной реализации.

Первый этап на языке Fortran (для наглядности упростим вложенные конструкции):

   Do i j k = 1 1 1, Nx Ny Nz
       r = 0.0D0
       Do i0 j0 k0 = 1 1 1, Mx My Mz
           x = i0*dx
           y = j0*dy
           z = k0*dz
           r = r + dx*dy*dz*f(i0,j0,k0)*Sin(Pi*i*x/Lx)*Sin(Pi*j*y/Ly)*Sin(Pi*k*z/Lz)
       End Do
       c(i, j, k) = 8*r / (Pi**2*(Ly*Lz/Lx*i**2 + Lx*Lz/Ly*j**2 + Lx*Ly/Lz*k**2))
   End Do

Второй этап:

   Do i0 j0 k0 = 1 1 1, M_x M_y M_z
       x = i0*dx
       y = j0*dy
       z = k0*dz
       s = 0.0D0
       Do i j k = 1 1 1, N_x N_y N_z
           s = s + c(i, j, k)*Sin(Pi*i*x/Lx)*Sin(Pi*j*y/Ly)*Sin(Pi*k*z/Lz)
       End Do
       u(i0, j0, k0) = s
   End Do

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

Вычислительную сложность первого варианта алгоритма можно оценить как [math]\displaystyle{ O(N_x N_y N_z M_x M_y M_z) }[/math] операций с плавающей точкой.