Разработка численных методов Для решения задач сейсмо и электроразведки (стр. 2 )

,

,

,

.

И следующие начальные условия:

,

.

В 4 показано, что тестовая система (2.3.1) − (2.3.6) с граничными условиями (2.3.10) − (2.3.13) и с начальными условиями (2.3.14), (2.3.15) имеет следующее аналитическое решение:

,

,

,

,

,

.

Также в приложении Б.4 приведен вывод матриц и векторов (2.2.27) − (2.2.31) и (2.2.38) − (2.2.45). Они имеют следующий вид:

,

, , ,

, , ,

,

, , ,

,

, , ,

,

.

4 Интерполяция на неструктурированных треугольных и тетраэдральных сетках

В ряде задач механики деформируемого тела, к примеру в задачах сейсмической разведки, используется сеточно-характеристический метод [2–4, 15] на неструктурированных треугольных для двумерных задач и тетраэдральных для трехмерных задач сетках. Как правило, при выборе таких типов сеток ограничиваются методами первого порядка аппроксимации, так как вычисляются только значения в вершинах, а для определения функции, более сложной, чем линейная, требуются значения в большем числе точек [15].

При использовании сеточно-характеристического метода возникает потребность определять решение на пересечении характеристики с предшествующим временным слоем исходя из значений в ближайших узлах. Поскольку координаты точек пересечения зависят от выбора координатных осей, от шага интегрирования, точного вида сетки в окрестности точки, то желательно уметь интерполировать [1] их в любой точке. Также процедуры интерполяции могут быть востребованы в случае использования деформируемых сеток, при введении новых треугольников (тетраэдров), значения в вершинах которых необходимо проинициализировать.

Ниже приведены процедуры интерполяции полиномами до пятого порядка включительно на треугольных и тетраэдральных сетках. Если у интерполируемой функции имеются разрывы, интерполяция полиномами степеней выше первой приводит к нефизичным флуктуациям интерполянта. В [15-16] введено понятие полной вариации (total variation) дискретной и континуальной функций соответственно:

,

.

Так как в сеточно-характеристическом методе интегрирование заменяется копированием решения с предыдущего временного слоя на новый слой, полная вариация решения на новом слое будет не больше полной вариации интерполянта предыдущего слоя. Поэтому если метод интерполяции не увеличивает полную вариацию, то полученная схема будет TVD-схемой [15, 16]. Использование порядка интерполяции выше первого не обеспечивает сохранения полной вариации. Поэтому возникает необходимость в ограничителе полинома [15].

Так же другим способом устранения осцилляций является поиск безусловных экстремумов, путем составления системы уравнений, приравнивающей градиент интерполянта к нулю. Но этот способ плох, так как при порядке интерполяции выше второго требует решения нелинейной системы уравнений и не учитывает наличие условных экстремумов. В случае наличия осцилляций в интерполянте предлагается пользоваться в данном треугольнике (тетраэдре) кусочно-линейной интерполяцией по всем имеющимся опорным точкам[3]. Достоинство кусочно-линейной интерполяции в том, что интерполянт не может иметь экстремумов ни в каких точках, кроме опорных.

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

4.1 Интерполяция высоких порядков на треугольных сетках

Будут приведены процедуры полиномиальной, кусочно-линейной и монотонной интерполяций, формулы для вычисления градиента интерполяционного полинома по значениям в опорных точках треугольника, а также процедуры монотонной интерполяции до третьего порядка включительно. Алгоритм построения этих интерполянтов, а также вычисления производных их полиномов изложен в [15].

4.1.1 Интерполяция полиномами на треугольных сетках

Для определения полиномиального поля степени , зависящего от x и y, требуется знать значения в точках, которые в дальнейшем будут называться опорными. Предлагается следующий способ расстановки опорных точек. В треугольнике проводятся прямые, параллельные его сторонам и делящие каждую на равных частей.

Каждой прямой присваивается номер. Прямые, содержащие вершины, имеют номера, равные ; прямые, содержащие стороны, − номера, равные . Каждой опорной точке присваивается номер вида по номерам прямых, на пересечении которых она лежит. Вначале ставится номер прямой, параллельной стороне, противоположной вершине , затем ставится номер прямой, параллельной стороне, противоположной вершине и т. д.

Когда речь идет о треугольных сетках, под векторами будем понимать двумерные векторы вида , когда речь идет о тетраэдральных сетках – трехмерные.

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

Обозначим векторы вершин треугольника , , . Для каждой точки вводятся следующие величины:

,

,

.

Здесь под обозначением понимается следующая операция:

.

Для правильной работы процедуры интерполяции нужно, чтобы выполнялось следующее строгое неравенство:

.

Для каждого по приведенным ниже формулам вычисляются веса опорных точек:

.

Значение полинома в искомой точке определяется по формуле

.

где − значение интерполируемой функции в опорной точке .

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

,

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

Ниже приводятся формулы (3.1.9) − (3.1.63) весов опорных точек для интерполяции полиномами различных степеней: от первой до пятой включительно. А также рисунки 3.1.1 − 3.1.5, поясняющие положение опорных точек в треугольнике. Вывод этих формул дается в 1.

Линейная интерполяция

Рисунок 3.1.1. Опорные точки линейной интерполяции в треугольнике

Веса опорных точек определяются следующими формулами:

,

,

.

Квадратичная интерполяция

Рисунок 3.1.2. Опорные точки квадратичной интерполяции в треугольнике

Веса опорных точек определяются следующими формулами:

,

,

,

,

,

.

Кубическая интерполяция

Рисунок 3.1.3. Опорные точки кубической интерполяции в треугольнике

Веса опорных точек определяются следующими формулами:

,

,

,

,

,

,

,

,

,

.

Интерполяция полиномом четвертого порядка

Рисунок 3.1.4. Опорные точки интерполяции полиномом четвертого порядка в треугольнике

Веса опорных точек определяются следующими формулами:

,

,

,

,

,

,

,

,

,

,

,

,

,

,

.

Интерполяция полиномом пятого порядка

Рисунок 3.1.5. Опорные точки интерполяции полиномом пятого порядка в треугольнике

Веса опорных точек определяются следующими формулами:

,

,

,

,

,

,

,

,

,

,

,

,

,

,

,

,

,

,

,

,

.

4.1.2 Кусочно-линейная интерполяция на треугольных сетках

Для каждого заданно порядка полинома параллельные прямые, определяющие положения опорных точек, разделяют треугольник на подобных ему малых треугольников. Суть процедуры кусочно-линейной интерполяции заключается в следующем. Вначале определяется малый треугольник, в который попадает точка, значение интерполянта в которой надо определить. Потом это значение вычисляется путем линейной интерполяции в малом треугольнике. Веса вершин малых треугольников вычисляются по уже готовым формулам.

Ниже приведены процедуры кусочно-линейной интерполяции в треугольнике по опорным точкам, соответствующим полиномам порядка . Доказательства этих формул даются в 1.

Случай (Опорные точки квадратичной интерполяции).

Рисунок 3.1.6

1) Если , то

.

2) Если , то

.

3) Если , то

.

4) Если , то

.

Случай (Опорные точки кубической интерполяции).

Рисунок 3.1.7

1) Если , то

.

2) Если , то

.

3) Если , то

.

4) Если , то

.

5) Если , то

.

6) Если , то

.

7) Если , то

.

8) Если , то

.

9) Если , то

.

Случай (Опорные точки интерполяции четвертого порядка).

Рисунок 3.1.8

1) Если , то

.

2) Если , то

.

3) Если , то

.

4) Если , то

.

5) Если , то

.

6) Если , то

.

7) Если , то

.

8) Если , то

.

9) Если , то

.

10) Если , то

.

11) Если , то

.

12) Если , то

.

13) Если , то

.

14) Если , то

.

15) Если , то

.

16) Если , то

.

4.1.3 Гибридная интерполяция на треугольных сетках

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

1) Определяем значение пробной функции в этой точке с помощью полиномиальной интерполяции порядка , пусть оно равно .

2) Определяем, в каком из малых треугольников лежит точка . Этот треугольник обозначим за .

3) Сравниваем с и .

3.1) Если , то значение интерполянта в точке : .

3.2) Если , то значение интерполянта в точке : .

3.3) Если , то значение интерполянта в точке : .

Алгоритм построения гибридной интерполяции на основе полиномиальной и кусочно-линейной интерполяции путем поиска безусловных экстремумов пробно полиномиального интерполянта путем решения системы уравнений изложен в

В 1 приведен результат вычисления градиентов полиномов степеней от первой до пятой включительно, а в 3 приведены соответствующие выкладки.

4.1.4 О хранении значений в опорных точках для треугольной сетки

Для топологии опорных точек, использующихся для интерполяции полиномом степени в каждой вершине треугольника должно храниться одно значение интерполируемой функции, на каждой стороне должно храниться значений и внутри каждого треугольника – значений.

При этом передачи интерполятору элементов всех этих массивов, осуществляемой в заранее определенном порядке, достаточно. Т. е. хранить ориентацию самих треугольников и их сторон не нужно.

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

4.2 Интерполяция высоких порядков на тетраэдральных сетках

В этом пункте пойдет речь об интерполяции высоких порядков на неструктурированных тетраэдральных сетках. Будут приведены процедуры полиномиальной, кусочно-линейной и монотонной интерполяций, формулы для вычисления градиента интерполяционного полинома по значениям в опорных точках тетраэдра, а также процедуры монотонной интерполяции до третьего порядка включительно.

4.2.1 Интерполяция полиномами на тетраэдральных сетках

Для определения полиномиального поля степени , зависящего от x, y и z, требуется знать значения в точках, которые в дальнейшем будут называться опорными. Предлагается следующий способ расстановки опорных точек. В тетраэдре проводятся плоскости, параллельные его граням и делящие каждое ребро на равных частей.

Каждой плоскости присваивается номер. Плоскости, содержащие вершины, имеют номера, равные ; плоскости, содержащие грани, − номера, равные . Каждой опорной точке присваивается номер вида по номерам плоскостей, на пересечении которых она лежит. Вначале ставится номер плоскости, параллельной грани, противоположной вершине , затем ставится номер плоскости, параллельной грани, противоположной вершине и т. д.

В каждой вершине тетраэдра находится одна опорная точка, внутри каждого ребра тетраэдра находится опорных точек, внутри каждой грани тетраэдра находится опорных точек и внутри самого тетраэдра находится опорных точек.

Обозначим векторы вершин тетраэдра , , , . Для каждой точки вводятся следующие величины:

,

,

,

.

Здесь под обозначением понимается следующая операция:

.

Для правильной работы процедуры интерполяции нужно, чтобы выполнялось следующее строгое неравенство:

.

Для каждого по приведенным ниже формулам вычисляются веса опорных точек:

.

Значение полинома в искомой точке определяется по формуле

.

где − значение интерполируемой функции в опорной точке .

Как будет видно из приведенных ниже формул для весов опорных точек, они симметричны относительно перестановок, переводящих тетраэдр в тетраэдр . Перестановки должны применяться ко всем индексам по следующему правилу:

,

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

Ниже приводятся формулы (3.2.10) − (3.2.134) весов опорных точек для интерполяции полиномами различных степеней: от первой до пятой включительно. А также рисунки 3.2.1 − 3.2.5, поясняющие положение опорных точек в тетраэдре. На них изображен вид снаружи всех четырех граней тетраэдра. Вывод формул дается в 2.