Пробные функции

Функция, которая описывает поведение на симплексе, называется пробной, вероятно потому, что на заре развития метода, ввиду слабости вычислительной техники, приходилось пробовать и подбирать функции, которые бы аппроксимировали решение. Существует огромное количество пробных функций, можно использовать и тригонометрические даже, но на практике распространение получили полиномиальные функции, что как бы намекает на разложение решения в ряд Тейлора, по аналогии с разложением в ряд Фурье, за тем лишь исключением, что разложение в ряд Тейлора мы осуществляем не на всей геометрии, а только лишь на симплексах, что позволяет нам использовать максимально простые полиномы. Выбирались, как правило, полиномы 3 степени, с 9-10 коэффициентами, и это требовало добавления дополнительных точек на симплекс, чтобы эти коэффициенты определить, что приводило к 9-10 СЛАУ на каждом симплексе, но решать их нужно было всего один раз, что с точки зрения вычислительной не проблема даже по тем временам. Конечно, это здравый подход, но я решил «попробовать» построить решение с помощью максимально простых полиномов, а именно полиномов вида:

υh(x,y,z)=i=1nqiϕi(x,y,z),\upsilon_h(x,y,z) = \sum_{i=1}^{n} q_i \cdot \phi_i(x, y, z),
(5.5)

где nn — количество узлов, qiq_i — коэффициенты, ϕi(x,y,z)\phi_i(x, y, z) — так называемые функции «крышки», которые равны 1 на узле симплекса и нулю в остальных узлах геометрии, и как раз-таки за счёт того, что функции «крышки» представляют из себя полиномы, пробные функции будут тоже полиномами.

ϕi(x,y,z)=ai+bix+ciy+diz,\phi_i(x,y,z) = a_{i} + b_{i} \cdot x + c_{i} \cdot y + d_{i} \cdot z,
(5.6)

где ii — номер узла, а не степень, ai,bi,ci,dia_{i}, b_{i}, c_{i}, d_{i} — коэффициенты, зависящие от геометрии конкретного узла ii на симплексе. После того, как геометрия будет разбита на симплексы, нужно будет их один раз вычислить для каждого узла симплекса. Функцию «крышку» я записал для трёхмерной системы координат, но дальше будем работать с двумерной или одномерной так же.

Давайте подробнее поговорим про функции «крышки» и начнём с одномерного случая.

Одномерные линейные функции крышки
Рис. 5.4. Функции «крышки» для симплексов-отрезков.

На рисунке Функции «крышки» для симплексов-отрезков приведено пять функций «крышек» ϕ0,,ϕ4\phi_0, …, \phi_4 и четыре симплекса: в общей сложности для одномерного пространства на nn точек приходится nn функций «крышек» и n1n - 1 симплексов соответственно. На краях, в узлах x0x_0 и x4x_4, функции имеют вид половинок «крышек». Примем i=3i = 3, тогда симплекс между точками xi1x_{i-1} и xix_i (на рисунке он выделен заливкой и ограничен пунктирными линиями) определяется следующим образом: υi1(x)=qi1ϕi1(x)+qiϕi(x)\upsilon_{i-1}(x) = q_{i-1} \cdot \phi_{i-1}(x) + q_i \cdot \phi_i(x) за счёт того, что остальные функции «крышки» не оказывают влияния на него ввиду того, что принимают нулевое значение на интервале (xi1,xi)(x_{i-1}, x_i). Функция ϕi1(x)\phi_{i-1}(x) проходит через точки (xi1,1)(x_{i-1}, 1) и (xi,0)(x_i, 0), а функция ϕi(x)\phi_i(x) через точки (xi1,0)(x_{i-1}, 0) и (xi,1)(x_i, 1). Тогда справедливо

{ai1+bi1xi1=1ai1+bi1xi=0{ai+bixi1=0ai+bixi=1\begin{cases} a_{i-1} + b_{i-1} \cdot x_{i-1} = 1\\ a_{i-1} + b_{i-1} \cdot x_i = 0 \end{cases} \begin{cases} a_{i} + b_{i} \cdot x_{i-1} = 0\\ a_{i} + b_{i} \cdot x_i = 1 \end{cases}
(5.7)

Решение систем соответственно

{ai1=xixi1xibi1=1xi1xi{ai=xi1xi1xibi=1xi1xi\begin{cases} a_{i-1} = \frac{\displaystyle -x_i}{\displaystyle x_{i-1} - x_i}\\ b_{i-1} = \frac{\displaystyle 1}{\displaystyle x_{i-1} - x_i} \end{cases} \begin{cases} a_{i} = \frac{\displaystyle x_{i-1}}{\displaystyle x_{i-1} - x_i}\\ b_{i} = \frac{\displaystyle -1}{\displaystyle x_{i-1} - x_i} \end{cases}
(5.8)

Если для одномерных симплексов функции «крышки» имеют вид линий на плоскости, то для двумерных симплексов уже будут поверхности в трёхмерном пространстве.

Рис. 5.5. Функции «крышки» для симплексов-треугольников.

Мы взяли все три точки треугольника t2t_2 с рисунка Симплексы-треугольники и построили три пирамиды, вершинами которых являются «поднятые» на величину 1 вершины исходного треугольника, а основаниями пирамид являются многоугольники с того же рисунка, содержащие внутри себя вершины треугольника. Так как нашему исследованию подвергается треугольник t2t_2 и мы планируем записывать пробную функцию для него, отбросим лишние части и перестроим рисунок, оставив только плоскости, проецируемые на треугольник.

Рис. 5.6. Части функций «крышек» для симплекса-треугольника.

На рисунке Части функций «крышек» для симплекса-треугольника для двумерного пространства, как и для одномерного, каждому из nn узлов сетки соответствует ровно одна функция «крышка», так что всего их nn; число же симплексов-треугольников зависит от конкретной триангуляции. Заметим также, что в пределах одного треугольника отличны от нуля ровно три функции «крышки» — по числу его вершин. Обозначим точки симплекса через (xi,yi),(xi+1,yi+1),(xi+2,yi+2)(x_i, y_i), (x_{i+1}, y_{i+1}), (x_{i+2}, y_{i+2}). В итоге пробная функция примет вид υi(x,y)=qiϕi(x,y)+qi+1ϕi+1(x,y)+qi+2ϕi+2(x,y)\upsilon_i(x, y) = q_i \cdot \phi_i(x, y) + q_{i+1} \cdot \phi_{i+1}(x, y) + q_{i+2} \cdot \phi_{i+2}(x, y). Для определения постоянных ai,bi,cia_{i},b_{i},c_{i}, ai+1,bi+1,ci+1a_{i+1},b_{i+1},c_{i+1}, ai+2,bi+2,ci+2a_{i+2},b_{i+2},c_{i+2} примем во внимание, что координата zz отлична от нуля только на вершинах пирамид, в итоге получим три системы уравнений

{ai+bixi+ciyi=1ai+bixi+1+ciyi+1=0ai+bixi+2+ciyi+2=0{ai+1+bi+1xi+ci+1yi=0ai+1+bi+1xi+1+ci+1yi+1=1ai+1+bi+1xi+2+ci+1yi+2=0{ai+2+bi+2xi+ci+2yi=0ai+2+bi+2xi+1+ci+2yi+1=0ai+2+bi+2xi+2+ci+2yi+2=1\begin{split} &\begin{cases} a_{i} + b_{i} \cdot x_i + c_{i} \cdot y_i = 1\\ a_{i} + b_{i} \cdot x_{i+1} + c_{i} \cdot y_{i+1} = 0\\ a_{i} + b_{i} \cdot x_{i+2} + c_{i} \cdot y_{i+2} = 0 \end{cases}\\ &\begin{cases} a_{i+1} + b_{i+1} \cdot x_i + c_{i+1} \cdot y_i = 0\\ a_{i+1} + b_{i+1} \cdot x_{i+1} + c_{i+1} \cdot y_{i+1} = 1\\ a_{i+1} + b_{i+1} \cdot x_{i+2} + c_{i+1} \cdot y_{i+2} = 0 \end{cases}\\ &\begin{cases} a_{i+2} + b_{i+2} \cdot x_i + c_{i+2} \cdot y_i = 0\\ a_{i+2} + b_{i+2} \cdot x_{i+1} + c_{i+2} \cdot y_{i+1} = 0\\ a_{i+2} + b_{i+2} \cdot x_{i+2} + c_{i+2} \cdot y_{i+2} = 1 \end{cases} \end{split}
(5.9)

Решение уравнений будет соответственно

{ai=(xi+1yi+2xi+2yi+1)/detbi=(yi+1yi+2)/detci=(xi+1+xi+2)/det{ai+1=(xiyi+2+xi+2yi)/detbi+1=(yi+yi+2)/detci+1=(xixi+2)/det{ai+2=(xiyi+1xi+1yi)/detbi+2=(yiyi+1)/detci+2=(xi+xi+1)/det\begin{split} &\begin{cases} a_{i} = (x_{i+1} \cdot y_{i+2} - x_{i+2} \cdot y_{i+1}) / \det\\ b_{i} = (y_{i+1} - y_{i+2}) / \det\\ c_{i} = (-x_{i+1} + x_{i+2}) / \det \end{cases}\\ &\begin{cases} a_{i+1} = ( - x_i \cdot y_{i+2} + x_{i+2} \cdot y_i) / \det\\ b_{i+1} = (-y_i + y_{i+2}) / \det\\ c_{i+1} = (x_i - x_{i+2}) / \det \end{cases}\\ &\begin{cases} a_{i+2} = (x_i \cdot y_{i+1} - x_{i+1} \cdot y_i) / \det\\ b_{i+2} = (y_i - y_{i+1}) / \det\\ c_{i+2} = (-x_i + x_{i+1}) / \det \end{cases} \end{split}
(5.10)

где det=xiyi+1xiyi+2xi+1yi+xi+1yi+2+xi+2yixi+2yi+1\det = x_i \cdot y_{i+1} - x_i \cdot y_{i+2} - x_{i+1} \cdot y_i + x_{i+1} \cdot y_{i+2} + x_{i+2} \cdot y_i - x_{i+2} \cdot y_{i+1} — определитель матрицы.

Сложность работы с трёхмерной системой координат заключается в том, что решение сложно визуализировать, как и функции «крышки», поэтому обойдёмся только формулами. Вообще, в данном руководстве основной акцент будет делаться на двумерные краевые задачи, так как они лишены тривиальности, в отличие от одномерных, и в то же время хорошо визуализируются. Итак, в трёхмерном пространстве тетраэдр служит основанием четырёхмерного пентахорона «крышки». По аналогии с формулами выше смеем предположить, что функция «крышка» будет иметь значение 1 в одной из точек тетраэдра, поднятой в четырёхмерное пространство перпендикулярно основанию и быть нулём в остальных точках тетраэдризации, при этом пробная функция будет иметь вид υh(x,y,z)=q1ϕ1(x,y,z)+q2ϕ2(x,y,z)+q3ϕ3(x,y,z)+q4ϕ4(x,y,z)\upsilon_h(x, y, z) = q_1 \cdot \phi_1(x, y, z) + q_2 \cdot \phi_2(x, y, z) + q_3 \cdot \phi_3(x, y, z) + q_4 \cdot \phi_4(x, y, z). Для определения постоянных ai,bi,ci,dia_{i},b_{i},c_{i},d_{i}, ai+1,bi+1,ci+1,di+1a_{i+1},b_{i+1},c_{i+1},d_{i+1}, ai+2,bi+2,ci+2,di+2a_{i+2},b_{i+2},c_{i+2},d_{i+2}, ai+3,bi+3,ci+3,di+3a_{i+3},b_{i+3},c_{i+3},d_{i+3} примем во внимание, что координата μ\mu четвёртого измерения отлична от нуля только на вершинах пентахоронов, в итоге получим четыре системы уравнений

{ai+bixi+ciyi+dizi=1ai+bixi+1+ciyi+1+dizi+1=0ai+bixi+2+ciyi+2+dizi+2=0ai+bixi+3+ciyi+3+dizi+3=0{ai+1+bi+1xi+ci+1yi+di+1zi=0ai+1+bi+1xi+1+ci+1yi+1+di+1zi+1=1ai+1+bi+1xi+2+ci+1yi+2+di+1zi+2=0ai+1+bi+1xi+3+ci+1yi+3+di+1zi+3=0{ai+2+bi+2xi+ci+2yi+di+2zi=0ai+2+bi+2xi+1+ci+2yi+1+di+2zi+1=0ai+2+bi+2xi+2+ci+2yi+2+di+2zi+2=1ai+2+bi+2xi+3+ci+2yi+3+di+2zi+3=0{ai+3+bi+3xi+ci+3yi+di+3zi=0ai+3+bi+3xi+1+ci+3yi+1+di+3zi+1=0ai+3+bi+3xi+2+ci+3yi+2+di+3zi+2=0ai+3+bi+3xi+3+ci+3yi+3+di+3zi+3=1\begin{split} &\begin{cases} a_{i} + b_{i} \cdot x_i + c_{i} \cdot y_i + d_{i} \cdot z_i = 1\\ a_{i} + b_{i} \cdot x_{i+1} + c_{i} \cdot y_{i+1} + d_{i} \cdot z_{i+1} = 0\\ a_{i} + b_{i} \cdot x_{i+2} + c_{i} \cdot y_{i+2} + d_{i} \cdot z_{i+2} = 0\\ a_{i} + b_{i} \cdot x_{i+3} + c_{i} \cdot y_{i+3} + d_{i} \cdot z_{i+3} = 0 \end{cases}\\ &\begin{cases} a_{i+1} + b_{i+1} \cdot x_i + c_{i+1} \cdot y_i + d_{i+1} \cdot z_i = 0\\ a_{i+1} + b_{i+1} \cdot x_{i+1} + c_{i+1} \cdot y_{i+1} + d_{i+1} \cdot z_{i+1} = 1\\ a_{i+1} + b_{i+1} \cdot x_{i+2} + c_{i+1} \cdot y_{i+2} + d_{i+1} \cdot z_{i+2} = 0\\ a_{i+1} + b_{i+1} \cdot x_{i+3} + c_{i+1} \cdot y_{i+3} + d_{i+1} \cdot z_{i+3} = 0 \end{cases}\\ &\begin{cases} a_{i+2} + b_{i+2} \cdot x_i + c_{i+2} \cdot y_i + d_{i+2} \cdot z_i = 0\\ a_{i+2} + b_{i+2} \cdot x_{i+1} + c_{i+2} \cdot y_{i+1} + d_{i+2} \cdot z_{i+1} = 0\\ a_{i+2} + b_{i+2} \cdot x_{i+2} + c_{i+2} \cdot y_{i+2} + d_{i+2} \cdot z_{i+2} = 1\\ a_{i+2} + b_{i+2} \cdot x_{i+3} + c_{i+2} \cdot y_{i+3} + d_{i+2} \cdot z_{i+3} = 0 \end{cases}\\ &\begin{cases} a_{i+3} + b_{i+3} \cdot x_i + c_{i+3} \cdot y_i + d_{i+3} \cdot z_i = 0\\ a_{i+3} + b_{i+3} \cdot x_{i+1} + c_{i+3} \cdot y_{i+1} + d_{i+3} \cdot z_{i+1} = 0\\ a_{i+3} + b_{i+3} \cdot x_{i+2} + c_{i+3} \cdot y_{i+2} + d_{i+3} \cdot z_{i+2} = 0\\ a_{i+3} + b_{i+3} \cdot x_{i+3} + c_{i+3} \cdot y_{i+3} + d_{i+3} \cdot z_{i+3} = 1 \end{cases} \end{split}
(5.11)

Решение уравнений будет соответственно

{ai=(xi+1yi+2zi+3+xi+1yi+3zi+2+xi+2yi+1zi+3xi+2yi+3zi+1xi+3yi+1zi+2+xi+3yi+2zi+1)/detbi=(yi+1zi+2yi+1zi+3yi+2zi+1+yi+2zi+3+yi+3zi+1yi+3zi+2)/detci=(xi+1zi+2+xi+1zi+3+xi+2zi+1xi+2zi+3xi+3zi+1+xi+3zi+2)/detdi=(xi+1yi+2xi+1yi+3xi+2yi+1+xi+2yi+3+xi+3yi+1xi+3yi+2)/det{ai+1=(xiyi+2zi+3xiyi+3zi+2xi+2yizi+3+xi+2yi+3zi+xi+3yizi+2xi+3yi+2zi)/detbi+1=(yizi+2+yizi+3+yi+2ziyi+2zi+3yi+3zi+yi+3zi+2)/detci+1=(xizi+2xizi+3xi+2zi+xi+2zi+3+xi+3zixi+3zi+2)/detdi+1=(xiyi+2+xiyi+3+xi+2yixi+2yi+3xi+3yi+xi+3yi+2)/det{ai+2=(xiyi+1zi+3+xiyi+3zi+1+xi+1yizi+3xi+1yi+3zixi+3yizi+1+xi+3yi+1zi)/detbi+2=(yizi+1yizi+3yi+1zi+yi+1zi+3+yi+3ziyi+3zi+1)/detci+2=(xizi+1+xizi+3+xi+1zixi+1zi+3xi+3zi+xi+3zi+1)/detdi+2=(xiyi+1xiyi+3xi+1yi+xi+1yi+3+xi+3yixi+3yi+1)/det{ai+3=(xiyi+1zi+2xiyi+2zi+1xi+1yizi+2+xi+1yi+2zi+xi+2yizi+1xi+2yi+1zi)/detbi+3=(yizi+1+yizi+2+yi+1ziyi+1zi+2yi+2zi+yi+2zi+1)/detci+3=(xizi+1xizi+2xi+1zi+xi+1zi+2+xi+2zixi+2zi+1)/detdi+3=(xiyi+1+xiyi+2+xi+1yixi+1yi+2xi+2yi+xi+2yi+1)/det\begin{split} &\begin{cases} a_{i} = (-x_{i+1} \cdot y_{i+2} \cdot z_{i+3} + x_{i+1} \cdot y_{i+3} \cdot z_{i+2} + x_{i+2} \cdot y_{i+1} \cdot z_{i+3} - x_{i+2} \cdot y_{i+3} \cdot z_{i+1} \\- x_{i+3} \cdot y_{i+1} \cdot z_{i+2} + x_{i+3} \cdot y_{i+2} \cdot z_{i+1}) / \det\\ b_{i} = (y_{i+1} z_{i+2} - y_{i+1} z_{i+3} - y_{i+2} z_{i+1} + y_{i+2} z_{i+3} + y_{i+3} z_{i+1} - y_{i+3} z_{i+2}) / \det\\ c_{i} = (-x_{i+1} z_{i+2} + x_{i+1} z_{i+3} + x_{i+2} z_{i+1} - x_{i+2} z_{i+3} - x_{i+3} z_{i+1} + x_{i+3} z_{i+2}) / \det\\ d_{i} = (x_{i+1} y_{i+2} - x_{i+1} y_{i+3} - x_{i+2} y_{i+1} + x_{i+2} y_{i+3} + x_{i+3} y_{i+1} - x_{i+3} y_{i+2}) / \det \end{cases}\\ &\begin{cases} a_{i+1} = (x_i \cdot y_{i+2} \cdot z_{i+3} - x_i \cdot y_{i+3} \cdot z_{i+2} - x_{i+2} \cdot y_i \cdot z_{i+3} + x_{i+2} \cdot y_{i+3} \cdot z_i \\+ x_{i+3} \cdot y_i \cdot z_{i+2} - x_{i+3} \cdot y_{i+2} \cdot z_i) / \det\\ b_{i+1} = (-y_i z_{i+2} + y_i z_{i+3} + y_{i+2} z_i - y_{i+2} z_{i+3} - y_{i+3} z_i + y_{i+3} z_{i+2}) / \det\\ c_{i+1} = (x_i z_{i+2} - x_i z_{i+3} - x_{i+2} z_i + x_{i+2} z_{i+3} + x_{i+3} z_i - x_{i+3} z_{i+2}) / \det\\ d_{i+1} = (-x_i y_{i+2} + x_i y_{i+3} + x_{i+2} y_i - x_{i+2} y_{i+3} - x_{i+3} y_i + x_{i+3} y_{i+2}) / \det \end{cases}\\ &\begin{cases} a_{i+2} = (-x_i \cdot y_{i+1} \cdot z_{i+3} + x_i \cdot y_{i+3} \cdot z_{i+1} + x_{i+1} \cdot y_i \cdot z_{i+3} - x_{i+1} \cdot y_{i+3} \cdot z_i \\- x_{i+3} \cdot y_i \cdot z_{i+1} + x_{i+3} \cdot y_{i+1} \cdot z_i) / \det\\ b_{i+2} = (y_i z_{i+1} - y_i z_{i+3} - y_{i+1} z_i + y_{i+1} z_{i+3} + y_{i+3} z_i - y_{i+3} z_{i+1}) / \det\\ c_{i+2} = (-x_i z_{i+1} + x_i z_{i+3} + x_{i+1} z_i - x_{i+1} z_{i+3} - x_{i+3} z_i + x_{i+3} z_{i+1}) / \det\\ d_{i+2} = (x_i y_{i+1} - x_i y_{i+3} - x_{i+1} y_i + x_{i+1} y_{i+3} + x_{i+3} y_i - x_{i+3} y_{i+1}) / \det \end{cases}\\ &\begin{cases} a_{i+3} = (x_i \cdot y_{i+1} \cdot z_{i+2} - x_i \cdot y_{i+2} \cdot z_{i+1} - x_{i+1} \cdot y_i \cdot z_{i+2} + x_{i+1} \cdot y_{i+2} \cdot z_i \\+ x_{i+2} \cdot y_i \cdot z_{i+1} - x_{i+2} \cdot y_{i+1} \cdot z_i) / \det\\ b_{i+3} = (-y_i z_{i+1} + y_i z_{i+2} + y_{i+1} z_i - y_{i+1} z_{i+2} - y_{i+2} z_i + y_{i+2} z_{i+1}) / \det\\ c_{i+3} = (x_i z_{i+1} - x_i z_{i+2} - x_{i+1} z_i + x_{i+1} z_{i+2} + x_{i+2} z_i - x_{i+2} z_{i+1}) / \det\\ d_{i+3} = (-x_i y_{i+1} + x_i y_{i+2} + x_{i+1} y_i - x_{i+1} y_{i+2} - x_{i+2} y_i + x_{i+2} y_{i+1}) / \det \end{cases} \end{split}
(5.12)

где определитель матрицы:

det=xiyi+1zi+2xiyi+1zi+3xiyi+2zi+1+xiyi+2zi+3+xiyi+3zi+1xiyi+3zi+2xi+1yizi+2+xi+1yizi+3+xi+1yi+2zixi+1yi+2zi+3xi+1yi+3zi+xi+1yi+3zi+2+xi+2yizi+1xi+2yizi+3xi+2yi+1zi+xi+2yi+1zi+3+xi+2yi+3zixi+2yi+3zi+1xi+3yizi+1+xi+3yizi+2+xi+3yi+1zixi+3yi+1zi+2xi+3yi+2zi+xi+3yi+2zi+1\begin{aligned} \det ={}& x_i \cdot y_{i+1} \cdot z_{i+2} - x_i \cdot y_{i+1} \cdot z_{i+3} - x_i \cdot y_{i+2} \cdot z_{i+1} \\ &+ x_i \cdot y_{i+2} \cdot z_{i+3} + x_i \cdot y_{i+3} \cdot z_{i+1} - x_i \cdot y_{i+3} \cdot z_{i+2} \\ &- x_{i+1} \cdot y_i \cdot z_{i+2} + x_{i+1} \cdot y_i \cdot z_{i+3} + x_{i+1} \cdot y_{i+2} \cdot z_i \\ &- x_{i+1} \cdot y_{i+2} \cdot z_{i+3} - x_{i+1} \cdot y_{i+3} \cdot z_i + x_{i+1} \cdot y_{i+3} \cdot z_{i+2} \\ &+ x_{i+2} \cdot y_i \cdot z_{i+1} - x_{i+2} \cdot y_i \cdot z_{i+3} - x_{i+2} \cdot y_{i+1} \cdot z_i \\ &+ x_{i+2} \cdot y_{i+1} \cdot z_{i+3} + x_{i+2} \cdot y_{i+3} \cdot z_i - x_{i+2} \cdot y_{i+3} \cdot z_{i+1} \\ &- x_{i+3} \cdot y_i \cdot z_{i+1} + x_{i+3} \cdot y_i \cdot z_{i+2} + x_{i+3} \cdot y_{i+1} \cdot z_i \\ &- x_{i+3} \cdot y_{i+1} \cdot z_{i+2} - x_{i+3} \cdot y_{i+2} \cdot z_i + x_{i+3} \cdot y_{i+2} \cdot z_{i+1} \end{aligned}