Разрывный метод Галёркина

Разрывный метод Галёркина (англ.  discontinuous Galerkin method, сокращенно DGM) — метод решения операторных уравнений, в основном дифференциальных уравнений. Является развитием классического метода конечных элементов (МКЭ), основанного на вариационной постановке Галёркина.

История метода

Разрывный метод Галёркина был впервые предложен и в начале 70-ы годов XX века, как метод решения дифференциальных уравнений в частных производных, в 1973 году Ридом (англ.) и Хиллом был предложен вариант метода для решения гиперболического уравнения переноса нейтронов. Первая формулировка метода для решения эллиптических задач не может быть определена единственной публикацией, однако на развитие метода оказали сильное влияние Иво Бабушка (англ.) и Жак-Луи Лионс (англ.). Для уравнений 4 порядка вариант метода представил Бэйкер в 1977 году. Так же своим развитием метода обязан публикациям Арнольди, Бреззи, Кокбёрна и Марини.

Основные определения

Конечным элементом [math]\displaystyle{ \Sigma }[/math] называется тройка пространств [math]\displaystyle{ (K,~ P, ~L) }[/math] , где:

• [math]\displaystyle{ K }[/math] — геометрический конечный элемент — множество точек пространства (часто, говоря «конечный элемент», подразумевают именно это множество);
• [math]\displaystyle{ P }[/math] — множество координатных функций — базис пространства, которое будет аппроксимировать решение;
• [math]\displaystyle{ L }[/math] — множество степеней свободы — базис сопряженного пространства. [math]\displaystyle{ L = P^* }[/math].

Идея метода и отличие от МКЭ

Рассмотрим идею метода для решения ДУ второго порядка в области [math]\displaystyle{ \Omega }[/math]. В отличие от метода Галёркина, где выполняется постановка в слабой форме, в DGM выполняется постановка в слабой слабой форме (ультра слабой форме, англ. ultra weak variational formulation).
Представим исходное уравнение в виде двух уравнений первого порядка. В зависимости от природы уравнений это можно сделать несколькими способами, которые приведут к различным вариационным постановкам. Далее на расчётной области строим сетку [math]\displaystyle{ K = \left \{ K_i \right \}: K_i \cap K_j = \varnothing, \Omega = \bigcup K_i }[/math], выполняем вариационную постановку Галёркина для каждой подобласти при этом будет использоваться четыре пространства: два пространства(координатное и проекционное) для самой функции и два для её производной. После этого уравнения суммируются по всей области и из получившихся системы из двух уравнений каким-нибудь способом исключается одно.
Данное описание является весьма общим и неоднозначным, поскольку метод всегда подстраивается под конкретные задачи и получение ультра слабой вариационной постановки зависит от природы процесса и цели решения уравнения.

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

Примеры для конкретных видов уравнений

Уравнение теплопроводности

Рассмотрим простейший случай стационарного уравнения теплопроводности:

[math]\displaystyle{ -\Delta u = f~ in~ \Omega }[/math]
[math]\displaystyle{ u = g ~on ~\partial \Omega }[/math]

[math]\displaystyle{ \lambda }[/math] — коэффициент теплопроводности, [math]\displaystyle{ f = f(x,y,z) }[/math] — правая часть уравнения.
Выполним замену [math]\displaystyle{ \nabla u = \sigma }[/math] и тем самым сведём уравнение второго порядка к двум уравнениям первого порядка:

[math]\displaystyle{ \left \{ \begin{array}{lcr} \nabla u = \sigma \\ - \nabla \cdot \sigma = f \\ \end{array} \right. }[/math]

На расчётной области введём пространство Лебега [math]\displaystyle{ L_2 (\Omega) }[/math] с соответствующим ему скалярным произведением: [math]\displaystyle{ (u, v) = \int_\Omega uv d\Omega }[/math]. И соответствующие ему конечноэлементные пространства:
[math]\displaystyle{ V_h = \left \{ v \in L_2 (\Omega) : v \in P_i, \forall K_i \in K \right \} }[/math] — пространство скалярных функций, для аппроксимации решения
[math]\displaystyle{ \Sigma_h = \left \{ \tau \in [L_2 (\Omega)]^3 : \tau \in \Sigma_i, \forall K_i \in K \right \} }[/math] — пространство векторных функций для аппроксимации градиента решения
Введённые пространства являются пространствами Соболева (скалярным и векторным) с соответствующей нормой. Из этих пространств выберем тестовые функции [math]\displaystyle{ v, \tau }[/math] и для каждого уравнения выполним постановку Галёркина на отдельном элементе, получим систему уравнений в слабой форме:

[math]\displaystyle{ \left \{ \begin{array}{lcr} \int_{K_i} \sigma \cdot \tau dx \ = - \int_{K_i} u \nabla \cdot \tau dx \ + \int_{\partial K_i} \hat u_{K_i} (n_{K_i} \cdot \tau) ds\ ~\forall \tau \in \Sigma_i\\ \int_{K_i} \sigma \cdot \nabla v dx\ = \int_{K_i} fvdx\ + \int_{\partial K_i} (\hat \sigma_{K_i} \cdot n_{K_i}) v ds\ ~\forall v \in V_i\\ \end{array} \right. }[/math]

Функции [math]\displaystyle{ \hat u, \hat \sigma }[/math] — это численные потоки, которые могут быть определены по-разному (что ведёт к различным методам) и должны удовлетворять следующим условиям:

• Численные потоки однозначны на границе [math]\displaystyle{ \partial \Omega }[/math]
• Численные потоки представимы в виде: [math]\displaystyle{ \hat u (v) = \Bigl. v \Bigr|_{\partial \Omega},~ \hat \sigma (v, \nabla v) = \Bigl. \nabla v \Bigr|_{\partial \Omega} }[/math], где [math]\displaystyle{ v }[/math] — гладкая функция удовлетворяющая заданным первым краевым условиям.

Для упрощения записи вводят оператор среднего и оператор скачка, определяющие поведение функций на границе элементов:

Операторы среднего и скачка
Оператор среднего	Оператор скачка	Область действия
[math]\displaystyle{ \left \langle v \right \rangle = (v_i + v_j)/2 }[/math]	[math]\displaystyle{ [v] = v_i n_i + v_j n_j }[/math]	[math]\displaystyle{ \Gamma_{ij} = \partial K_i \cap \partial K_j }[/math]
[math]\displaystyle{ \left \langle v \right \rangle = v_i }[/math]	[math]\displaystyle{ [v] = v_i n_i }[/math]	[math]\displaystyle{ \partial \Omega }[/math]
[math]\displaystyle{ \left \langle \tau \right \rangle = (\tau_i + \tau_j)/2 }[/math]	[math]\displaystyle{ [\tau] = \tau_i \cdot n_i + \tau_j \cdot n_j }[/math]	[math]\displaystyle{ \Gamma_{ij} }[/math]
[math]\displaystyle{ \left \langle \tau \right \rangle = \tau_i }[/math]	[math]\displaystyle{ [\tau] = \tau_i \cdot n_i }[/math]	[math]\displaystyle{ \partial \Omega }[/math]

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

[math]\displaystyle{ \left \{ \begin{array}{lcr} \int_{\Omega} \sigma \cdot \tau dx \ = - \int_{\Omega} u \nabla \cdot \tau dx \ + \sum_{K_i \in K}\int_{\partial K_i} \hat u_{K_i} (n_{K_i} \cdot \tau) ds\ ~\forall \tau \in \Sigma_h\\ \int_{\Omega} \sigma \cdot \nabla v dx\ = \int_{\Omega} fvdx\ + \sum_{K_i \in K} \int_{\partial K_i} (\hat \sigma_{K_i} \cdot n_{K_i}) v ds\ ~\forall v \in V_h\\ \end{array} \right. }[/math]

Воспользуемся свойством: [math]\displaystyle{ \sum_{K_i} \int_{\partial K_i} v(\tau \cdot n_{K_i}) ds \ = \sum_{K_i} \int_{\partial K_i} ([v] \cdot \left \langle \tau \right \rangle + \left \langle v \right \rangle [\tau]) ds \ }[/math] и в результате получим ультра слабую вариационную постановку для исходного уравнения:

[math]\displaystyle{ \int_{\Omega} \nabla u \cdot \nabla v dx\ + \sum_{K_i} \int_{\partial K_i}{ [\hat u - u] \cdot \left \langle \nabla v \right \rangle ds} + \sum_{K_i} \int_{\partial K_i}{\left \langle \hat u - u \right \rangle [\nabla v] ds} = \int_{\Omega} {fv dx} + \sum_{K_i} \int_{\partial K_i}{ [v] \left \langle \hat \sigma \right \rangle ds} + \sum_{K_i} \int_{\partial K_i}{ \left \langle v \right \rangle [\sigma] ds},~\forall v \in V_h }[/math]

Осталось определить численные потоки. Определение численных потоков связанно с задачей и требованиями к решению и ведёт к различным методам, например:

Функция и область	IP-метод	Стабилизированный IP-метод	NIPG
[math]\displaystyle{ \hat u }[/math] на [math]\displaystyle{ \partial \Omega }[/math]	[math]\displaystyle{ \hat u = 0 }[/math]
[math]\displaystyle{ \hat u }[/math] на [math]\displaystyle{ \partial K_i }[/math]	[math]\displaystyle{ \hat u = \langle u \rangle }[/math]		[math]\displaystyle{ \hat u = \langle u \rangle + n_{K_i} \cdot [u] }[/math]
[math]\displaystyle{ \hat \sigma }[/math] на [math]\displaystyle{ \partial \Omega }[/math] и на [math]\displaystyle{ \partial K_i }[/math]	[math]\displaystyle{ \hat \sigma = \langle \nabla u \rangle }[/math]	[math]\displaystyle{ \hat \sigma = \left \langle \nabla u \right \rangle + c[u], ~ c- const }[/math]

Уравнения Максвелла в гармоническом режиме

Подход к построению ультра слабой вариационной постановки для уравнений Максвелла может быть различным: систему уравнений первого порядка можно получить напрямую из самих уравнений Максвелла или сведя эти уравнения к уравнению Гельмгольца, а затем сделав замену, аналогичную замене для уравнения теплопроводности, получить систему первого порядка. В данном случае мы воспользуемся первым способом. Система уравнений Максвелла в гармоническом режиме с частотой [math]\displaystyle{ \omega }[/math], в одном из простейших случаев выглядит как:

[math]\displaystyle{ \left \{ \begin{array}{lcr} -i\omega \epsilon \mathbf{E} - \nabla \times \mathbf{H} = 0\\ -i\omega \mu \mathbf{H} + \nabla \times \mathbf{E} = 0 \\ \end{array} \right. }[/math]

Оба уравнения выполняются в расчётной области [math]\displaystyle{ \Omega }[/math]. Краевые условия: [math]\displaystyle{ \mathbf{E} \times \mathbf{n} = -2\mathbf{g} }[/math]. Домножим скалярно оба уравнения на пробные функции [math]\displaystyle{ \mathbf{\xi^{K_i}}, \mathbf{\psi^{K_i}} }[/math], определённые на соответствующем элементе [math]\displaystyle{ K_i \in K }[/math]. Функции из этого же пространства будут использоваться как базисные. Для их определения используем присоединённую систему уравнений Максвелла:

[math]\displaystyle{ \left \{ \begin{array}{lcr} i\omega \epsilon \mathbf{\xi^{K_i}} + \nabla \times \mathbf{\psi^{K_i}} = 0\\ i\omega \mu \mathbf{\psi^{K_i}} - \nabla \times \mathbf{\xi^{K_i}} = 0 \\ \end{array} \right. }[/math]

Оба уравнения этой системы записаны для одного элемента [math]\displaystyle{ K_i }[/math]. Домножив каждое уравнение системы на пробную функцию, преобразовав их, используя аналог формулы Грина и сложив, получим следующее выражение:

[math]\displaystyle{ \int_{K_i} { \left( \mathbf {E}\cdot \overline{i\omega \epsilon \mathbf{\xi^{K_i}} + \nabla \times \mathbf{\psi^{K_i}}} + \mathbf{H} \cdot \overline{i\omega \mu \mathbf{\psi^{K_i}} - \nabla \times \mathbf{\xi^{K_i}}} \right) dx} = \int_{\partial K_i} { \left( \mathbf{n^{K_i}} \times \mathbf{H} \cdot \mathbf{\xi^{K_i}} - \mathbf{n^{K_i}} \times \mathbf{E} \cdot \mathbf{\psi^{K_i}} \right) ds} }[/math]

Учитывая систему уравнений для пробных функций данное выражение упрощается до:

[math]\displaystyle{ \int_{\partial K_i} { \left( \mathbf{n^{K_i}} \times \mathbf{H} \cdot \mathbf{\xi^{K_i}} - \mathbf{n^{K_i}} \times \mathbf{E} \cdot \mathbf{\psi^{K_i}} \right) ds} = 0 }[/math]

Введём обозначения:

Вектора	Матрицы
[math]\displaystyle{ \mathbf{u^{K_i}} = \binom{\Bigl.\mathbf{E} \Bigr|_{K_i}}{\Bigl.\mathbf{H} \Bigr|_{K_i}} ~~\mathbf{\phi^{K_i}} = \binom{\mathbf{\xi^{K_i}}}{\mathbf{\psi^{K_i}}} }[/math]	[math]\displaystyle{ Z^{K_i} = \begin{pmatrix} 0 & n_3^{K_i} & -n_2^{K_i} \\ -n_3^{K_i} & 0 & n_1^{K_i} \\ n_2^{K_i} & - n_1^{K_i} & 0 \end{pmatrix} }[/math]
[math]\displaystyle{ \chi^{K_i} = \Bigl. L^{K_i, +} \mathbf{u^{K_i}} \Bigr|_{\partial K_i} ~~ \gamma^{K_i} = \Bigl. L^{K_i, +} \mathbf{\phi^{K_i}} \Bigr|_{\partial K_i} }[/math] [math]\displaystyle{ F_{K_i}(\gamma^{K_i}) = \Bigl. - L^{K_i, -} \mathbf{\phi^{K_i}} \Bigr|_{\partial K_i} }[/math]	[math]\displaystyle{ D^{K_i} \begin{pmatrix} 0 & (Z^{K_i})^T \\ Z^{K_i} & 0 \end{pmatrix} }[/math] [math]\displaystyle{ L^{K_i, \pm} = \frac{1}{\sqrt{2\sigma}}\left( Z^{K_i},~ \pm \sigma (Z^{K_i})^2 \right) }[/math]

Теперь задача ставится как нахождение векторов [math]\displaystyle{ \chi^{K_i} }[/math] для всех элементов, удовлетворяющих следующим уравнениям :

[math]\displaystyle{ \int_{\partial K_i} {\chi^{K_i} \cdot \overline{\gamma^{K_i}} ds} - \sum_{\Gamma_{ij} = \partial K_i \cap \partial K_j \neq \varnothing} \int_{\Gamma_{ij}} { \chi^{K_j} \cdot \overline{F_{K_i}(\gamma^{K_i})}ds} - \sum_{\Gamma_{i} = \partial K_i \cap \partial \Omega \neq \varnothing} \int_{\Gamma_{i}} { \chi^{K_i} \cdot \overline{F_{K_i}(\gamma^{K_i})}ds} = \int_{\partial K_i} {\mathbf{g} \cdot \overline{F_{K_i}(\gamma^{K_i})} ds} ~~ \forall \gamma^{K_i} }[/math]

При наличии у исходных уравнений правой части в итоговой ультра слабой постановке появились бы дополнительные слагаемые в виде интегралов по самому конечному элементу. Особенностями метода является то, что после получения решения системы необходимо решить ещё одну, чтобы получить вектор [math]\displaystyle{ \mathbf{u} }[/math], однако найдя его мы сразу узнаём значениях двух составляющих электромагнитного поля: [math]\displaystyle{ \mathbf{E} }[/math] и [math]\displaystyle{ \mathbf{H} }[/math]. Данную постановку можно ещё преобразовать, получив сразу уравнение относительно вектора [math]\displaystyle{ \mathbf{u} }[/math].

Литература

• Arnold D.N., Brezzi F., Cocburn B., Mariri D. «Unified analysis of DCM for elliptic problems», 2002