Ctrl+K

Теория метода конечных элементов (FEM)

Теория метода конечных элементов (FEM)#

Мотивация#

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

Нужен более универсальный метод - таким и является Finite Element Method.

Общая формулировка метода конечных элементов (FEM)#

Постановка математической задачи#

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

\[ \hat{L}(u)=f(\vec x) \]

где функция \(f(\vec x)\) - не зависит от искомого решения \(u\), а оператор \(\hat{L}\) в общем случае может быть даже нелинейным. Решение необходимо найти на области \(\Omega\).

Примеры:

  • Линейный оператор по координате (Задача Пуассона)

\[ \hat{L}(u) = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \]

  • Линейный оператор c эволюцией по времени с переменными коэффициентами

\[ \hat{L}(u) = \frac{\partial u}{\partial t} - \sin y \cdot \frac{\partial^2 u}{\partial x^2} - x\cdot\frac{\partial^2 u}{\partial y^2} \]

  • Линейный оператор с высшими производными и переменными коэффициентами

\[ \hat{L}(u) = \frac{\partial^3 u}{\partial x^3} - \frac{\partial^4 u}{\partial y^2 \partial x^2} \]

  • Нелинейный оператор

\[ \hat{L}(u) = \Delta u+u-u^3 \]

Алгоритм Finite Element Method#

Опишем «скелет» метода FEM:

  1. Искомую функцию записываем в виде линейной комбинации известных функций (условный базис) \(\{\phi_i(x)\}_{i=1}^N\):

\[ u(\vec x) \approx u^N(\vec x)=\sum_{i=1}^N u_i \cdot \phi_i(\vec x), \quad u_i \in \mathbb R \]

  1. Подставляем \(u^N(\vec x)\) в исходное уравнение и получаем невязку:

\[ r^N(\vec x)=\hat L \left(u^N\right) - f(\vec x) \]

  1. Записываем условие ортогональности невязки и тестовых функций - так называемую слабую форму. Для простоты в качестве тестовых функций возьмём наш «базис» (условие Голёркина - часто приводит к нефизичным осцилляциям решения):

\[ \int_{\Omega} r^N(\vec{x}) \cdot \phi_i(\vec{x}) \cdot d \Omega=0, \quad i \in \overline{1, \ldots, N} \]

  1. Интегрируя, получаем \(N\) алгебраических уравнений относительно неизвестных коэффициентов разложения \(u_i\). Если исходное УРЧП линейно, то и эта система будет линейной.

Примечание. Основное предположение метода заключается в том, что если невязка \(r^N(\vec x)\) ортогональна всем векторам из хорошего базиса, то она обязана быть близка к нулю \(\rightarrow\) мы нашли решение практически удовлетворяющее исходному УРЧП.

Обратите внимание, что в данном подходе нигде явно не фигурирует дискретизация - она как бы «зашита» в выбор базисных функций.

Базисные функции не обязаны быть ортогональными, главное, чтобы коэффициенты разложения определялись единственным образом. Пример базиса в 1D - треугольники на сетке:

Basis

Image

$\(\phi_j(x)= \begin{cases}0 & \text { if } x \notin\left[x_{j-1}, x_{j+1}\right] \\ \frac{x-x_{j-1}}{h} & \text { if } x \in\left[x_{j-1}, x_j\right] \\ \frac{x_{j+1}-x}{h} & \text { if } x \in\left[x_j, x_{j+1}\right]\end{cases}\)$

image

Для любой сеточно-заданной функции \(u^h\) такой набор функций будет базисом (на краю можно добавить половинку треугольника), хотя строго говоря они относятся к разным пространствам: \(u^h\) - к сеточным функциям, \(\phi_i(x)\) - к непрерывным над \(\Omega\).

Определение. Сетка или mesh - это представление области \(\Omega\) в виде конечного набора однотипных простых множеств. Ранее мы рассматривали только прямоугольные сетки одинакового размера, но в FEM удобнее использовать неравномерные треугольные/прямоугольные сетки для более точного описания геометрии \(\Omega\).

original image

В общем случае (на произвольном mesh), базисные функции будут похожи на одномерный случай - главное выполнение следующего свойства:

\[\begin{split} \phi_j\left(x_i\right)= \begin{cases}1 & i=j \\ 0 & i \neq j\end{cases} \end{split}\]

Таким образом, количество базисных функций должно равняться числу узлов, на которых ищем сеточную функцию \(u^h\).

Определение. Конечный элемент или finite element - это тройка \(\left(K, \, P_K, \, \Sigma\right)\) где

  • \(K\) - элемент сетки (ex. треугольничек или интервал);

  • \(P_K\) - конечномерное линейное пространство полиномов, определенное на \(K\) (две прямые в примере выше);

  • \(\Sigma\) - набор степеней свободы (DoF), т.е. значения полинома в вершинах \(K\).

Эти три составляющие определяют единственный полином на элементе \(K\).

Пример. Как-то так выглядит решение одномерного уравнения Пуассона с \(f=-1\) при разных сетках \(h\). Ясно видим отличительную особенность метода FEM - вне сетки решение получается нефизичным, но в узлах сетки сходимость крайне быстрая.

original image

Пример. То же самое уравнение, но с другим базисом. На этот раз используются параболы. За деталями см. дополнительные материалы.

original image

Аппроксимация и устойчивость#

Аппроксимация

Если базисные функции - полиномы р степени, а решение нашей задачи \(\mathrm{p}+1\) раз дифференциируемо, то

\[ \left\|u-u_h\right\|_{L^2(\Omega)} \leqslant c \cdot h^{p+1}\cdot \|u\|_{H^{p+1}(\Omega)}, \]

где

\[ \|f\|_{H^p}=\left(\int_{\Omega}|f|^2+|\nabla f|^2+\ldots+\left|\nabla^p f\right|^2\right)^{\frac{1}{p}} \]

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

В случае с уравнением Пуассона получаем, что \(\mathrm{p}=1\) и

\[ \left\|u-u_h\right\|_{L^2(\Omega)} \leqslant c \cdot h^2 \cdot \|u\|_{H^2(\Omega)}, \]

Устойчивость

Устойчивость решения задачи определяется устойчивостью решения системы линейных уравнений. При плохо выбранном базисе мы получим плохо обусловленную матрицу.

Содержание