Вычисление частных производных

[sergeik]

В Декартовой плоскости задано множество точек с координатами $(x_i,y_i)$ (частицы жидкости) и значения функции в них $u(x_i,y_i)$ и $v(x_i,y_i)$ ($(u,v)$ - вектор скорости). Необходимо численно посчитать значения производных $\frac{\partial}{\partial x}$ и $\frac{\partial}{\partial y}$ от функций $(u,v)$. Можно ли решить данную задачу с помощью аппроксимации функции сплайнами? В классических учебниках по сплайнам говорится о прямоугольном разбиении области (равномерном и не равномерном), а как быть в данном случае, когда точки накиданы произвольным образом? Если кто-то знает решение данной задачи, сообщите, пожалуйста.

[Zai]

Вам необходимо провести разделение области на конечные ячейки (треугольники), так чтобы вершины ячеек были только в заданных Вами точках. После этого в каждой треугольной ячейке однозначно определяются градиенты скоростей на основе линейной интерполяции. Градиент скорости в заданных Вами точках определяются как среднее с весом площади всех ячеек, прилегающик к данной точке
$\frac {\partial} {\partial x}= \frac {\sum \limits _{i} \frac {\partial} {\partial x}_i s_i} {\sum \limits _{i} S_i}
$\frac {\partial} {\partial y}= \frac {\sum \limits _{i} \frac {\partial} {\partial y}_i s_i} {\sum \limits _{i} S_i}
Использование сплайновой интерполяции, например квадратичными полиномами, существенно усложнит алгоритм. Вам нужно будет построить интерполяцию в треугольниках, так чтобы сохранялась непрерывность производных на ребрах ячеек.

[sergeik]

Zai писал(а):

Вам необходимо провести разделение области на конечные ячейки (треугольники), так чтобы вершины ячеек были только в заданных Вами точках. После этого в каждой треугольной ячейке однозначно определяются градиенты скоростей на основе линейной интерполяции. Градиент скорости в заданных Вами точках определяются как среднее с весом площади всех ячеек, прилегающик к данной точке
$\frac {\partial} {\partial x}= \frac {\sum \limits _{i} \frac {\partial} {\partial x}_i s_i} {\sum \limits _{i} S_i}
$\frac {\partial} {\partial y}= \frac {\sum \limits _{i} \frac {\partial} {\partial y}_i s_i} {\sum \limits _{i} S_i}
Использование сплайновой интерполяции, например квадратичными полиномами, существенно усложнит алгоритм. Вам нужно будет построить интерполяцию в треугольниках, так чтобы сохранялась непрерывность производных на ребрах ячеек.

Спасибо Вам за предложенный способ вычисления производных. Возникло несколько вопросов:
1) На сколько сильно точность вычисления данным способом будет зависеть от триангуляции области? При малом расстоянии между 2 частицами триангуляция будет плохой (очень острый угол появляется в треугольнике).
2) Под линейной интерполяцией понимается интерполяция как в методе конечных элементов вида $f=ax+by+c$? Тогда производная на элементе будет постоянной. И при вырожденном треугольном элементе
3) Данный способ позволяет считать производные только внутри области? А на границе?
4) Отличается ли по точности данный способ от аналогичного вычисления через теорему о среднем и формулу Грина?
5) Если можно, не подскажите ли ссылки на литературу, в которой бы описывался метод подсчета через сплайны.

[VLarin]

На мой взгляд:
1. Сильно, так как при усреднении в качестве веса используются площади треугольников. Если к узлу прилегают треугольники с разными площадями, то значение градиента будет близко к значениям на треугольниках с большими площадями. Хотя максимальные значения внутри КЭ могут иметь место на треугольниках с малыми площадями. Здесь можно такде попробовать в качестве весовых коэф. использовать расстояние до центра масс треугольников.
2. При линейной интерполяции градиенты будуь постоянны внутри треугольника. Но вся суть предложенного Zai ранее метода в том, что мы переходим от постоянных градиентов внутри каждого треугольника к значениям градиентов в узлах. После этого можно опять-таки возпользоваться линейными полиномами для интерполяции градиентов внутри треугольников.
3. После того, как были найдены значения градиентов, можно воспользоваться линейной интерполяцией в пределах каждого треугольника, и таким образом можно вычислить значения градиентов в любых точках, в том числе и на границах.
4. Использование усреднения приводит к ошибкам, ИМХО применение т. Грина будет точнее.
5. Посмотрите литературу по определению градиентов в МКЭ -

Наиболее полные обзоры существующих методов можно найти в:
1. Omeragic D., Silvester P.P. Progress in differentiation of approximate data // IEEE Antennas and Propagation Magazine. – 1996. – Vol. 38, N 1, February.
2. Omeragic D., Silvester P.P. Numerical differentiation in magnetic field postprocessing // Int. J. Numer. Model. – 1996, N 9, P. 99 – 113.

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

[sergeik]

VLarin писал(а):

На мой взгляд:
1. Сильно, так как при усреднении в качестве веса используются площади треугольников. Если к узлу прилегают треугольники с разными площадями, то значение градиента будет близко к значениям на треугольниках с большими площадями. Хотя максимальные значения внутри КЭ могут иметь место на треугольниках с малыми площадями. Здесь можно такде попробовать в качестве весовых коэф. использовать расстояние до центра масс треугольников.
2. При линейной интерполяции градиенты будуь постоянны внутри треугольника. Но вся суть предложенного Zai ранее метода в том, что мы переходим от постоянных градиентов внутри каждого треугольника к значениям градиентов в узлах. После этого можно опять-таки возпользоваться линейными полиномами для интерполяции градиентов внутри треугольников.
3. После того, как были найдены значения градиентов, можно воспользоваться линейной интерполяцией в пределах каждого треугольника, и таким образом можно вычислить значения градиентов в любых точках, в том числе и на границах.
4. Использование усреднения приводит к ошибкам, ИМХО применение т. Грина будет точнее.
5. Посмотрите литературу по определению градиентов в МКЭ -

Наиболее полные обзоры существующих методов можно найти в:
1. Omeragic D., Silvester P.P. Progress in differentiation of approximate data // IEEE Antennas and Propagation Magazine. – 1996. – Vol. 38, N 1, February.
2. Omeragic D., Silvester P.P. Numerical differentiation in magnetic field postprocessing // Int. J. Numer. Model. – 1996, N 9, P. 99 – 113.

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

Спасибо за ответ.
В том то и дело, что в области есть треугольная сетка (хотя местами очень плохая из-за того, что точки где-то сгустились), также есть диаграмма Вороного. Но вот что-то нет мыслей, как точно посчитать значение производных. Эти производные нужны для пересчета конвективного члена в уравнениях движения жидкости $u \frac{\partial u}{\partial x}$, $u \frac{\partial u}{\partial y}$, $v \frac{\partial v}{\partial x}$, $v \frac{\partial v}{\partial y}$. Причем считать его желательно как можно точнее. В связи с этим и возникла мысль об аппроксимации сплайнами, но не нашел способа, как аппроксимировать на плоскости по произвольно заданному набору точек.
А с границей все равно будут трудности, например, в угловой точке области (типа прямоугольника - просто для нее практически не будет соседей).
К сожаление, статей, указанных авторов, в открытой библиотеке не нашел. Попробую еще по МКЭ посмотреть.

[Zai]

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

[sergeik]

Zai писал(а):

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

В том то и дело, что с треугольными элементами не очень хорошо получается. Мы решаем нестационарные задачи лагранжевым условно бессеточным методом (сетка (диаграмма Вороного) нужна только для поиска соседей для построения функций формы), а в ходе решения без особых трудностей можем получить триангуляцию (не очень хорошую, т.к. следим за движением расчетных узлов - в следствии чего и не можем применить конечные разности). Строить какую-нибудь другую сетку не хочется, чужую программу брать тоже. Вот и хотелось бы на том что есть посчитать производные. Встречали работы, в которых вычисление производных производится через аппроксимацию сплайнами, но в этих работах только упоминается об этом, а как делать - вопрос. Вот и заострили свое внимание на сплайнах, может кто-нибудь что-нибудь делал. А на треугольниках, пожалуй, кроме Вами предложенного способа и через теоремы Грина и о среднем, больше ничего и нет. Классические книги по сплайнам (Завьялов, Квасов, Бур) применяют везде прямоугольную сетку, а нам нужно по произвольному набору точек, который получается после каждого шага по времени.

[Zai]

Попробуйте квадратичную интерполяцию на треугольнике с промежуточными узлами. Узлов у Вас много и это вполне возможно.
$f(x,y)=a_{11}x^2+2a_{12}xy+a_{22}y^2+b_1x+b_2y+c
Шесть коэффициентов в квадратичной форме однозначно определяются по шести значениям скоростей. Такая интерполяция уже несколько десятилетий используется в МКЭ. В регулярной зоне она очень эффективна. На границе иногда бывает потеря аппроксимации. Вместе с тем Вам все равно необходимо усреднение по ячейкам. Данная квадратичная интерполяция дает непрерывность производной на ребре, насколько точно это я не знаю.

[sergeik]

Zai писал(а):

Попробуйте квадратичную интерполяцию на треугольнике с промежуточными узлами. Узлов у Вас много и это вполне возможно.
$f(x,y)=a_{11}x^2+2a_{12}xy+a_{22}y^2+b_1x+b_2y+c
Шесть коэффициентов в квадратичной форме однозначно определяются по шести значениям скоростей. Такая интерполяция уже несколько десятилетий используется в МКЭ. В регулярной зоне она очень эффективна. На границе иногда бывает потеря аппроксимации. Вместе с тем Вам все равно необходимо усреднение по ячейкам. Данная квадратичная интерполяция дает непрерывность производной на ребре, насколько точно это я не знаю.

Спасибо за совет. Обязательно попробуем запрограммировать и сравнить с 2 предыдущими способами. Только вопрос возникает. Для аппроксимации нужно брать дополнительно по одной точке в центре ребра треугольника. Значение скорости в этой точке брать как среднее арифметическое по скорости на концах ребра? А затем действовать по алгоритму предложенному Вами ранее.

[Zai]

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

[sergeik]

И еще один вопрос возник. Как посчитать производную в узлах некоторой кривой (она может быть и прямой), если в них задана функция? По сути, задача таже самая, только вдоль некоторой линии. Есть алгоритм автора Longuet-Higgiens 5 порядка точности, но он работает для $\frac{\partial u}{\partial s}$. Когда граница горизонтальная или вертикальная прямая, то отлично работает. Не встречался ли кто-нибудь с подобной задачей?

[VLarin]

Цитата:

В том то и дело, что в области есть треугольная сетка (хотя местами очень плохая из-за того, что точки где-то сгустились), также есть диаграмма Вороного. Но вот что-то нет мыслей, как точно посчитать значение производных. Эти производные нужны для пересчета конвективного члена в уравнениях движения жидкости , , , . Причем считать его желательно как можно точнее. В связи с этим и возникла мысль об аппроксимации сплайнами, но не нашел способа, как аппроксимировать на плоскости по произвольно заданному набору точек.
А с границей все равно будут трудности, например, в угловой точке области (типа прямоугольника - просто для нее практически не будет соседей).
К сожаление, статей, указанных авторов, в открытой библиотеке не нашел. Попробую еще по МКЭ посмотреть.

Ну тогда можно еще оригинал одного из метода посмотреть:
1. Zienkiewicz O.C., Zhu J.Z. The superconvergence patch recovery and a posteriori error estimates. Part 1: The recovery technique // Int. j. numer. methods eng. – 1992. – N 33, P. 1331–1364.
2. Zienkiewicz O.C., Zhu J.Z. The superconvergence patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity // Int. j. numer. methods eng. – 1992. – N 33, P. 1365–1382.

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

[sergeik]

VLarin писал(а):

Цитата:

Ну тогда можно еще оригинал одного из метода посмотреть:
1. Zienkiewicz O.C., Zhu J.Z. The superconvergence patch recovery and a posteriori error estimates. Part 1: The recovery technique // Int. j. numer. methods eng. – 1992. – N 33, P. 1331–1364.
2. Zienkiewicz O.C., Zhu J.Z. The superconvergence patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity // Int. j. numer. methods eng. – 1992. – N 33, P. 1365–1382.

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

Спасибо за совет. К сожалению, не могу получить доступ к статьям Зенкевича (этот журнал закрыт с нашего логина на elibrary), а вот статью его соавтора Zhu J.Z. нашел "Superconvergence of The Derivative Patch Recovery Technique and A Posteriori Error Estimation". Поразбираюсь, возможно, это то, что нужно. Извините за нескромность, а Вы не могли бы выслать эти статьи по email?

[VLarin]

sergeik писал(а):

Извините за нескромность, а Вы не могли бы выслать эти статьи по email?

Нет:(
Есть только Сильвестер, в виде ксерокопий и то на работе, а окажусь я там еще не скоро:(