CC BY
67
УДК 519.652
АЛГОРИТМЫ ПОСТРОЕНИЯ МОНОТОННОГО ВЕСОВОГО КУБИЧЕСКОГО СПЛАЙНА
М.М. Ромаданова
Рассматриваются алгоритмы построения весового кубического сплайна, сохраняющего монотонность данных. Первый алгоритм основан на разбиении данных на интервалы возрастания и убывания. Следующие два алгоритма используют масштабирование данных. В статье показано как данное преобразование позволяет построить сплайн сразу на всем интервале изменения данных, а также может упростить алгоритм выбора весовых параметров. Результаты работы предложенных алгоритмов показаны на примере интерполяции геофизических данных.
Ключевые слова: весовой кубический сплайн, монотонная интерполяция, зональная составляющая скорости ветра, меридиональная составляющая скорости ветра.
Введение. При решении многих прикладных задач, связанных с интерполяцией данных, используются кубические сплайны. Классические кубические сплайны являются одним из наиболее распространенных средств вычислительной математики. Описание построения таких сплайнов можно найти во многих монографиях, например, в [1, 2]. Однако при решении практических задач классические кубические сплайны, которые являются дважды непрерывно дифференцируемой функцией, могут не сохранять геометрические свойства исходных данных, такие как положительность, монотонность или выпуклость. Для выполнения требуемых свойств функции было разработано множество модификаций классических кубических сплайнов. Однако при этом может не сохраняться какое-либо из достоинств классического сплайна: гладкость, точность, порядок приближения, простота реализации и др. [3].
Возможность монотонной интерполяции обычными кубическими сплайнами была впервые рассмотрена в работе Мирошниченко В. Л. [4]. Его исследования основаны на доказанной лемме о неотрицательном решений трехдиагональной системы линейных уравнений.
Проблемы построения монотонных и выпуклых классических кубических сплайнов были рассмотрены Квасовым Б.И. в работе [2], а в [5] им был предложен весовой кубический сплайн, который позволяет влиять на значение второй производной интерполируемой функции, однако не сохраняет её непрерывность. Весовой кубический сплайн также можно использовать, когда требуется сохранить свойства монотонности и выпуклости функции. Алгоритмы построения монотонного и выпуклого весового кубического сплайна подробно изложены в [5].
180
В данной статье мы рассмотрим и подробно опишем несколько алгоритмов построения весового кубического сплайна, сохраняющего монотонность данных. Один из алгоритмов будет основан на разбиении интервала задания функции на подынтервалы. Другие алгоритмы - на масштабировании данных по осям х и у, что упрощает построение сплайна и выбор весовых параметров.
Численное тестирование полученных алгоритмов будет рассмотрено на примере интерполяции гидрометеорологических данных. В метеорологии и гидрологии задача интерполяции данных играет большую роль, т.к. от точности интерполяции исходных данных зависят последующие прогнозы [6]. Возможности использования сплайнов при анализе гидрометеорологической информации рассмотрены во множестве работ [7, 8]. В работе [9] авторами показано, что весовой кубический сплайн дает хорошую аппроксимацию профилей скорости ветра, которая может быть использована для вычисления коэффициента турбулентности в пограничном слое атмосферы.
Постановка задачи. Пусть на отрезке [а,Ь] в узлах сетки а = хо < Х1 < к < = Ь заданы значения некоторой функции
(х /), I=0,к,М +1. (1)
Введем обозначение
/Ь,хI+1] = /г+\ / , ИI = хI+1 - х^, I = 0,к,N.
Данные (1) будем называть монотонными, если /[х^, хг+1 ]> 0, I = 0,...,N.
Задача формосохраняющей интерполяции [5] состоит в построении достаточно гладкой функции £ такой, что £(х^ ) = /, I = 0,...,N +1, и £ монотонна на участках монотонности исходных данных.
Весовым кубическим сплайном с весами >0, I = 0,...,N называется функция £ такая, что
1) между узлами сетки А функция £ является кубическим многочленом;
2) £е Ск [а,Ь], к > 1;
3) wi-1 £"(х-)= wl£*х+), I = 1,.,N.
Будем считать, что кубический сплайн £ удовлетворяет условиям интерполяции
) = /, I = 0,..., N +1. (2)
Для однозначного определения сплайна нам также потребуются краевые условия. Наиболее часто используемыми являются ограничения следующих типов.
I. Задание граничных значений первой производной: £'(а) = /0, ^ \Ъ ) = +1.
II. Задание граничных значений второй производной: £ '(а ) = /0,
г (ь )=л+1.
Построение весового кубического сплайна. Введем обозначения
£ '(х, ) = т,, I = 0,..., N +1. (3)
Сплайн £ можно рассматривать как эрмитов кубический сплайн [1], удовлетворяющий условиям (2)-(3). Тогда для хе [х,,х,+1 ] получаем
£(х)= /г (1 - г)2 (1 + 2*)+ / +1 г2 (3 - 2*)+ т* г(1 - г)2 - тг+1*4 г2 (1 - г),
Х Х,
где г =--.
Кубический сплайн, представленный в таком виде на каждом из промежутков [х,,х,+1 ], непрерывен вместе со своей первой производной
всюду на
[а, Ь]. Используя условие гладкости -1 £"(х, )= £"(х+), I = 1, к, N, получаем
1 ¡ту-1 + 2т,- +тгтг+1 = с,, I = 1,..., N, (4)
где
с, = 3( т/[х|, х,+1 ] +11/[х,-1, х, ]) (5)
и
1 =—^А , т = 1 -1. (6)
wi-ф, + щ*,-1
К уравнениям (4) необходимо добавить уравнения, вытекающие из краевых условий. Таким образом, получается система N + 2 уравнений для определения N + 2 неизвестных т,, I = 0,., N +1. В случае краевых условий I типа находим
2т0 = С0
1 ¡т-1 + 2т + +1 = с,, -=. ^N (7)
2mN+1 = СN+1
где С0 = 2/0, СN+1 = 2/N+1 и коэффициенты с,, 1, и т, имеют вид (5) и (6) соответственно.
Рассмотрим пример построения кубического сплайна при интерполяции гидрометеорологических данных. Для численного моделирования были выбраны данные профилей скорости ветра и и V в пределах высот 0-3000 м. для многолетнего июля в 12 часов по станции Нижний Новгород за 1964-2010 гг. [10, 11]. На рис. 1 (а) приведены значения зональной составляющей скорости ветра и, а также построенный по ним С -кубический сплайн. Все коэффициенты Wi были взяты равными 1.
На рис. 1, б выделенный фрагмент данного графика показывает, что при построении сплайна на этом интервале не сохраняется монотонность данных.
На рис. 2, а показаны значения меридиональной составляющей скорости ветра V и построенный по ним С2 -кубический сплайн, и на рис. 2, б, в, г - выделенные фрагменты, показывающие, что при построении сплайна также не везде сохраняется монотонность данных.
а
1.4
1200 1400 Высота z(m)
б
Рис. 1. Наблюдаемые и проинтерполированные кубическим сплайном
значения зонального ветра u
а
1600 1700 1800 Высота г (м) в
б
1900 2000 2100 Высота г (м) г
Рис. 2. Наблюдаемые и проинтерполированные кубическим сплайном значения меридионального ветра V
Условия монотонности весового кубического сплайна. В работах [2, 5] были получены достаточные условия монотонности интерполяционного С -кубического сплайна.
Теорема. Пусть весовой кубический сплайн £ с краевыми условиями £ '(хо) = /о и £'(.хN+1) = /ы+1 интерполирует монотонные данные {/} , I = 0,. к N +1. Если выполнены условия
0 £ /0 £3/[хо,Х1 ], 0 £ /Ы+1 £ 3/[хы, хы+1],
11/[х1 -1, X] £ (1 +11)/ [хг, х1+1], (8)
184
т,/[х, , х,+1 ]<(1 + т,)/ [х,-1, х, ], -=1,..., N, (9)
то £ (х)> 0 для всех х е [а, Ь], то есть функция £ монотонная на [а, Ь].
Учитывая формулы для коэффициентов 1, и т, вида (6), неравенство (8) можно переписать в виде:
Щ *-1 > /[х,-l,Х,] - 2
Wi-1 * / [х,, х,+1] (10) и неравенство (9) - в виде
Щ-1 * > / [х,, Х,+1] - 2
Щ, *-1 / [х,-1, х, ] . (11)
Алгоритмы построения весового кубического сплайна, сохраняющего монотонность данных. Опишем три алгоритма построения весового кубического сплайна и выбора весовых коэффициентов щ, , позволяющих сохранять монотонность данных. Для этого рассмотрим следующий выбор весовых параметров [3, 5]:
=(1 + Сг/[хг,х,+1 ]2)-Рг, С, > 1, рг >0, , = 0,...,N. (12)
Заметим, что если выбрать Ь, = 0 для всех I = 0,.. N, то все параметры щ,
2
будут равны 1, и мы получим С -кубический сплайн. Алгоритм 1
1. Выделим интервалы возрастания и убывания функции.
2. Зададим первый весовой параметр Щ0 по формуле (12), выбрав какие-либо коэффициенты С0 > 1 и Р0 >0.
/ - /
3. Зададим значение первой производной: при , = 0: /0 = —-0,
х1 - х0
при , = N +1: /^+1 = ^^+1—. Значение первой производной на грани-
ХN+1 - ХN
цах интервалов монотонности возьмем равным /^ = 0.
4. На каждом участке монотонности данных выберем первый параметр щ,-1. При этом для стыковки сплайнов первый параметр щ,-1 берем равным последнему на предыдущем участке. Если функция на интервале является монотонно убывающей, то будем считать, что мы работаем с - /, .
Далее выбираем щ, равное предыдущему значению щ,-1. Если неравенства (10) и (11) при этом выполняются, то полагаем щ, = щ,-1. Если одно из неравенств не выполняется, то вычисляем щ, из того неравенства, которое не выполняется при щ, = щ,-1, заменяя знак неравенства на знак равенства. Таким образом, если не выполняется неравенство (10), то
Щ = Щ-1
/-1
/ [х/-1, х/ ] / [х/ , X/+1]
2
и если не выполняется неравенство (11), то
Щ/ = Щ
/ [х/,х/+1]
Л-1
2
/ -11/ [х/-l, х/]
Если на каком-либо шаге получается, что щ < е, то полагаем щ = е, где £ - достаточно малое положительное число. Таким образом, находим все весовые параметры щ на каждом интервале возрастания или убывания функции.
5. Находим решение системы уравнений (7) и строим весовой кубический сплайн £ отдельно на каждом интервале монотонности.
Проиллюстрируем выполнение этого алгоритма для интерполяции гидрометеорологических данных, приведенных на рис. 1 и 2. На рис. 3 (а) построен график значений зональной составляющей скорости ветра и и построенный по ним весовой с '-кубический сплайн по первому алгоритму с адаптивным выбором весовых параметров щ. На рис. 3 (б) - вырезанный фрагмент этого графика. Сравнивая его с рис. 1 (б), мы видим, что в этом случае монотонность данных сохраняется.
На рис. 4, а построены значения меридиональной составляющей
скорости ветра V и построенный по ним весовой С1-кубический сплайн. На рис. 4, б, в, г также показаны приближенные фрагменты этого графика, аналогичные рис. 2, б, в, г. Как видно на рисунках, построенные по данному алгоритму сплайны сохраняют монотонность данных, в отличие от классического кубического сплайна.
а б
Рис. 3. Иллюстрация построения весового кубического сплайна значений зонального ветра и по 1-му алгоритму
_2. (М)_
Проинтерполированные значения v
Меридиональная составляющая скорости ветра м
а
Высота г(м)
б
Рис. 4. Иллюстрация построения весового кубического сплайна значений меридионального ветра V по 1-му алгоритму
Для упрощения этого алгоритма сделаем преобразование координат так, чтобы сетка стала равномерной по х и значения функции возрастали, т.е. / [х!, х!+1 ]> 0. Этого можно добиться, выполнив преобразование координат по осям х и у . Таким образом, мы сможем построить весовой кубический сплайн сразу на всем интервале. Затем, выполнив обратное преобразование координат, получим значения сплайна в исходной системе координат.
Алгоритм 2
1. Выполним преобразование данных по оси х так, чтобы они стали равномерными, например, X^ = 1,2, к.
2. Выполним преобразование данных по оси у так, чтобы они стали возрастающими, т.е. f [х^,х^+1 ]> 0, / = 0,....N.
3. Зададим первый весовой параметр w0 по формуле (12), выбрав какие-либо коэффициенты Со - 1 и bo -0.
4. Для i = 1,...,N сначала выбираем Wj равное предыдущему значению Wi_i. Аналогично предыдущему алгоритму, если неравенства (10) и (11) при этом выполняются, то полагаем Wj = WjЕсли одно из неравенств не выполняется, то вычисляем Wi из того неравенства, которое не выполняется при Wj = Wj_1, заменяя знак неравенства на знак равенства. Если на каком-либо шаге получается, что Wj < e, то полагаем Wj = e, где £ - достаточно малое положительное число. Таким образом, вычисляем все весовые параметры Wj.
5. Находим решение системы уравнений (7) и строим весовой кубический сплайн S по преобразованным данным.
6. Выполним обратное преобразование узловых и промежуточных точек сплайна.
Покажем построение весового кубического сплайна по этому алгоритму на том же численном примере. Сначала выполняется масштабирование данных так, чтобы они стали возрастающими. На рис. 5 (а) и (б) точками показаны преобразованные данные профилей скорости ветра u(z) и v(z) - Yu (X) и Yv (X) соответственн°. Далее по преобразованным данным строится весовой кубический сплайн на всем интервале. Так как данные приведены к возрастающим, то задав начальный весовой параметр W0 , остальные Wj j = 1,...,N берем либо равными Wj_1, либо находим из одного из неравенств (10) или (11). На рис. 5 (а) и (б) сплошной линией показаны построенные весовые кубические сплайны. После этого выполняется обратное преобразование значений узловых и промежуточных точек сплайна, результат которого показан на рис. 5 (в) и (г).
Для описания следующего алгоритма будем предполагать, что сетка является равномерной по x и f [xj, Xj+1 ]> 1 для всех j. Эти условия также могут быть обеспечены путем масштабирования данных по x и y. Если при этом f [xj, Xj+1 ]> f [xj_1, Xj ], то неравенство (10) будет выполняться для всех j. Неравенство (11) с учетом выбора весовых параметров по формуле (12) и того, что данные по оси x являются равномерными, примет вид
(13)
Обозначим а/ = /[х/, х/+1 ]. В частном случае, когда С/-1 = С/ = 1 и Ь/-1 = Ь/ = 1, неравенство (13) примет вид
188
1 + а 2
х а
>
а;
2.
1 + а/2_1 а I-1
б
1000 2000 3000 0 1000 2000 3000
г (м) т. (м)
в г
Рис. 5. Иллюстрация построения весового кубического сплайна
по 2-му алгоритму
Выполнение этого неравенства будет эквивалентно выполнению условия
g (а I _1, а I) = (1 + а 2) а | _1 - (1 + а 2_1) а | + 2а | _1 (1 + а 2_1 )> 0. Функция g(а| _1, а|) принимает минимальное значение по а |, когда
/ * 1 ( ga = 0. Решением этого уравнения будет а| = — а|_1 +
1 21
g (а | _1, а |) в точке минимума принимает значение
а
и функция
I-1)
1
g(а/-1,а/ )= 4 7а3-1 + 10а/-1
1
а,
> 0.
V -1 у
Тогда неравенство (13) также выполняется, и весовой кубический сплайн будет сохранять монотонность данных.
Таким образом, данное преобразование данных и выбор весовых параметров щ по формуле (12) при С/ = 1 и Ь/ = 1 гарантирует построение монотонного сплайна.
Алгоритм 3
1. Выполним так же, как во 2-м алгоритме, преобразование данных по оси х так, чтобы они стали равномерными X/ = 1,2, к.
2. Выполним преобразование данных по оси у так, чтобы выполнялись условия / [хо, х1 ]> 1 и / [х/, х/+1 ]> / [х/-1, х/ ], / = 1,..., N.
3. Возьмем коэффициенты С/ = 1 и Ь/ = 1 и найдем все весовые параметры Щ/, / = 0,...,N.
4. Построим весовой кубический сплайн £ по преобразованным данным.
5. Выполним обратное преобразование узловых и промежуточных точек сплайна.
а
б
1000 2000 г (м)
3000
в г
Рис. 6. Иллюстрация построения весового кубического сплайна
по 3-му алгоритму
190
Проиллюстрируем реализацию третьего алгоритма на тех же численных данных профилей скорости ветра. Аналогично предыдущему примеру на рис. 6 (а) и (б) показаны преобразованные данные профилей скорости ветра и(г) и у(г) - Уи (X) и Уу (X) соответственно, таким образом, чтобы выполнялись пункты 1 и 2 третьего алгоритма. Для построения весового кубического сплайна выбираем коэффициенты С| = 1 и Ь = 1. Так как для преобразованных данных неравенства (10) и (11) будут выполняться, то легко находим все весовые параметры , I = 0,...,N по формуле (12) и строим весовой кубический сплайн, который показан на этих рисунках сплошной линией. Далее выполняем обратное преобразование значений узловых и промежуточных точек сплайна, результат которого показан на рис. 6 (в) и (г).
Заключение. В работе построены три алгоритма интерполяции весовым кубическим сплайном, сохраняющим монотонность данных.
Первый алгоритм основан на том, что выделяются интервалы возрастания и убывания функции. Значения весовых параметров на границах интервалов при этом выбираются так, чтобы _1 £"(х_)= м^Б*(х+). Далее отдельно на каждом интервале монотонности строится весовой кубический сплайн таким образом, чтобы в результате сохранялась непрерывность первой производной.
Во втором алгоритме исходные данные масштабируются таким образом, чтобы они стали возрастающими. Это позволяет не разбивать данные на промежутки возрастания и убывания, а строить сплайн сразу на всем интервале.
В третьем алгоритме данные не только приводятся к возрастающим, но и добавляется условие, чтобы f [х0, х1 ]> 1 и на каждом следующем интервале выполнялось f[х,х^+1 ]>f[х_1,х^ ], I = 1,...,N. Это позволяет упростить вычисление всех весовых параметров и не разбивать данные на промежутки возрастания и убывания.
Список литературы
1. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн функций. М.: Наука, 1980. 352 с.
2. Квасов Б.И. Методы изогеометрической аппроксимации сплайнами. М.: Физматлит, 2006. 360 с.
3. Волков Ю.С. О монотонной интерполяции кубическими сплайнами // Вычислительные технологии. 2001. Т. 6. № 6. С. 14 - 24.
4. Мирошниченко В.Л. Достаточные условия монотонности и выпуклости для интерполяционных кубических сплайнов класса С2 // Вычислительные системы: Сб. науч. тр. / АН СССР. Сиб. отд-ние. ИМ. 1990. Вып. 137: Приближение сплайнами. С. 31 - 57.
191
5. Kvasov B.I. Monotone and convex interpolation by weighted cubic splines // Computational Mathematics and Mathematical Physics. 2013. V. 53. № 10. P. 1428 - 1439.
6. Петрова А.В., Вагер Б.Г. Весовые кубические и бикубические сплайны при анализе гидрометеорологических наблюдений // Труды Главной геофизической обсерватории им. А.И. Воейкова. 2016, № 583. С. 149 -161.
7. Марчук Г.И. Математическое моделирование в проблеме окружающей среды. М.: Наука, 1982. 319 с.
8. Вагер Б.Г., Серков Н.К. Сплайны при решении прикладных задач метеорологии и гидрологии // Л.: Гидрометеоиздат, 1987. 160 с.
9. Ромаданова М.М., Вагер Б.Г. Методы обработки экспериментальных данных при моделировании геофизических процессов // Системы. Методы. Технологии. 2018. № 2 (38). C. 70 - 75.
10. Алдухов О. А. О методах расчета и контроля данных по пограничному слою атмосферы // Труды ФГБУ «ВНИИГМИ-МЦД». 2012. Вып. 176. С. 257 - 286.
11. Алдухов О.А., Черных И.В. Методы анализа и интерпретации данных радиозондирования атмосферы. Т.1. Контроль качества и обработка данных / М-во природных ресурсов и экологии Российской Федерации, Федеральная служба по гидрометеорологии и мониторингу окружающей среды (Росгидромет), Всероссийский науч.-исследовательский ин-т гидрометеорологической информ. Мировой центр данных. Обнинск, Калужская обл.: ВНИИГМИ-МЦД. 2013. 306 с.
Ромаданова Мария Михайловна, канд. физ.-мат. наук, доцент, Romadanova@yandex. ru, Россия, Санкт-Петербург, Санкт-Петербургский государственный архитектурно-строительный университет
ALGORITHMS TO CONSTRUCT THEMONOTONIC WEIGHTED CUBIC SPLINE
M.M. Romadanova
We consider algorithms for constructing a weighted cubic spline that preserves the monotonicity of the data. The first algorithm is based on splitting the data into intervals of increasing and decreasing. The following two algorithms use data scaling. This transformation allows to construct a spline on the entire interval of data change. It also simplifies the algorithm of choosing weight parameters. The application of proposed algorithms is demonstrated on the example of interpolation of geophysical data.
Key words: weighted cubic spline, monotonic interpolation, zonal component of wind speed, meridional component of wind speed.
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
