Построение кубических сплайнов, сохраняющих выпуклость данных Текст научной статьи по специальности «Математика»

УДК 519.652

ПОСТРОЕНИЕ КУБИЧЕСКИХ СПЛАЙНОВ, СОХРАНЯЮЩИХ ВЫПУКЛОСТЬ ДАННЫХ

М.М. Ромаданова

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

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

Введение. Задача формосохраняющей интерполяции данных возникает при решении многих прикладных задач. При этом наиболее используемым способом интерполяции является построение кубических сплайнов. Способы построения таких сплайнов можно найти во многих монографиях, таких как [1-3]. Под формосохраняющей интерполяцией обычно понимают сохранение у сплайнов геометрических свойств, таких как положительность, монотонность и выпуклость [4]. Построение сплайнов с сохранением данных свойств рассматривалось многими авторами. Например, возможность монотонной интерполяции кубическими сплайнами была впервые рассмотрена в работе Мирошниченко В.Л. [5]. Его исследования основаны на доказанной лемме о неотрицательном решений трехдиагональной системы линейных уравнений. В работе [6] Квасовым Б.И. был предложен весовой кубический сплайн и алгоритмы, позволяющие подобрать весовые параметры так, чтобы сохранялась монотонность и выпуклость данных.

В работах [7, 8] были рассмотрены три алгоритма построения весового кубического сплайна, сохраняющего монотонность данных, и его применение при анализе геофизических данных. В статье [7] показано, что масштабирование данных по осям х и у упрощает построение монотонного весового кубического сплайна и выбор весовых параметров. Также при выборе определенных весовых коэффициентов масштабирование данных может гарантировать построение монотонного сплайна. Задача построения к -монотонного кубического сплайна рассматривалась Волковым Ю.С. и др. в работе [9]. Авторами был предложен единый подход для получения к -монотонного сплайна для к = 0,...,4, где к = 0 соответствует неотрицательности данных, к = 1 - монотонности данных, к = 2 - выпуклости данных.

В данной статье рассматриваются алгоритмы построения кубических сплайнов, сохраняющих выпуклость данных. Сначала приводится построение весового кубического сплайна и адаптивный алгоритм выбора весовых параметров, предложенный Квасовым Б.И. в [6]. Затем предлагается улучшенный алгоритм построения кубического сплайна, использующий минимизацию по модулю вторых производных. В конце приводится несколько численных примеров, на которых сравнивается реализация предложенного алгоритма с простым и весовым кубическим сплайном.

Построение весового кубического сплайна с сохранением выпуклости данных. Пусть на отрезке [а, Ь\ в узлах сетки а = хо < Х1 < к < хN+1 = Ь заданы значения некоторой функции

{хг Л), г = о,к,N +1. (1)

Введем обозначения

Ащ, х+1 ]=, Ъ = хм - х,, ,=0,...,N

и

8,/ = /[х/,х,-+1] - /[х,-1, х,], / =1,..., ^.

Данные (1) будем называть монотонными, если /[х,,х,+1 ]> 0, / = 0,...,N, и выпуклыми, если 8// > 0, / = 1,...,N.

Весовым кубическим сплайном с весами щ, > 0, / = 0,...,N называется функция $ такая, что

1) между узлами сетки А функция $ является кубическим многочленом;

2) $е Ск [а,Ь], к > 1;

3) щ-1 Г(х")= / = 1,..., N.

Будем считать, что кубический сплайн $ удовлетворяет условиям интерполяции

$(х, ) = /,, / = 0,..., N +1. (2)

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

I. Задание граничных значений первой производной:

$,(а) = /o, $'(Ь) = /к+1.

II. Задание граничных значений второй производной:

$ "(а) = /0, $ "(Ь) = +1.

Алгоритм построения весового кубического сплайна $ был предложен в работе [6]. Автор работы Квасов Б.И. также предложил адаптивный алгоритм выбора весовых параметров щ,, позволяющих сохранить выпуклость данных.

Введем обозначения:

М, = щ,-1 г(х")= щгГ(х+), , = 1,..., N,

М 0 = щ $ "(х0), MN+1 = wNS№ (х-+1).

Для х е [х,, х,+1 ] имеем

Ъ 2

$ (х) = / (1 - г) + /+1/ - г(1 - г)-!- [(2 - г )М, + (1 + г )М,+1 ],

6Wj

X Xj

где t = ■ '

hi

Так как S'(x- )= S '(x+), i = 1,..., N, получаем -Mi-i + 2

hi-1 . hi-1 , hiЛ

wi-1

^ wi-1 wi

h

Mi + -LMi+1 = 6dif, i =N. (3)

Wi

виде:

Уравнения (3) вместе с краевыми условиями II типа могут быть переписаны в

2M о = do

m-Mi-1 + 2Mi +1 iMi+1 = di, i = 1,..., N

2MN+1 = dN+1 178

где ¡тг и %г имеют вид

1 =

™г-Л + щЬ-1

1 =, ¡г = 1 -1г, (4)

и правые части

^ = 2^о/о , dN+1 = 2wNfN+1, ^ = ^ -1^г .

^Л-1 + ^-1И1

Далее для построения алгоритма выбора весовых параметров wг предположим, что данные являются выпуклыми, т.е. 8г/ >о, г = 1,...,N. Если выполняются следующие ограничения на краевые условия:

о</о'<^. о</N+1 <^

ко hN

и на весовые параметры

^о < 6 8/ -1, Щ к о~ h о /о ' 1 8г/ < ^-1 ^ + 1£ 2 8г/

+1 <_—, г = 2, к, N, (5)

21 -1 8г^ к-1 ^г-1 8г-1/

WN hN-1 < 6 8 Nf - 1,

wN-1 hN +1

то S"(х)> о для всех хе [а,Ь\, т.е. S(х) является выпуклой функцией на отрезке [а,Ь\. Однако, выполнение этих неравенств возможно только для слабо меняющихся данных [6].

Предположим, что на участке выпуклости данных известны первые два параметра wг - 2 и wг-1. Обычно они задаются по формуле [4, 6]:

wг =(1 + Сг/[хг,хг+1\2)-Рг, Сг > 1, Рг >о, г = о,...,N.

Если неравенства (5) выполняются для wг = Wг-1, то полагаем wг = Wг-1. В противном случае wг находим из невыполненного неравенства (5) путем замены в нем

знака неравенства на знак равенства. Для wг < е (wг > е 1) полагаем wг = е (wг = е 1), где е - достаточно малое положительное число. Алгоритм начинаем с wо и wl и находим все параметры wг, обеспечивающие выпуклость весового кубического сплайна для любых выпуклых данных.

Построение кубического сплайна с минимизацией по модулю вторых производных. Пусть на отрезке [а, Ь\ в узлах сетки а = хо < х1 < ... < XN+1 = Ь заданы значения некоторой функции (хг-,/г), г = о,к, N + 1. Будем считать, что на каждом подот-резке [хг-, х/+1 \ заданы значения функции (хг-, /г) и (хг-+1, /г+1), и введем обозначения

М 2г = ^ '(*+) и М 2г+1 = 5 '(х"+1), г = о, к, N: г = о: (хо,/о), (хь/1), Мо, М1 г = 1: (х1, /1), (х2, /2), М2, М3

г: (хг, /г), (хг+l, /г+1), М2г, М2г+1

г = N: (xN, ^), (xN+1, ^+1), М2N, М2N+1

179

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

Г(х)° #(х)= М 2, (1 - г) + М 2,+1г, х е [х,, х,+1 ], (6)

г = , ъ = х,+1 - х,, , = N. Ъ,

Повторное интегрирование формулы (6) и использование условий интерполяции (2) (х, ) = /, и (х,+1 ) = /,+1 дает форму кубического сплайна на подотрезке

[х,, х,+1]:

Ъ 2

$ (г) = /, (1 - г) + /,+1г + г(1 - г[(г - 2М2, -(1 + гМ2,+1 ] , = 0- ■ -N

6 '

Для нахождения неизвестных коэффициентов М2, и М2,+1, , = 0,...,N будем использовать непрерывность первой производной сплайна. Дифференцирование формулы (6) дает:

+ +1 + Ъ Ъ Ъ 6 2'у ' ;6

Из условия £у_1(х, - 0) = $(х, + 0), , = 1,...,N получаем систему линейных алгебраических уравнений

АМ = В, (7)

где матрица коэффициентов системы А ~ NX 2^ +1) имеет вид:

гЪ-0 Ь-0 Ъ-1 Ь-1

(г) = -/ + /+1 + ^М2, (1 - 3(1 - г)2) + ~~М2,+1 (зг2 -1)

А =

6 3 3 6

Ъ Ь-1 ^2 ^2

0 0 ^ ^ ^ ^

63 0 0 0 0

0000 00

36

¿2 ¿2 Ъ3 Ъ3

0 0 0

0 0 0 0 0

6336

Ъ Ъ Ъ+1

63

3

0 0 0

Ъ+1 6

0 0 0

0

0 0 0 0 0 0 0

%-1 %-1 % % 3 6 у

63

вектор неизвестных М ~ 2^ + 1)х 1 и вектор правых частей В ~ N X1 имеют вид:

(

М =

С М 0 Л М1 М 2 М 3

М 2 N

М 2 N+1,

В =

У0 Ъ0

У1 ¿1

с

г

1 1

— + —

Ъ0 Ъ1

Л

У1 +

У2

Л

Ъ

1

(

11 — + —

¿1 ¿2

У + У3

У2 +-Г-Ъ2

11

— + ■

У,

Ъ VЪ, Ъ,+1

Л

УN-1 %-1

11

-+

У,+1 +

Л

V ^-1

Ъ

N

УN +

У,+2 Ъ+1

yN+1

Ъ

N

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

2N+1 2

задачу минимизации ^ (Mj )2 ® min. С помощью псевдообратной матрицы можно

i=0

найти минимальное по норме решение системы уравнений [10].

Для определенности также будем считать, что данные выпуклые, т.е. все M j,

j = 0,...,2N +1 должны быть положительны. Решая систему уравнений (7), получим вектор вторых производных M . Если в решении какие-либо Mj получились отрицательными, то добавляем дополнительные ограничения, состоящие в том, что M j = 0 .

Таким образом, в матрицу системы A добавляются дополнительные строки, в которых в j столбце стоит единица, а все остальные ноли. Соответственно в вектор правых частей B добавляются нулевые строки. Далее находится минимальное по норме решение дополненной системы уравнений. Если в решении какие-либо M j снова получились

отрицательными, то аналогично добавляем ограничения M j = 0 и находим решение

системы уравнений. Процесс итераций заканчивается, когда все Mj > e,

j = 0,k,2N +1, где e - достаточно маленькое число.

Численные примеры. Для первого примера рассмотрим построение кубического сплайна для средних значений температуры в пределах высот 0-3000 м. для многолетнего января в 12 часов по станции Нижний Новгород за 1964-2010 гг. [11, 12]. На рис. 1 приведен построенный по данным значениям кубический сплайн по предложенному алгоритму. Весовой кубический сплайн дает практически такой же результат.

Рис. 1. Наблюдаемые и проинтерполированные кубическим сплайном с минимизацией вторых производных средние значения температуры

Для более наглядного примера рассмотрим интерполяцию функции

ехр(100 х)~1 х е[од]

f (x) = 1 --

exi

p(100)-1 '

на равномерной сетке xi

i

10

, i = 0,1,...,10. Данный численный пример рассматривал-

ся во многих работах, где приводилось сравнение построения обычного кубического сплайна с рациональным сплайном [1], с весовым кубическим сплайном, сохраняющим монотонность и выпуклость данных [6, 13]. Здесь, на рис. 2 (а) и (б), приведем построение кубического сплайна (весовые параметры = 1, краевые условия /0 = 0, fN+1 = 0) и сплайна, построенного с минимальными вторыми производными по предложенному алгоритму.

а

б

Рис. 2. Построение (а) С -кубического сплайна и (б) сплайна с минимальными

вторыми производными

Для следующего примера рассмотрим интерполяцию средних значений площади снежного покрова р. Невы за 2000-2009 гг. Как было показано в работе [14], при интерполяции кубическим сплайном на интервале [129, 193] функция принимает отрицательные значения, что противоречит физическому смыслу. В работе показано, что применение весового кубического сплайна позволяет сохранить положительность функции. На рис. 3 для тех же экспериментальных данных аналогично приведем сравнение построения кубического сплайна с построением сплайна с выбором минимальных вторых производных.

а

2

б

Рис. 3. Построение (а) С -кубического сплайна и (б) сплайна с минимальными вторыми производными для средних значений площади снежного покрова р.

Невы за 2000-2009 гг.

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

182

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

Список литературы

1. Завьялов Ю.С., Квасов Б.И., Мирошниченко В. Л. Методы сплайн функций. М.: Наука, 1980. 352 с.

2. Квасов Б.И. Методы изогеометрической аппроксимации сплай-нами. М.: Физматлит, 2006. 360 с.

3. Вагер Б.Г., Серков Н.К. Сплайны при решении прикладных задач метеорологии и гидрологии. Л.: Гидрометеоиздат, 1987. 160 с.

4. Волков Ю.С. О монотонной интерполяции кубическими сплайнами // Вычислительные технологии. 2001. Т. 6. № 6. С. 14 - 24.

5. Мирошниченко В. Л. Достаточные условия монотонности и выпуклости для интерполяционных кубических сплайнов класса C2 // Вычислительные системы: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ИМ. 1990. Вып. 137: Приближение сплайнами. С. 31 - 57.

6. Kvasov B.I. Monotone and convex interpolation by weighted cubic splines // Computational Mathematics and Mathematical Physics. 2013. V. 53. № 10. P. 1428 -1439.

7. Ромаданова М.М. Алгоритмы построения монотонного весового кубического сплайна // Известия Тульского государственного университета. Технические науки. 2018. Вып. 9. С. 180 - 192.

8. Ромаданова М.М., Вагер Б.Г. Методы обработки эксперимен-тальных данных при моделировании геофизических процессов // Системы. Методы. Технологии. 2018. № 2 (38). C. 70 - 75.

9. Волков Ю.С., Богданов В.В., Мирошниченко В.Л., Шевалдин В.Т. Формосо-храняющая интерполяция кубическими сплайнами // Матем. заметки. 2010. Т. 88. Вып. 6. С 836 - 844.

10. Алберт А. Регрессия, псевдоинверсия и рекуррентное оценивание. М.: Наука, 1977. 224 с.

11. Алдухов О.А. О методах расчета и контроля данных по пограничному слою атмосферы // Труды ФГБУ «ВНИИГМИ-МЦД». 2012. Вып. 176. С. 257 - 286.

12. Алдухов О.А., Черных И.В. Методы анализа и интерпретации данных радиозондирования атмосферы. Т.1. Контроль качества и обработка данных. М-во природных ресурсов и экологии Российской Федерации, Федеральная служба по гидрометеорологии и мониторингу окружающей среды (Росгидромет), Всероссийский науч.-исследовательский ин-т гидрометеорологической информ. Мировой центр данных. Обнинск, Калужская обл.: ВНИИГМИ-МЦД. 2013. 306 с.

13. Kim Tae-wan, Kvasov B.I. A shape-preserving approximation by weighted cubic splines // J. Comput. Appl. Math. 2012. V. 236. P. 4383 - 4397.

14. Петрова А.В., Вагер Б.Г. Весовые кубические и бикубические сплайны при анализе гидрометеорологических наблюдений // Труды Главной геофизической обсерватории им. А.И. Воейкова. 2016. № 583. С. 149 - 161.

Ромаданова Мария Михайловна, канд. физ.-мат. наук, доцент, Romadanova@yandex. ru, Россия, Санкт-Петербург, Санкт-Петербургский государственный архитектурно-строительный университет

CONSTRUCTION OF CUBIC SPLINES PRESERVING CONVEXITY OF THE DATA

M.M. Romadanova

The algorithms for constructing cubic splines that preserve the convexity of the data are considered in this article. We construct a weighted cubic spline and an adaptive algorithm for selecting weight parameters. Then, an improved algorithm is proposed. This algorithm for constructing a cubic spline is based on minimization of sum of squared second order derivatives. The efficiency of proposed algorithm in comparison with simple and weighted cubic splines is demonstrated on several examples.

Key words: shape-preserving interpolation, convex interpolation, weighted cubic spline, geophysical experimental data.

Romadanova Mariya Mikhailovna, candidate of physico-mathematical sciences, do-cent, Romadanova@yandex.ru, Russia, Saint Petersburg, Saint Petersburg State University of Architecture and Civil Engineering

УДК 519.8

АЛГОРИТМ ОБРАБОТКИ ТРАЕКТОРНОЙ ИНФОРМАЦИИ УТОЧНЕНИЯ

ПАРАМЕТРОВ ДВИЖЕНИЯ ОБЪЕКТОВ ПО ОДНОВРЕМЕННЫМ ДВУКРАТНЫМ ИЗМЕРЕНИЯМ ДВУМЯ КОСМИЧЕСКИМИ АППАРАТАМИ

РОТОРНОГО ТИПА

Е.П. Минаков, Р.П. Власов

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

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

Необходимость комплексного решения задачи контроля космического пространства, получения и уточнения данных о параметрах движения объектов космического мусора (ОКМ) техногенного и естественного происхождения, для предотвращения воздействия на функционирующие космические аппараты (КА), делают актуальным решение задачи разработки алгоритма уточнения параметров движения ОКМ по одновременным двукратным измерениям двумя КА роторного типа (РТ). Это подтверждается тем, что анализ существующих подходов к решению задач, возникающих при решении указанной проблемы, продемонстрировал их вербальный характер, отсутствие системности, глубокой научной и тем более инженерной проработки [1, 6, 7].

Постановка задачи и циклограмма применения КА РТ

Пусть кластер КА РТ представляет собой совокупность двух КА РТ, оснащенных аппаратурой измерения параметров ОКМ.