Об одном подходе к интерполяции значений неизвестных в задачах улучшения качества расчетной сетки Текст научной статьи по специальности «Математика»

МАТЕМАТИКА

УДК 532.5: 519.652

ОБ ОДНОМ ПОДХОДЕ К ИНТЕРПОЛЯЦИИ ЗНАЧЕНИЙ НЕИЗВЕСТНЫХ В ЗАДАЧАХ УЛУЧШЕНИЯ КАЧЕСТВА РАСЧЕТНОЙ СЕТКИ

Т. С. Рейн, С. Н. Карабцев

AN APPROACH TO INTERPOLATION OF UNKNOWN FUNCTIONS IN THE PROBLEMS OF CALCULATION MESH QUALITY IMPROVEMENT

T. S. Reyn, S. N. Karabtsev

Работа выполнена в рамках задания № 2014/64 на выполнение государственной работы «Организация проведения научных исследований».

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

The paper presents an approach to solving the task of calculating the value of some function in a set point basing on its values set on some fixed system of points. The approach is based on the use of interpolation functions of Sibson’s and Laplace’s shape-functions. The properties of both functions are considered, the interpolation accuracy is defined by solving a Poisson equation with the method of trial functions.

Ключевые слова: интерполяция Сибсона, интерполяция Лапласа, функции формы, бессеточные методы, численное моделирование, комплексы программ, перестроение расчетной сетки.

Keywords: Sibsonian interpolation, Laplace interpolation, shape function, mesh-free methods, numerical simulation, mesh refinement.

Введение

Для моделирования неустановившихся физических явлений, описываемых различными уравнениями в частных производных, создан целый класс сеточных методов, для которых на каждом временном шаге требуется знание информации о межузловой связности. Наиболее популярные представители этой группы - метод конечных элементов (МКЭ) [5] и метод контрольных объемов (МКО) [6]. Основным недостатком МКЭ и МКО является то, что для каждого временного шага сетка, на которой строится решение, не теряет свою узловую связность, что, в свою очередь, при больших деформациях жидкости может быстро приводить к вырожденности сетки. Для корректного разбиения области течения на треугольные элементы могут использоваться различные способы. При этом точность интерполяции зависит от качества триангуляции; например, в работе [8] показано, что для МКЭ якобиан преобразования треугольных и четырехугольных элементов зависит от минимального угла элемента. Дополнительно, неединственность триангуляции приводит к структуре соседей, которая может резко меняться при малом изменении координат узлов и, тем самым, приводить к изменениям результата интерполяции.

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

В настоящее время развитие получили численные методы, которые аппроксимируют уравнения в част-

ных производных, основываясь только на наборе узлов, без знания дополнительной информации о структуре сетки. Данные методы известны как методы частиц или бессеточные методы (БМ) [13]. К группе методов, зарекомендовавших себя в задачах гидродинамики, относятся метод сглаженных частиц (Smoothed Particle Hydrodynamics) [14], полунеявный метод движущихся частиц (Moving Particle Semi-implicit) [12], метод Лагранжево-Эйлеровых частиц [9]. К недостаткам этих бессеточных методов можно отнести сравнительно невысокую точность и трудность введения граничных условий. Дополнительно, в стандартных бессеточных методах для построения интерполяционных функций требуется обеспечить узловую связность.

Эти обстоятельства заставили исследователей искать новые, условно-бессеточные методы, сочетающие в себе идеи и возможности бессеточного подхода, но, вместе с тем, обладающие достоинствами сеточных методов. Первыми из бессеточных методов нового поколения появились бессеточный метод конечных элементов (MFEM) и метод естественных соседей (Natural Element Method) [10; 15].

Особенность данных методов в том, что для стационарных задач они является обычными (классиче-скимии) методами Галеркина, то есть являются сеточными. Для нестационарных задач применяется подход Лагранжа к описанию среды. А именно, на каждом шаге по времени по найденному на предыдущем шаге положению узлов строится сетка, определяющая новую структуру соседей для каждой узловой точки области; на вновь построенной сетке снова решается методом Галеркина система уравнений. В силу этого методы NEM и MFEM сохраняют некоторые преимущества классического метода Галеркина, а

Т. С. Рейн, С. Н. Карабцев

57

МАТЕМАТИКА

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

Основное отличие классического метода конечных элементов от бессеточного, помимо интерполяционных коэффициентов для приближения значений неизвестных, - необходимость перестройки расчетной сетки на каждом временном шаге. Каждый раз на вновь построенной сетке строятся новые интерполяционные функции Сибсона и Лапласа [17], вычисляются матрицы жесткости, строится результирующая система уравнений, и решение многомерной задачи на дробных шагах сводится к решению отдельных уравнений методом сопряженных градиентов [1; 2].

Для определения структуры соседей и соседних узлов [7], а также для вычисления функций формы Сибсона и Лапласа, используется разбиение области ячейками Вороного [4]. Среди преимуществ такого разбиения можно отметить его единственность, что, в свою очередь, гарантирует единственность значений коэффициентов и интерполяции. Сама расчетная сетка необходима для численного интегрирования при вычислении элементов матрицы жесткости результирующей системы уравнений (рассматривается слабая форма уравнений движения среды).

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

1. Интерполяционные функции Сибсона

Алгоритмы решения задачи о вычислении в заданной точке x0 = (x0, y0 ) значения некоторой

функции f (x) по ее значениям f }• заданным на фиксированной системе точек \xk }, сводятся тем или иным способом к формуле вида:

f0 = T;aif , ECi = 1 C > ° i = 1К. (1)

i=1 i=1

Здесь CCi - коэффициенты интерполяции, зависящие от расположения системы точек и не зависящие от функции f , K - количество точек, по которым производится интерполяция функции f . Дополни-

тельным требованием на коэффициенты C i является требование точного выполнения равенства (1) для случая линейной функции f (интерполяция первого

K

порядка). При этом условие Е Ci = 1 обеспечивает

i =1

точное выполнение равенства (1) для функции f ,

к

равной константе. Условия Ci > 0 и Ес = 1

i=1

обеспечивают ограниченность нормы результата интерполяции:

||fo|| ^ Ес (maif 11)=maX f I Ec=maX f\.

i=1 1 1 i=1 1

Здесь и далее рассматривается пространство En измеримых на отрезке [a, b\ функций с интегрируемыми по Лебегу квадратами.

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

1.1. Ячейки Вороного первого и второго порядка

Для заданной точки Pi, i = 1 ..n, n - число узлов

области, многоугольником (ячейкой) Вороного называется геометрическое место точек на плоскости, которые находятся к Pi ближе, чем к любой другой заданной точке P , i Ф j.

J 7 J

Диаграммой Вороного заданного множества точек {Pj,.. Pn } называется совокупность всех многоугольников Вороного заданной системы точек. Строгое определение многоугольника Вороного T вводится следующим образом:

T = {* = (x,y) е R2 : d(x,xt) < d(x,x.), Vi Ф j}.

Диаграмма Вороного второго порядка множества

узлов N есть подразбиение плоскости на ячейки T..,

j

где каждая область Tij ассоциируется с узловой парой соседей (ni , n j ) , такой, что Tij является геометриче-

ским местом всех точек, которые имеют nt в качестве ближайшего соседа, и n j в качестве второго ближайшего соседа [11]. Подчеркнем, что ячейка является непустой тогда и только тогда, когда ni и n j -

соседи.

Строгое определение ячейки Вороного второго порядка для точки x = (x, y) можно представить в следующем виде:

x = (x, у) е R2:

Tj =\: d (x, x ) < d (x, x.) < d (x, xk ), >.

Vk Ф i, j

58 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

МАТЕМАТИКА

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

естественных соседей точки x и их координаты. Понятие соседних узлов является расширением и обобщением определения естественных соседей.

Рис. 1. а) диаграмма Вороного с введенной в нее точкой;

б) к построению интерполяции Сибсона;

в) к построению интерполяции Лапласа

На рис. 1а показана точка x, расположенная внутри диаграммы Вороного множества N . Если x рассматривать как узел вместе с множеством узлов N, то естественные соседи для x - это те узлы, которые в новом разбиении формируют треугольник, содержащий точку x . Видно, что x имеет четыре естественных соседа, а именно, узлы 1, 2, 3, и 4. Можно отметить, что точка x принадлежит ячейке Вороного первого порядка для узла с номером 3 (рис. 1а).

1.2. Определение функций формы Сибсона Функции формы Сибсона, основанные на понятии естественных соседей, базируются на диаграммах Вороного первого и второго порядков и определяются через отношение площадей многоугольников в двумерном случае [16]:

Рг (x) = Ax)., A(x) = J Aj (xX (2)

A(x)

где A(x) - площадь ячейки Вороного первого порядка, в которую входит точка x, Ai (x) - площадь пересечения ячейки Вороного второго порядка точки x и xj (рис. 1,б). Если точка x ^ xi, тогда

Р, (x) = 1 и все другие функции формы равны нулю.

Производные функций формы Сибсона можно получить дифференцированием уравнения (2):

Ai, j (x) ~Рг (x) A,; (x)

-w =---------A(x)-------=

Рг,, (x) =

j = 1,2.

Строгое определение интерполяции Сибсона можно сформулировать в следующем виде. Пусть

{xi}, i = 0,..,N - точки евклидова пространства En

и

T = j x = (x, У) e En : |

k d(x, xk) < d(x, xm), Vk Ф m\

x = (x , y) e En :

Tkm =j: d (x, xk ) < d (x, xm ) < d (x, xi ), >

Vl Ф m, k

4.

где d - евклидово расстояние, Tk - ячейка Вороного для xk . Тогда Tk = U m*kTkm . Если |Tk| < , то

J\Tkm\xm = \Tk\xk , где И - мера Лебега в E" . Оп-

m^k

ределенные таким образом |Tkm | называются коэффициентами интерполяции Сибсона. При этом pm вво-

T I/ _

дится как pm = 1 kmV I (запись для точки xk) [3].

/ \Tk

Ниже приводятся основные свойства сибсонов-ской интерполяции:

- для любого i Ф j pt (x j ) = Sj - символ Кро-некера;

n

- для i = 1..П, p(x) > 0, Jp(x) = 1;

i=1

- x = Jp (x)xi.

i

Доказательства приведенных выше свойств интерполяционных функций Сибсона приводятся в работе [15].

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

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 59

МАТЕМАТИКА

кроме узловых значений. В узлах - непрерывность

Г’О

класса C .

1.3. Вычисление функций формы Сибсона

Для вычисления функций формы Сибсона необходимо искать пересечение ячеек Вороного первого и второго порядков, что является очень трудоемкой задачей вычислительной геометрии. Д. Уотсоном [18] был предложен алгоритм поиска площадей пересечения многоугольников путем разбиения многоугольника ориентированными треугольниками (в каждом из треугольников задается порядок обхода вершин).

2. Интерполяция Лапласа

Построение интерполяции Лапласа также опирается на определение соседей посредством разбиения области ячейками Вороного.

(4) интерполяция Лапласа точно воспроизводит функции равные константе.

Линейная полнота. Несибсоновская интерполяция удовлетворяет следующему свойству для локальных

координат: X = (X)Xi, где xt - естественные

i

соседи точки x .

В соответствии с приведенными свойствами интерполяция Лапласа может точно аппроксимировать линейную функцию:

f(Х) = Т • Х + t = Т • ) +1 •(^Ф ) =

i i

= ХФ (ТХг + 1) = TjVf (Х )'

2.1. Определение функций формы Лапласа Пусть точка X = (х, у) принадлежит многоугольнику Вороного с числом сторон равным N. Обозначим длины сторон многоугольника через

s, i = 1,.., N, а высоты, опущенные из точки

X = (x, у) на si, - через kt. Тогда (1) примет следующий вид:

f = ^f , Ф =

i=1

=(si 1k) Z s-1 kj

V j =

,i = 1..N.

(3)

Такой способ определения коэффициентов (pt

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

Обозначим (X) = st (X) I kt (X) . Тогда произ-

водные функций формы Лапласа могут быть получены дифференцированием уравнения (3):

a, j(Х) -Фг(x)a,}

фi, j = N

Za(Х)

j=1

(Х)

j = 1,2.

2.2. Свойства интерполяции Лапласа

Ниже приведем основные свойства интерполяционных функций Лапласа в соответствии с общими требованиями к интерполирующим функциям, представленным в выражении (1).

Разложение единицы. Интерполяционные функции формы Лапласа представляют разложение единицы:

П

Z^' (Х) = 1. (4)

i=1

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

2.3. Вычисление функций формы Лапласа Алгоритм вычисления функций формы Лапласа для случая, когда ячейка Вороного имеет произвольное число узлов, то есть когда точка X = (X, у) имеет

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

am (Х) = sm (Х) 1 km (Х) = \Гш (Х) - lm (Х)[

где

r (Х) = (Хт - Xm+1 )(Хт+1 - х) + (Ут - Ут+1 )(Ут+1 - У) , (5) т (Хт - Х)(Ут+1 - У) - (Хт+1 - х)(Ут - У)

t (Х) = (Хт - Хт-1 )(Хт-1 - х) + (Ут - Ут-!)(Ут-1 - У) . (6) т (Хт - х)(Ут-1 - У) - (Хт-1 - х)(Ут - У)

Используя соотношения (5) и (6), можно построить функции формы Лапласа.

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

3. Носитель функций формы Сибсона и Лапласа. Расширенная интерполяция Лапласа

Различные алгоритмы интерполяции отличаются

способами выбора коэффициентов (р{ и методами

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

60 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

МАТЕМАТИКА

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

а)

Рис. 2. а) функция формы Сибсона; б) функц

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

б)

формы расширенной интерполяции Лапласа

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

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

Рис. 3. а) носитель функции формы Сибсона; б) носитель функции формы Лапласа

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 61

МАТЕМАТИКА

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

и соседями точки Х0, введенной в первоначальное

разбиение, являются вершины EDT, которой Х0 принадлежит. Носителем функции формы Лапласа является многоугольник расширенной триангуляции Делоне, содержащий точку (рис. 3,б).

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

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

4. Сравнение интерполяций Сибсона и Лапласа на решении уравнения Пуассона

4.1. Постановка задачи

Точность использования интерполяционных функций Сибсона и Лапласа в плоском случае проверяется на решении уравнения Пуассона в области D при заданном значении правой части (метод пробных функций): Au = b , с граничными условиями:

Дирихле: u = и на участке Ги границы Г; du ~

Неймана: q = дп = q на участке Гч границы

Г, где п = (nx, пу) - внешняя единичная нормаль к поверхности; и, q - заданные значения функции и ее нормальной производной на границе Г = Ги + ГЧ.

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

4.2. Краткое описание метода Галеркина

Метод Галеркина является частным случаем метода взвешенных невязок, в котором весовые функции совпадают с базисными.

Пусть дана система уравнений:

L(u) - p = 0, Х gQ , F(u) - q = 0, Х e Г

N

и аппроксимирующая функция и = ^ Okuk ,

k=1

удовлетворяющая граничным условиям

F(u) = q(x), Х g Г.

Здесь uk - узловые значения искомой функции, ОСк - функции формы. Невязка

N

s = L(E°kuk) - Р

k=1

ортогонализируется по отношению к базисным функ-

циям (Oi:

Отсюда:

(е,а) = 0, i = 1,..., N.

(7)

Л

{[L(Z°kuk) -p oidQ=ai = l,..,n .

□ V k=1 У

Если L - линейный оператор, то система (7) переходит в систему линейных алгебраических уравне-

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

4.3. Численные результаты

Тесты проводились в областях с различной геометрией и с различным числом расчетных узлов. Первая группа тестов была проведена на линейных функциях. При этом оказалось, что для таких функций максимальная относительная погрешность интерполяции Сибсона (метод NEM) составила 10 6, а интерполяции Лапласа - 4 -10 5 (метод MFEM). Значение относительной погрешности As вычислялось следующим образом:

max u - u

As =------'-r-r-1. (8)

max| u

Здесь u - аналитическое решение уравнения Пуассона, записанное согласно методу пробных функций, u - найденное численное решение.

Далее приводятся результаты для функции, сильно отличающейся своими значениями в узлах сетки. Рассмотрим область D, верхняя граница Г1 которой задается уравнением y = 0,5sin(х), х e [0,2*], а

боковые и нижняя границы Г2 , Г3, Г4 являются прямыми линиями, а именно: Г2 : Х = 0, y e [- 1,0]; Г3 : у = -1, x e [0,2*]; Г4 : x = 2*, y e [-1,0].

В описанной выше области решалось уравнение

2

Пуассона (8) для функции u = Ху с граничными

условиями Дирихле на Г1 и Неймана на границах

Г Г Г :

± 2’ А 3’ А 4 •

du

u = xy2, (х,у) e Г!; — = -у2,(x,у) еГ2,

дп

du _ , . „ du 2 , . „

— = -2xy,(x,у) еГз; — = у ,(x,у) e Г4.

dn дп

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

62 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

МАТЕМАТИКА

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

рования nint, а именно: nint = 7, nint = 13,

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

Рис. 4. а) погрешность решения (функции формы Сибсона); б) время расчета (функции формы Сибсона)

На рис. 4а приведены кривые погрешности решения уравнения Пуассона в описанной области методом пробных функций. В качестве интерполяционных функций выбирались функции формы Сибсона (метод NEM). Кривые под номерами 1, 2, 3 соответствуют семи, тринадцати и двадцатипетиточечной схемам интегрирования соответственно. Погрешность As вычислялась по формуле (8).

На рис. 4б представлены графики зависимости времени расчета (секунды) от числа узлов области также для nint = 7, nint = 13, nint = 25.

Из представленных графиков видно, что оптимальным в смысле выбранной точности и временных затрат являются результаты расчетов, соответствующие nint = 13. При nint = 25 показаны более точные

результаты, однако в этом случае существенно увеличивается время расчета для числа узлов области большего 4500.

В таблицах 1 и 2 для тех же значений точек интегрирования приводится относительная погрешность значений функции и время расчета для различного числа узлов расчетной области при использовании интерполяции метода NEM.

Относительная погрешность интерполяции Сибсона A(u)

Таблица 1

1209 3122 5552 6972 8668 10500

7 1,7942 • 10-3 8,8874•Ш 4 7,1326 • 10-4 5,9086 40-4 5,3996 •Ш-4 4,355840 4

13 1,1387 •Ш 3 8,0157 • 10 4 6,4594 -10 4 5,9679 40-4 5,1295 40-4 3,7801 •Ю-4

25 1,0760 •Ш 3 4,4341 • 10-4 2,838140-4 2,460040 4 2,1035 •Ш-4 1,7780 40-4

Таблица 2

Время решения уравнения Пуассона с использованием интерполяции Сибсона

1209 3122 5552 6972 8668 10500

7 0,810 0,343 0,875 1,297 1,875 2,687

13 0,125 0,468 1,126 1,625 2,297 3,235

25 0,203 0,719 1,640 2,313 3,218 4,406

Ниже на рис. 5 и рис. 6 представлены результаты аналогичных расчетов, полученные при использовании в качестве интерполяционной функции функции

формы Лапласа. А именно: на рис. 5 приведены погрешность решения и время расчета для функции формы Лапласа, носитель которой совпадает с носи-

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 63

МАТЕМАТИКА

телем функции формы Сибсона (расширенная интерполяция Лапласа). На рис. 6 интерполяционная функция строилась на многоугольниках расширенной триангуляции Делоне (метод MFEM).

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

новскими функциями формы. При этом для интерполяции бессеточного метода конечных элементов (рис. 6) наблюдается зависимость погрешности решения от числа узлов, а именно: график погрешности не является монотонно убывающим, в диапазоне узлов от 7500 до 9500 виден рост ошибки интерполяции.

Рис. 5. а) погрешность решения (расширенная интерполяция Лапласа); б) время расчета (расширенная интерполяция Лапласа)

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

му для дальнейших расчетов было выбрано число точек интегрирования wint = 13.

В таблицах 3 - 6 приведены данные расчетов (погрешность интерполяции и время расчетов) для не-сибсоновской интерполяции и расширенной интерполяции Лапласа.

Рис. 6. Функции формы Лапласа а) погрешность решения; б) время расчета

Относительная погрешность расширенной интерполяции Лапласа А(и)

Таблица 3

1209 3122 5552 6972 8668 10500

7 1,7941-10 3 9,2896 -10 4 6,9703-10 4 5,8920 -10 4 5,2443 -10 4 4,3567 -10 4

13 1,1387 -10 3 6,0531-10 4 5,0360-10 4 5,3779-10~4 5,2752 -10~4 4,6223 -104

25 1,0760 -10 3 4,4341-10 4 2,7591-10 4 2,3885 -10 4 2,1057 -104 2,0521-10 4

64 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1

МАТЕМАТИКА

Таблица 4

Время решения уравнения Пуассона с использованием расширенной интерполяции Лапласа t [c]

1209 3122 5552 6972 8668 10500

7 0,094 0,375 0,985 1,4530 2,140 2,485

13 0,125 0,531 1,359 1,953 2,797 3,235

25 0,219 0,843 2,078 2,906 4,047 4,765

Относительная погрешность интерполяции Лапласа А(и)

Таблица 5

1209 3122 5552 6972 8668 10500

7 9,8544-10 4 7,2259 -10~4 1,5713-10 3 2,0601-10 3 2,1315-10 3 1,8822 -10 3

13 1,0495 -10 3 5,1352 -104 4,9359 -10 4 4,9666 -10 4 6,6723 -10 4 7,7108 -10 4

25 9,7954 -10 4 4,9590 -10 4 4,3861-10 4 3,8822-10~4 8,0551-10 4 1,0665 -10 3

Таблица 6

Время решения уравнения Пуассона с использованием интерполяции Лапласа

1209 3122 5552 6972 8668 10500

7 0,063 0,312 0,844 1,266 1,891 5,020

13 0,930 0,406 1,078 1,594 2,343 3,422

25 0,140 0,609 1,562 2,281 3,313 7,001

Для анализа приведем сравнение относительной погрешности (рис. 7,а) и времени расчета (рис. 7б) трех видов интерполяции для числа узлов интегрирования Hint = 13. Кривая с номером 1 соответствует

сибсоновским функциям формы, с номером 2 - функциям формы расширенной интерполяции Лапласа, с номером 3 - несибсоновским функциям формы (функциям формы Лапласа).

Рис. 7. Сравнение алгоритмов интерполяции. а) погрешность решения; б) время расчета

Из приведенного на рис. 7 сравнения различных алгоритмов интерполяций можно сделать вывод, что при одинаковом числе узлов области несибсоновская интерполяция требует меньше времени, чем интерполяция естественных соседей и расширенная интерполяция Лапласа. Это объясняется тем, что методы построения функций формы Сибсона включают в себя алгоритмы поиска площадей пересечения много-

угольников, что является трудоемкой задачей вычислительной геометрии.

Заключение

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

Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1 65

МАТЕМАТИКА

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

формы Лапласа. А так как размер результирующей матрицы системы уравнений ограничен размером оперативной памяти ЭВМ, то, несмотря на некоторые преимущества во времени построения несибсонов-ских функций формы, выгоднее использовать для интерполяции неизвестных при перестройке расчетной сетки функции формы Сибсона.

Литература

1. Афанасьев К. Е., Карабцев С. Н., Рейн Т. С., Стуколов С. В. Численное моделирование течений жидкости со свободными границами бессеточными методами // Труды X Всероссийского съезда по фундаментальным проблемам теоретической и прикладной механики. Нижний Новгород, 2011.

2. Афанасьев К. Е., Рейн Т. С. Моделирование задач гидродинамики вязкой несжимаемой жидкости со свободными границами бессеточным методом естественных соседей // «Вычислительные технологии». Т. 13. № 4.

2008. С. 7 - 24.

3. Беликов В. В., Иванов В. Д., Конторович В. К., Корытник С. А., Семенов А. Ю. Несибсоновская интерполяция - новый метод интерполяции значений функции на произвольной системе точек // Вычислительная математика и математическая физика. 1997. Т. 37. № 1. С. 11 - 17.

4. Карабцев С. Н., Стуколов С. В. Эффективный алгоритм генерации конечно-элементной сетки для метода естественных соседей // Материалы III Международной научной летней школы «Гидродинамика больших скоростей и численное моделирование». Кемерово: ИНТ, 2006. С. 401 - 409.

5. Коннор Дж., Бреббия К. Метод конечных элементов в механике жидкости. Л.: Судостроение, 1979. 264 с.

6. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат,

1984. 152 с.

7. Скворцов А. В. Триангуляция Делоне и ее. Томск: ТГУ, 2002. 128 с.

8. Терентьев А. Г., Афанасьев К. Е. Численные методы в гидродинамике: учебное пособие. Чебоксары: Чуваш. ун-т., 1987. 80 с.

9. Франк А. М. Дискретные модели несжимаемой жидкости. М.: Физматлит, 2001. 208 с.

10. Facundo P. The meshless finite element method applied to a lagrangian particle formulation of fluid flows // Partial Fulfillment of the Requirements for the Degree Doctor of Philosophy. Instituto de Desarrollo tecnologico para la industria quimica (INTEC) universidad nacional del litoral noviembre. 2003. 157 p.

11. Fortune S. J. A sweepline algorithm for Voronoi diagrams // Journal Algorithmica. 1987. № 2. P. 153 - 174.

12. Koshizuka S., Tamako H., Oka Y. A particle method for incompressible viscous low with fluid fragmentation // Computational Fluid Dynamics Journal. 1995. Vol. 4. № 1. Р. 29 - 46.

13. Liu G. R. Mesh free methods: moving beyond the finite element method. London: CRC Press, 2003. 693 p.

14. Monaghan J. J., Thompson M. C., Hourigan K. Simulation of free surface flows with SPH // Journal of computational physics. 1994. Vol. 110. Р. 399 - 406.

15. Onate E., Idelsohn S. R., Zienkiewicz O. C. A finite point method in computational mechanics. Applications to convective transport and fluid flow // International Journal for Numerical Methods in Engineering. 1996. Vol. 39.

P. 3839 - 3866.

16. Sibson R., Barnett In V. A brief description a natural neighbor interpolation // Interpret multivariate data. Chichester: John Wiley, 1981. P. 21 - 36.

17. Sukumar N., Moran B., Belytschko T. The natural element method in solid mechanics // International journal of numerical methods in engineering. 1998. Vol. 43. № 5. Р. 839 - 887.

18. Watson D. F., Philip G. M. Neighborhood-based interpolation // Geobyte. 1987. Vol. 2. № 2. P. 12 - 16.

Информация об авторах:

Рейн Татьяна Сергеевна - кандидат физико-математических наук, доцент кафедры ЮНЕСКО по НИТ, КемГУ, [email protected].

Tatiana S. Reyn - Candidate of Physics and Mathematics, Assistant Professor at the UNESCO Department for New Information Technologies, Kemerovo State University.

Карабцев Сергей Николаевич - кандидат физико-математических наук, доцент кафедры ЮНЕСКО по НИТ КемГУ, [email protected].

Sergey N. Karabtsev - Candidate of Physics and Mathematics, Assistant Professor at the UNESCO Department for New Information Technologies, Kemerovo State University.

Статья поступила в редколлегию 09.10.2014 г.

66 Вестник Кемеровского государственного университета 2014 № 4 (60) Т. 1