Пример. Уравнение Пуассона

Примечание. Следующий код рекомендуется выполнять на google colab - феникс достаточно прожорливый по месту на диске + скорее всего его невозможно запустить на windows (гуглите wsl).

\[\begin{split} \begin{aligned} -\Delta u=f & \text { in } \Omega, \\ u=u_0 & \text { on } \partial \Omega . \end{aligned} \end{split}\]

1. Перевести уравнение в слабую форму

Давайте домножим уравнение на функцию \(v\), которая равна нулю на границе, и проинтегрируем его. Эта функция называется тестовой.

Определение 1. \(H_1(\Omega)\) - Пространство функций с конечным интегралом от \(v^2 u|\nabla v|^2\).

Единичка внизу говорит про количество производных, интеграл от квадрата которых конечен.

Определение 2. Тестовая функция - \(v \in H_1(\Omega): v(\partial \Omega)=0\).

\[ -\int_{\Omega}(\Delta u) v \mathrm{~d} x=\int_{\Omega} f v \mathrm{~d} x \]

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

\[ -\int_{\Omega}(\Delta u) v \mathrm{~d} x=\int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x-\int_{\partial \Omega} \frac{\partial u}{\partial n} v \mathrm{~d} s \]

Последний член занулятеся, так как \(\mathrm{v}=0\) на границе. И мы приходим к уравнению в слабой форме.

\[ \int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x=\int_{\Omega} f v \mathrm{~d} x \]

Определение 3. Слабая форма - уравнение вида \(a(u, v)=L(v)\), где \(a\) - билинейная форма, \(L\) - линейная форма.

Наша основная задача - получить уравнение в таком виде. Каждая из этих функций (\(u\), \(v\)) принадлежит бесконечномерному пространству функций с конечным интегралом от квадрата производной. Компьютер умеет работать с конечными простраствами, поэтому решения мы будем искать на конечномерных подпространствах. Но как же их выбрать?

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

Однако такое есть не всегда.

2. Выбрать пространство функций, на котором ищется ответ

Начиная с этого пункта, за нас всё будет делать fenics. Приводим сами вычисления только для понимания происходящего.

Далее мы будем рассматривать одномерную задачу на отрезке \([0,1]\), где \(\mathrm{f}=1, \; \mathrm{u}(0)=\mathrm{a}, \; \mathrm{u}(1)=\mathrm{b}\).

Разобьем отрезок на \(\mathrm{N}+2\) части. И выберем простой набор функций, так чтобы как можно больше из них были ортогональны. Ортогональность значительно упрощает вычисления, и это мы увидим ниже.

Пусть \(\phi_j\) будет линейной функцией на отрезках \([(\mathrm{j}-1) \mathrm{h}, \mathrm{jh}]\) и \([\mathrm{jh},(\mathrm{j}+1) \mathrm{h}]\), где \(\mathrm{h}=1 /(\mathrm{N}+1)\). Отдельно рассмотрим \(\phi_0\) и \(\phi_{N+1}\). Эти функции нужны, чтобы задать начальные условия.

Итак, мы ищем решение в виде \(u^h=\sum_{j=0}^{N+1} u_j \phi_j\). Буква \(\mathrm{h}\) значит, что это решение для конечномерного пространства.

Давайте зададим пространство пробных функций почти также. Мы знаем, что всякая тестовая функция зануляется на границе, поэтому коэффициенты перед \(\phi_0, \phi_{N+1}\) равны 0 . Значит, \(v^h=\sum_{j=1}^N v_j \phi_j\)

Вспомним нашу задачу.

\[ \int_{[0,1]} \nabla u \cdot \nabla v \mathrm{~d} x=\int_{[0,1]} f v \mathrm{~d} x . \]

Скалярное произведение линейно по каждому аргументу, поэтому нам не надо проверять его для произвольной функции \(v\). Достаточно убедится, что оно выполняется на базисных функциях.

3. Составить и решить систему линейных уравнений на коэффициенты

\[\begin{split} \begin{aligned} & \left(u_h^{\prime}, \phi_h^{\prime}\right)=\left(f, \phi_h\right) \quad \forall \phi_h \\ & \Leftrightarrow\left(\left(\sum_{j=1}^n u_j \phi_j\right)^{\prime}, \phi_h^{\prime}\right)=\left(f, \phi_h\right) \quad \forall \phi_h \\ & \Leftrightarrow \sum_{j=1}^n u_j\left(\phi_j^{\prime}, \phi_h^{\prime}\right)=\left(f, \phi_h\right) \quad \forall \phi_h \\ & \Leftrightarrow \sum_{j=1}^n u_j\left(\phi_j^{\prime}, \phi_i^{\prime}\right)=\left(f, \phi_i\right) \quad \forall 1 \leq i \leq n \\ & \underbrace{\left(\begin{array}{cccc} \left(\phi_1^{\prime}, \phi_1^{\prime}\right) & \left(\phi_1^{\prime}, \phi_2^{\prime}\right) & \cdots & \left(\phi_1^{\prime}, \phi_n^{\prime}\right) \\ \left(\phi_2^{\prime}, \phi_1^{\prime}\right) & \left(\phi_2^{\prime}, \phi_2^{\prime}\right) & \cdots & \left(\phi_2^{\prime}, \phi_n^{\prime}\right) \\ \vdots & \vdots & \ddots & \vdots \\ \left(\phi_n^{\prime}, \phi_1^{\prime}\right) & \left(\phi_n^{\prime}, \phi_2^{\prime}\right) & \cdots & \left(\phi_n^{\prime}, \phi_n^{\prime}\right) \end{array}\right)}_{=A} \underbrace{\left(\begin{array}{c} u_1 \\ u_2 \\ \vdots \\ u_n \end{array}\right)}_{=U}=\underbrace{\left(\begin{array}{c} \left(f, \phi_1\right) \\ \left(f, \phi_2\right) \\ \vdots \\ \left(f, \phi_n\right) \end{array}\right)}_{=F} \\ & \end{aligned} \end{split}\]

Заметим, что здесь нет \(\phi_0, \phi_{N+1}\), потому что они не принадлежат пространству тестовых функций. Они определяются с помощью начальных условий.

\[\begin{split} \begin{aligned} \int_{\Omega} \phi_i^{\prime}(x) \phi_i^{\prime}(x) d x & =\int_{x_{i-1}}^{x_{i+1}} \phi_i^{\prime}(x) \phi_i^{\prime}(x) d x \\ = & \int_{x_{i-1}}^{x_i} \frac{1}{h^2} d x+\int_{x_i}^{x_{i+1}}\left(-\frac{1}{h}\right)^2 d x \\ =\frac{h}{h^2} & +\frac{h}{h^2}=\frac{2}{h} \\ \int_{\Omega} \phi_i^{\prime}(x) \phi_{i \pm 1}^{\prime}(x) d x & =\int_{x_{i-1}}^{x_{i+1}} \phi_i^{\prime}(x) \phi_{i \pm 1}^{\prime}(x) d x \\ & =\int_{x_i}^{x_{i+1}}\left(-\frac{1}{h}\right) \cdot \frac{1}{h} d x \\ & =-\frac{h}{h^2}=-\frac{1}{h} \end{aligned} \end{split}\]

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

Добавим уравнения на \(\phi_0, \phi_{N+1}\). В итоге получим:

\[\begin{split} \begin{aligned} & A=\frac{1}{h}\left(\begin{array}{ccccc} 1 & 0 & & & 0 \\ -1 & 2 & -1 & & \\ & \ddots & \ddots & \ddots & \\ & & -1 & 2 & -1 \\ 0 & & & 0 & 1 \end{array}\right) \\ & F=h \cdot(a, 1, \ldots, 1, b)^T \end{aligned} \end{split}\]

# Данный код протестирован только в Google Colab
# Уточните, как устанавливать FEniCS на вашу систему

try:
    import google.colab  # noqa: F401
except ImportError:
    import ufl_legacy
    import dolfin
else:
    try:
        import ufl_legacy
        import dolfin
    except ImportError:
        !wget "https://fem-on-colab.github.io/releases/fenics-install-real.sh" -O "/tmp/fenics-install.sh" && bash "/tmp/fenics-install.sh"
        import ufl_legacy
        import dolfin

import fenics as fe
import matplotlib.pyplot as plt
import numpy as np

N_ELEMENTS = 20
mesh = fe.UnitSquareMesh(N_ELEMENTS, N_ELEMENTS) #Определяем сетку на которой будем искать решение

lagr_polinomial_first_order = fe.FunctionSpace(
        mesh,
        "Lagrange",
        2,
)#Определяем пространство функций. Оно общее для тестовых функций и для нашего приближенного решения.
#Здесь я использую линейные функции, которые описал в файле.

#зададим начальное условие.
u0 = fe.Expression("1 + 0.2*x[0]*x[0] + 0.1*x[1]*x[1]", degree = 1)
u0 = fe.interpolate(u0, lagr_polinomial_first_order)

#Эта функция нужна для скрытых методов сетки. on_boundary булева переменная,
#которая приходит из автоматической проверки точки на принедлежность границе.
#можно переписать функию и не использовать on_boundary.
def on_boundary(x, on_boundary):
    return on_boundary

#Граничные условия
bc = fe.DirichletBC(
    lagr_polinomial_first_order,
    u0,
    on_boundary,
)

#задаем функции для составления основного уравнения
u = fe.TrialFunction(lagr_polinomial_first_order)
v = fe.TestFunction(lagr_polinomial_first_order)
u_solution = fe.Function(lagr_polinomial_first_order)


#Определим f
sigma = 0.3
x0 = 0.2
y0 = 0.2

f = fe.Expression("7*exp(-0.5*(pow((x[0] - x0)/sigma, 2)) "\
                    " - 0.5*(pow((x[1] - y0)/sigma, 2)))",
                   x0=x0, y0=y0, sigma=sigma, degree = 2)

#Замечу, что константы определяются так, а не через expression. Можно попробовать подставить f1 в решение.
f1 = fe.Constant("0.0")


#составим уравнение
A = fe.dot(fe.grad(u), fe.grad(v)) * fe.dx
L = f * v * fe.dx

fe.solve(A == L, u_solution, bc)

pl = fe.plot(u_solution)
#fe.plot(mesh) #сетка
plt.colorbar(pl)
plt.show()
#fe.interactive()