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

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

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

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

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

Конечным элементом $\Sigma$ называется тройка пространств $(K,~P,~L)$ , где:

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

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

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

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

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

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

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

$-\Delta u=f~in~\Omega$
$u=g~on~\partial \Omega$

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

$\left\{{\begin{array}{lcr}\nabla u=\sigma \\-\nabla \cdot \sigma =f\\\end{array}}\right.$

На расчётной области введём пространство Лебега $L_{2}(\Omega )$ с соответствующим ему скалярным произведением: $(u,v)=\int _{\Omega }uvd\Omega$. И соответствующие ему конечноэлементные пространства:
$V_{h}=\left\{v\in L_{2}(\Omega ):v\in P_{i},\forall K_{i}\in K\right\}$ — пространство скалярных функций, для аппроксимации решения
$\Sigma _{h}=\left\{\tau \in [L_{2}(\Omega )]^{3}:\tau \in \Sigma _{i},\forall K_{i}\in K\right\}$ — пространство векторных функций для аппроксимации градиента решения
Введённые пространства являются пространствами Соболева (скалярным и векторным) с соответствующей нормой. Из этих пространств выберем тестовые функции $v,\tau$ и для каждого уравнения выполним постановку Галёркина на отдельном элементе, получим систему уравнений в слабой форме:

$\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 vdx\ =\int _{K_{i}}fvdx\ +\int _{\partial K_{i}}({\hat {\sigma }}_{K_{i}}\cdot n_{K_{i}})vds\ ~\forall v\in V_{i}\\\end{array}}\right.$

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

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

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

Операторы среднего и скачка
Оператор среднего	Оператор скачка	Область действия
$\left\langle v\right\rangle =(v_{i}+v_{j})/2$	$[v]=v_{i}n_{i}+v_{j}n_{j}$	$\Gamma _{ij}=\partial K_{i}\cap \partial K_{j}$
$\left\langle v\right\rangle =v_{i}$	$[v]=v_{i}n_{i}$	$\partial \Omega$
$\left\langle \tau \right\rangle =(\tau _{i}+\tau _{j})/2$	$[\tau ]=\tau _{i}\cdot n_{i}+\tau _{j}\cdot n_{j}$	$\Gamma _{ij}$
$\left\langle \tau \right\rangle =\tau _{i}$	$[\tau ]=\tau _{i}\cdot n_{i}$	$\partial \Omega$

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

$\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 vdx\ =\int _{\Omega }fvdx\ +\sum _{K_{i}\in K}\int _{\partial K_{i}}({\hat {\sigma }}_{K_{i}}\cdot n_{K_{i}})vds\ ~\forall v\in V_{h}\\\end{array}}\right.$

Воспользуемся свойством: $\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\$ и в результате получим ультра слабую вариационную постановку для исходного уравнения:

$\int _{\Omega }\nabla u\cdot \nabla vdx\ +\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 }{fvdx}+\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}$

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

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

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

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

$\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.$

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

$\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.$

Оба уравнения этой системы записаны для одного элемента $K_{i}$. Домножив каждое уравнение системы на пробную функцию, преобразовав их, используя аналог формулы Грина и сложив, получим следующее выражение:

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

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

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

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

Вектора	Матрицы
$\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}}} }}$	$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}}$
$\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}}$ $F_{K_{i}}(\gamma ^{K_{i}})={\Bigl .}-L^{K_{i},-}\mathbf {\phi ^{K_{i}}} {\Bigr |}_{\partial K_{i}}$	$D^{K_{i}}{\begin{pmatrix}0&(Z^{K_{i}})^{T}\\Z^{K_{i}}&0\end{pmatrix}}$ $L^{K_{i},\pm }={\frac {1}{\sqrt {2\sigma }}}\left(Z^{K_{i}},~\pm \sigma (Z^{K_{i}})^{2}\right)$

Теперь задача ставится как нахождение векторов $\chi ^{K_{i}}$ для всех элементов, удовлетворяющих следующим уравнениям :

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

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

Литература

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