CC BY
78
Задача изогеометрической сплайн-интерполяции А.Д. Гомонов
НТФ "'Комплексные системы "
Аннотация. В данной статье анализируются перспективы и преимущества применения методов сплайн-интерполяции в прикладной математике. Более подробно рассмотрена метод изогеометрической интерполяции и его отличия от других методов. Дана постановка задачи и пути решения. Приведены практические примеры, демонстрирующие его преимущества перед другими методами.
Abstract. Prospects and advantages of spline-interpolation methods' application in applied mathematics have been analysed in the work. The task of izogeometric interpolation and its differences from other methods has been thoroughly considered. The paper contains the general formulation of the task and methods of its solution. The practical examples demonstrating its advantages have been shown.
1. Введение
Без сплайнов немыслимо создание современных систем автоматизированного проектирования и изготовления изделий сложной геометрической формы, например, поверхностей летательных аппаратов, корпусов судов и т.д. Именно эти задачи послужили первоначальным толчком к началу исследований по сплайнам (Стечкин, Субботин, 1976). Однако практическое применение методов сплайн-интерполяции было осложнено тем, что для быстро осциллирующих функций или для функций с большими градиентами непосредственное применение параболической или кубической интерполяции приводило к нежелательным "выбросам". Для решения данной проблемы был разработан метод изогеометрической интерполяции.
2. Постановка задачи
Помимо вышеуказанных задач, по нашему мнению, сплайн-интерполяцию можно использовать и для моделирования (построения) температурных полей поверхности океана по определенному массиву замеров. Это поможет более точно прогнозировать нахождение косяков рыбы и, как следствие этого, более эффективно использовать ресурсы флота.
Однако, если применять для решения данной задачи классические интерполяционные квадратические или кубические сплайны, то мы можем наблюдать "выбросы", которые в реальности отсутствуют (рис. 1). Это существенный недостаток. Мы считаем, что для его устранения и более точного построения температурных полей целесообразно применять метод изогеометрической интерполяции (рис. 2). Иными словами, необходимо, имея массив замеров, решить задачу восполнения данных.
Задача восполнения данных в технических приложениях решается несколькими способами, выбираемыми в зависимости от целей расчетов и предъявляемых к ним требований.
Наиболее употребительными методами восполнения данных являются методы аппроксимации, которая в общей форме осуществляется следующим образом. Предположим, что аппроксимируемые и подлежащие восполнению данные представляют собой значения некоторой (в общем случае неизвестной) функции.
Для удобства дальнейшего изложения совокупность точек, в которых заданы эти значения, будем называть, как это принято в вычислительной математике, исходной сеткой, а каждую точку из указанной совокупности - узлом исходной сетки. В задачу аппроксимации включается нахождение такой функции, которая достаточно "хорошо" описывала бы данные исходной сетки (Воронин, 1982). Как правило, для аппроксимирующей функции необходимо, чтобы она была "хорошим" приближением. Разумеется, предстоит еще определить, что же представляет собой "хорошее" приближение.
В нашем случае необходимо, чтобы значения аппроксимирующей функции в узлах исходной сетки незначительное отличалась от соответствующих исходных значений и, кроме того, между узлами мало отклонялись от теоретически (или практически) возможных значений исходной функции. Очевидно, что и аппроксимируемая, и аппроксимирующая функции должны быть непрерывными и иметь непрерывные производные (Воронин, 1987).
Критерий, по которому отклонения могут считаться "малыми", определяется многими обстоятельствами: точностью измерения и представления данных, их полнотой и т.д. Часто условием интерполяции является не только хорошее приближение данных, но и их сглаживание.
Одним из наиболее часто применяемых и сравнительно простых методов является метод наименьших квадратов. Согласно ему, сначала, учитывая характер данных, эвристически подбирают сравнительно простую "хорошую" функцию (ее называют эмпирической, опытной) с небольшим числом параметров или коэффициентов. Затем отыскиваются такие значения указанных параметров, которые доставляют минимум
квадрату суммы отклонений значений эмпирической функции в узлах исходной сетки от исходных данных. Полученная таким образом функция позволяет однозначно определить недостающие промежуточные значения. Достоинства метода наименьших квадратов - сравнительная простота и возможность экстраполяции, т.е. определения не только промежуточных значений, но и значений, лежащих за пределами охваченного сеткой интервала. Но у метода есть и недостатки, обусловленные интуитивным выбором типа кривой и малым числом параметров. При большом (относительно параметров) числе узлов, где заданы исходные данные, представляется весьма затруднительным "натянуть" значения эмпирической функции одинаково хорошо на все исходные данные и тем самым обеспечить наиболее желательную точность (Корнейчук, 1984).
Поэтому в практике восполнения данных применяются преимущественно методы, более приспособленные к исходным данным. Это достигается большим, чем число задаваемых значений данных, набором параметров или функций, с помощью которых строится аппроксимация. Одним из подобных, наиболее эффективных методов приближения функций (в том числе интерполяции) является сплайн-аппроксимация. В качестве аппроксимирующей функции, называемой сплайном, строится цепочка полиномов (или других сравнительно простых типов функций). Соседние по цепочке полиномы "склеиваются" так, что их значения и значения их производных на стыках совпадают. При этом все исходные значения попадают на некоторые стыки. Таким образом, внутри промежутков, образованных узлами исходной сетки, оказывается по одному (реже по два или три) полиному из указанной цепочки. Дополнительно определяются значения производных на концах всего интервала. Внутренние стыки, не лежащие в узлах исходной сетки, определяются как дополнительные узлы. Описанные сплайны называются основными. Параметры всех полиномов основного сплайна завязаны друг с другом в систему линейных уравнений. Решая уравнения, находим параметры.
Другой вариант полиномиальной сплайн-аппроксимации - обобщение метода скользящей средней. Каждое значение аппроксимирующей функции (сплайна) складывается из близлежащих исходных значений с коэффициентами в виде некоторых полиномов. Так как одновременно подлежат согласованию не все используемые в этом методе полиномы, а только близлежащие, такого рода сплайны называют локальными или Б-сплайнами.
Признанными достоинствами основных и Б-сплайнов считаются высокая точность, достаточная гладкость и технологическая отработанность. Но проявляются эти достоинства в полной мере лишь на равномерной сетке и при относительно гладких данных. Для восполнения данных и обеспечения их сопоставимости с использованием сплайн-аппроксимации желательно сохранить достоинства сплайнов, но избавиться от вышеуказанных недостатков.
3. Постановка задачи изогеометрической интерполяции
Пус
Обозначим:
Пусть в узлах сетки А: а = х0 < х1 < ... < хп = Ь заданы значения функцииД(х) е С2[а, Ь]:Д =/(х), xi е А.
hi = хм - х„ 8i = (Л+1 - Л) / И, I = 0, ..., п - 1.
Допустим, что в узлах исходной сетки заданы значения первой и второй производной сплайна (например, первые и вторые производные аппроксимируемой функции):
Я = ад, Я"i = ^(х)^ = 0, ..., п.
Обозначим аналог разделенной разности для первой производной сплайна в узлах исходной сетки:
д> = - Я) / И, i = 0, ..., П - 1.
Интерполяцию Дх) и соответствующую аппроксимацию её производных/(х) и/'(х) на интервале [а, Ь] будем осуществлять локальными кубическими сплайнами с двумя дополнительными узлами, разбивающими каждый интервал [х, хт], i = 0,...,п - 1, на три равных промежутка:
ий = х1 + /3, и2 = х1 + 2 /3.
Предполагается, что третья производная может иметь разрыв только в местах склейки (т.е. в узлах исходной сетки и дополнительных узлах). Тогда аппроксимация может быть представлена следующим образом:
ал(х-х,)3 + Ьп(х-х,)2 + сй(х-х,) + йл, х е [х„ ип],
Б(х) = { ай(х-и,1)3 + Ьi2(x-uй)2+ С2(х-ил) + ёа, х е [ий, иа\, (1)
V aa(x-Ui2)ъ + Ьiз(x-ua)2 + еа(х-иа) + х е [ий, хт],
i = 0, 1,...,п - 1.
Предварительно следует рассчитать значения первых и вторых производных сплайна в дополнительных узлах, выразив их через значения разделенных разностей 8г и значения первых и вторых производных сплайна в узлах исходной сетки, так как это позволяет не только быстро вычислять коэффициенты а^, Ьу, вц, d1j, но и анализировать изогеометрическую ситуацию в каждом интервале.
5?(ий) = / 2 - /2 + И,(5", /3 +25",+ / 3) /12,
(иц) = (9$ - 65, - 35т) / И, - 55", / 6 + Я"»! / 3, (2)
5(и,2) = 3 $ / 2 - 5,/2 - И{ (25",■/3 + 5",+ / 3) /12, (2)
^'(и,2) = -(9$ - 35,- - 65-^1) / И, + 5",- / 3 - 55",+ / 6. Прежде всего отметим, что для того, чтобы первые производные сплайна в дополнительных узлах сохраняли знак разделенной разности 8 (как необходимое условие сохранения монотонности самого сплайна в каждом интервале), достаточно, чтобы первые и вторые производные сплайна в узлах исходной сетки сохраняли знак соответствующих разделенных разностей и не более, чем вдвое превосходили их по абсолютной величине:
8 5, > 0, 8 5+ > 0, причем, если 8 = 0, то & / = 5,+ = 0; (3)
8, > 0, 8,5'т > 0, причем, если 8, = 0, то 5", = 5',+ = 0; (4)
|5", | < 2 | 8 |, | ^ < 2 | 8\, 8,Ф 0; (5)
|5', | < 2 | 8, |, | 5'т| < 2 | 8,1 8, Ф 0. (6)
Нетрудно убедиться, что когда эти условия для какого-либо интервала выполняются, то знак первой производной сохраняется и внутри крайних промежутков в интервале. Для доказательства этого утверждения сформулируем простой критерий, позволяющий делать обоснованное заключение по поводу того, сохраняется знак в данном промежутке или нет. Пусть для промежутка [у1, у2] задана дуга параболы
Щ(х) = Л(х-у1)2 + Б(х-у1) + С, х е К у2],
такая, что ее концы имеют одинаковый знак: СЩу2 - у1) > 0.
Очевидно, если выполняется хотя бы одно из условий
б(2л(у2 - у1) + Б) > 0 или БС > 0,
то дуга параболы Щ(х) на всем интервале [у1, у2] не меняет знак. Если же эти условия не выполняются, то она не меняет знак, если
Я = (2б(у2 - у1) + С) С > 0. (7)
Применим критерий (5) для исследования поведения производной сплайна 5(х) на интервале [а, Ь]. Для этого сначала вычислим его значения для каждого промежутка [х,, ил], [ил, и,2], [и,2, хт], , = 0, ..., п - 1.
Я,1 = + У,И, /6, Яа = 3 % - 5', - 5+ + И(5ГМ - 5',)/ 9 = 3 $ - - 5,+ + И,(5',+1 - 5',)/ 9 = = 3$ - (5 , + 5' , И, /6 + 5 ,+1 - И, 5' ,+1/6) - И,(5'+ - 5' ,)/ 18, Я,3 = 5,+1- ,+1 И, /6.
Очевидно, Яй и Я,3 сохраняют знак разделенной разности 8, если для соответствующих крайних промежутков для первых и вторых производных сплайна в узлах исходной сетки выполняются условия (3-6).
Для сохранения знака и в центральном промежутке какого-либо интервала достаточно, чтобы выполнялись уже более жесткие условия для величины первых и вторых производных в узлах исходной сетки относительно разделенных разностей с>, и 8, соответственно, а именно:
| 5, | < 1,5 | $ |, | 5< 1,5 | Я|, 8,Ф 0, (8)
|5' ,| < 1,5 | |, | 5' ,+11 < 1,5 | ¿>,|, д1 Ф 0. (9)
Остается показать, что алгоритмы расчета первых и вторых производных в узлах исходной сетки, обеспечивающие выполнение для них условий (3, 4, 8, 9), существуют.
В качестве конструктивного доказательства приведем следующие формулы для расчета первой производной во внутренних узлах исходной сетки:
ад = { ^ ^ (10)
^ Ч,, иначе,
где ( 1,58,, 8, / 8,-1 < 1/3,
Ч, = { 28гЛ 8, / (8,-1 +8,), 1/3 < 8, / 8,-1 < 3, (11)
1,5^-1, 5-/ > 3.
Алгоритмы, аналогичные (10)-(11), можно предложить также и для вычисления вторых производных сплайна во внутренних узлах исходной сетки.
В свете изложенного очевидны ограничения на величину и знак производных на концах интервала [а, Ь], а именно:
^0 = у^0, 5\ = у8п -Ь 5'0 = У^0, 5'п = у8п _ь 0 < У < 1,5.
Следует отметить высокую степень сходства описанного сплайна с кусочно-линейной интерполяцией таблично заданной функции. Как и при аппроксимации подобных функций ломаными (для функций от одного аргумента) и плоскими прямоугольными участками (для функций от нескольких аргументов) все локальные экстремумы находятся в узлах исходной сетки, а не внутри интервалов. При аппроксимации же обычными (основными) сплайнами значения функций в узлах исходной сетки не всегда являются локальными экстремумами, а часто оказываются на "излете" кривой или поверхности сплайна.
4. Заключение
Таким образом, суммируя вышеизложенное, мы можем сформулировать следующие выводы:
Во-первых, по сравнению с традиционными схемами и алгоритмами аппроксимации и интерполяции функций, предложенная нами схема намного технологичнее, ибо для ее реализации не требуется держать в памяти компьютера все коэффициенты сплайнов для каждого промежутка.
Во-вторых, она хорошо приспособлена к разнообразию ситуаций и сохраняет наиболее важные изогеометрические свойства интерполируемых функций, зафиксированных табличной формой их представления: монотонность и знак.
В-третьих, сохраняя монотонность и знак интерполируемых функций, схема повышает точность аппроксимации, так как исключает на каждом интервале колебания сплайна и перемену знака (если значения функции на концах интервала не разного знака).
В-четвертых, обеспечивая единообразие в построении основных и дополнительных узлов, схема (за счет равномерной сходимости сплайна к интерполируемой функции и исключения критических ситуаций, требующих смещения и увеличения числа дополнительных узлов) позволяет строить многомерные сплайны для аппроксимации поверхностей и кривых в многомерных пространствах.
В-пятых, предложенная схема устойчива к краевым условиям, перепадам градиентов и неравномерности исходной сетки.
Теперь вернемся к задачам, перечисленным в начале статьи и опирающимся на алгоритмы сохранения изогеометрии.
Наиболее простой и естественной выглядит здесь задача восполнения данных. Для этой цели пригодна обычная схема интерполяции. Но в отличие от традиционных алгоритмов сплайн-интерполяции вышеописанная схема не только выдает более точные количественные результаты, но и сохраняет их структуру. Результаты могут быть представлены в форме гладкой кривой или поверхности.
Возможность построения поверхностей и кривых, сохраняющих на интервалах монотонность и знак интерполируемых функций, обеспечивает решение другой задачи - сохранения изогеометрической и топологической структуры данных.
Наконец, весьма привлекательной представляется возможность заменить все кусочно-линейные графики поверхностей и кривых на соответствующие гладкие.
5. Пример
Если рассмотреть графики сплайнов ^(х) и S2(x), интерполирующих функцию Д(х) = х2 по ее значениям в узлах сетки: х0 = 0.2, XI = 0.3, х2 = 1.0, х3 = 2.0, где ^(х) - классический интерполяционный квадратический сплайн, а 5"2(х) -изогеометрический сплайн (сохраняющий монотонность и выпуклость), то получим следующие графики (рис. 1).
Отметим, что несмотря на очень грубую сетку (всего 4 узла) результаты интерполяции сплайном £2(х) достаточно удовлетворительны.
Литература
Воронин В.Т. Построение сплайнов, сохраняющих изогеометрию. Препринт ВЦ СО АН СССР № 404.
Новосибирск, с.33-35, 1982. Воронин В.Т. Восстановление плотности распределения выборочных данных с помощью сплайнов, сохраняющих изогеометрию. Сплайн-функции в экономико-статистических исследованиях. Новосибирск, с.15-17, 1987. Стечкин С.Б., Субботин Ю.Н. Сплайны в вычислительной математике. М., Наука, с.77-78, 1976. Корнейчук Н.П. Сплайны в теории приближения. М., Наука, с.22-25, 1984.
Рис. 1. Интерполяция квадратическим (слева) и изогеометрическим (справа) сплайном
