Полулокальные сглаживающие сплайны седьмой степени Текст научной статьи по специальности «Математика»

УДК 519.6, 517.9

ПОЛУЛОКАЛЬНЫЕ СГЛАЖИВАЮЩИЕ СПЛАЙНЫ СЕДЬМОЙ СТЕПЕНИ

Д.А. Силаев, Ж.Г. Ингтем

SEMILOCAL SMOOTHIHG SPLINES OF SEVENTH DEGREE

D.A. Silaev, J.G. Ingtem

Полулокальные сглаживающие сплайны или S-сплайны были введены Д.А. Силаевым. Ранее рассматривались и применялись сплайны 3-й и 5-й степени. Настоящая работа посвящена построению сплайнов 7-й степени, доказаны теоремы существования и единственности, установлены условия устойчивости таких сплайнов.

Ключевые слова: аппроксимация, сплайн, численные методы

Semilocal smoothing splines or S-splines were initiated by D.A. Silaev. Earlier there were considers and applied the splines of third and fifth degree. This work is devoted to seventh degree splines construction. Uniqueness and existence theorems are proved. Stability and convergence conditions for these splines are established.

Keywords: an approximation, a spline, numerical methods

Введение

Рассматривается задача восстановления функции с помощью полулокального сглаживающего сплайна или S-сплайна, состоящего из полиномов седьмой степени. Здесь будут рассмотрены S-сплайны разных порядков гладкости: классов С0, С1, С2, С3 и С4.

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

1. Построение S'-сплайна седьмой степени

Рассмотрим на отрезке [а, Ь] равномерную сетку с узлами — а + kh, к =

0,1,... ,К, где h = ——--шаг сетки. Разобьем отрезок [а, Ь] на группы. Для этого введем

еще одну равномерную сетку {£z}i=o с узлами £; = а + Ш, где Н = mh, meN,l = 0,l,...,i. Таким образом, при переходе из одной группы в другую, будем осуществлять сдвиг системы координат и рассматривать каждый I-й полином на отрезке [0, if].

Пусть у = (уо, уъ ..., Ук) G значения приближаемой функции на сетке {xk}\z.Q .

Рассмотрим случаи I - V, когда S-сплайн принадлежит классу С0, С1, С2, С3,С*4.

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

Обозначим через

Р| = {и : и(х) = ао + а\х + а,2Х2 + а^х3 + &4Ж4 + а^х5 + а^х6 + а7х7}

множество полиномов 7-й степени с фиксированными коэффициентами ао, а*, а2 .(Здесь зафиксированы только три коэффициента, так как наш сплайн класса С2. В случаях I, II, IV и V фиксированными будут соответственно коэффициенты ао; ао и а\ \ ао, <*1, а2 и аз; ао, ах, а2, аз и 04.)

м 2

Рассмотрим функционал Ф 1(и) = ^ (^(6 + — Ут.1+к) • Будем искать в Р| такой

/с=0

полином д1, который минимизирует функционал Ф1

м

= + к1г^ Ут1+к)^ -> тт(а3,а4,а5,а6,а7)

к=о

и удовлетворяет условиям гладкой склейки:

9т=а10 = д^1(Н), 91(0) = а11=д11(Н), = а12 = (1)

причем 50(0) = а£ = дь^г{Н), д'0(0) = а? = д'ь_х{Н), = а12 = 9ь^ в периодиче-

ском случае.

у”

В непериодическом случае задаются начальные условия 1 : а[] = уо, а® = у'0, а2 = — • Здесь М + 1 - количество точек осреднения, т.е. необходимые точки для построения полинома, д1, а ш+1 - количество точек, входящих в область определения полинома Будем предполагать, что значения заданной функции у к, к = 0,1 известны с некоторой

погрешностью и с уменьшением шага точность измерения будет увеличиваться, то есть если функция / £ С*8[а, Ь] задана на сетке {хк}*к=ъ значениями ук, то |ук — /(хк)| ^ С/г8.

Определение 1. Б-сплайном седьмой степени будем называть функцию 8'1пМ, которая совпадает с полиномом седьмой степени д1 на отрезке [£ь &+!.]•

1Если функция задана таблицей, то у’0, ,..., у^ можно вычислить с помощью формул числен-

ного дифференцирования высокого порядка аппроксимации, например

(г) _ <і(г)ІУта(ж) Уо ~ (1хг

+ 0(1іп+1 г) при г = 1,... ,р,

х=0

где Мп(х)~ интерполяционный полином степени п, построенный по значениям уо, 2/1, ■ ■., уп- В форме Ньютона этот полином имеет вид: Ип{х) = уо + ]С"=і РвІУо, Уі, - ■ ■, Ув)х{х — Н)... (х — (з — 1)/г), где Р*ІУ0,Уі, ■ ■ ■ ,Уа) = Е;=о сЇУз-з/($'■№) - в-я разделенная разность.

При п — 8 мы получаем формулы вида:

У>0 = “ ШЙ[228%0 “ 67202/1 +1176%2 “ 156802/3 + 14700^-

-94082/5 + 39202/6 - 960у7 + 105»8] + 0(Н8),

-г [2953І2/0 - 138 5 282/1 + 3 1 29 84у2 - 4Ш72у3 +

+4353302/4 - 284256г/5 + 120 0 082/е - 29664^/7 + 3267г/8] + 0(И7).

Будем минимизировать функционал Ф1 по коэффициентам аз, 04, 05, щ, а7 (в случаях І, II, IV и V минимизируем по коэффициентам а\,... ,0,7] а2, ■ ■ ■,07; 04,...,07 и а^,...,а7 соответственно).

Для этого продифференцируем Фц по этим коэффициентам и приравняем результат к нулю. Получим:

а

,'3h3S6 + a{h4Sr + al5h5S8 + al6h6S9 + a\h7 Sw = al3h3S7 + al4h4S8 + al5h5S9 + al6hbS10 + al7h7Sn = a3h3Ss + a\h4Sg + al5h5Sw + aleh6Sn -j- al7h7 S12 — <^3 al3h3Sg + a{hASio + al5h5Sn + al6h6S\2 + al7h7 S13 = c\ al3h3S10 + a[h4Sn + al5h5Sl2 + al6h6S13 + al7h7Su = 4

2

J

c3

(2)

где

M M

Si = X! 4 = X [{уті+к - alo - Ahk - al2{hk)2) k2+3

(3)

к=0 к=О

Сделаем следующую замену переменных щ = щ1ъг, г = 0,1,... ,7. Из (2) и (3) получаем систему уравнений для определения коэффициентов при старших степенях:

ajjSs “Ь a2S4 + al2Sb + + a,\S7 + al^Sg + a$Sg + a7S\o — Pj

al0S4 + o,2S^ + al2Se + a3S7 + a\_Sg + a^Sg + ag/Sio + ays'll = al0S5 + al2Se + al2S7 + al3S% + a\Sg + al5Sio + аг65ц + alrSi2

+ a2S7 + a2Sg + сil3Sg + a\Sw + a^Su + a^Si2 + al7Si3 = P\ al0S7 + a2S8 + a2Sg + ag^io + a4Sn + agS'12 + а^із + a7Sii — P5

pi

Ґ2

- Pl

- ]3

(4)

M

к-0

Уравнения гладкой склейки (1) дают нам следующую систему:

-1-і

+ таг1”1 + т^а1^1 + т^а1^1 + rrilal±,~ + т°а l5~1 + rriaal^L + т'а7 al{~1 + 2та121 + Зт2^1 + Ат3а1^~1 + Бга4^-1 + 6 тьа1^1 + 7тьа17-1 = а\

21m5al7~1 = ак

-.2 = 2-1

,з~г-і

4-г—і

-.7^2-1

— ал

-1-1

+ Зта13 1 + 6 mzal4 1 + 10 m6aLb 1 + 15m4ag

,2~г-і

з~г-і

,4=1-1

(5)

В дальнейшем волну над переменной а\ , і = 0,1,..., 7 будем опускать.

Системы уравнений (4) и (5) составляют систему линейных алгебраических уравнений для определения всех коэффициентов полиномов й'-сплайна. Запишем эту систему в матричном виде. Введем следующие обозначения:

S3 • •• S5 56 ...

Aq = , Аг = 7

s7 ... Sg Sio • • ■ ^14

1 т 2 т Я А т т т т6 т7

В0 = 0 1 2т , В1 = 3 т2 4 т3 5т4 6 т5 7 т6

0 0 1 3 т 6 то2 Ют3 15 т4 21т5

Р1 =

и X

21

,Х

21+1

\п;

Е и 0 - единичная и нулевая матрицы.

( 4 N

п1 С* 4

\ а7 )

где I = 0,1,..., L — 1.

(6)

Систему для определения коэффициентов 5-сплайна записываем в виде следующей блочной матрицы, клеточные строки которой состоят по очереди из трех и пяти строк (аналогично столбцы):

-Е 0 0 0 . • Во Вг \ ( Xй \ ( 0 \

Ао Аг 0 0 . . 0 0 X1 р°

В0 Вг —Е 0 . . 0 0 X2 0

0 0 Ао Аг . . 0 0 X3 __ Р1

0 0 В0 Вх . . 0 0 X4 0

0 0 0 0 . • А0 ) V X21-1 / V в1'-' /

Матрицу этой системы обозначим через С? 2 В случае I клетки матрицы в имеют вид:

А0 =

р/

51 52 . • 58

> А\ —

57 58 . ■

Вп = 1

В\ = II т

т

/ Р[

V П

х21 = ( 4 ), х21+1

а\ \

где I = 0,1, ...,£ — 1.

В случае II иммеем:

52 53

Ао — > —

57 58

Р1 =

В случае IV имеем:

( Р[ \

54 ... 5д 5Э ... 514

Х21 = ) ,Х2г+1

В0 =

1 т О 1

, Вх

т ... ш 2т ... 7тв

/ ак

\ а7

где I = 0,1,..., Ь — 1.

54 - • 57 58 . ■ • 5ц

А0 = 57 • • 5ю > А\ — 5ц • 514

Ва =

1 ш ТВ2 т3 л т т5 т6 т7

0 1 2т 3 т2 4 т3 5 т4 6т5 7 т6

0 0 1 3 ш ) Вг — 6 ш2 Ют3 15т4 21 т5

0 0 0 1 4 т О 3 (О 20т3 35т4

( Р\

Х21 =

,х

21+1

, где / = 0,1,...,Ь- 1.

V 7

В случае V:

55 • . 59 5ю • ■ 512

Ао — 57 • • 5ц > — 512 . • 514

(7)

Для непериодического сплайна изменения в системе (7) происходят только в первой строке, которая имеет вид:

Е 0 0 0 ... О О

( Х° \

\ X21-1 /

У0 где ¥° =

У о

Ы

\н2у[

(8)

Преобразуем нечетные строки матрицы таким образом, чтобы избавиться от зависимости от нечетных неизвестных Х1,Х3,..., Х2!'~1. Для этого необходимо, чтобы определитель матрицы А\ был отличен от нуля.

сІеЦАї) =

1

-М5(М - 4)(2М + 1 ){М + 5)(М + 4)2(М - 3)2(М + З)3х

428072077674867225600000 х (М - 2)3(М + 2)4(М - 1)4(М + І)5 62747360064М + 34407503328М2 +

+ 191507803664М3 + 55638343352М4 - 121599386484М5-

- 109726189632М6 + 33394853149М7 + 94469535622М8+

+ 14160714141М9 - 34983098157М10 - 13778627821МП +

+ 3270586487М12 + 2372934151М13 + 146954773М14 - 114362976М15-

- 26623653М16 + 3671136М17 + 3749928М18 + 276444М19 - 323568М20-

-74844М21 + 4788М22 + 3024М23 + 252М24 + 14649189120^ .

Уравнение с1еі(А]) = 0 имеет 6 положительных вещественных корней:

{1; 1.474; 2; 2.471; 3; 4}. Так как М - целое положительное число, то при М ^ 5, существует обратная матрица3 А^1.

Чтобы избавиться от зависимости от нечетных неизвестных Х1,Х3,..., Х2^"1, из третьей строки вычтем вторую, умноженную на матрицу ВхА^1. Выпишем для наглядности первые три строки:

—Е

Ао

0 0 0 .. • В0 в 1

Аг 0 0 .. . 0 0

0 -Е 0 .. . 0 0

( Х° \

Х2Ь-1

/

(9)

Произведем аналогичные преобразования для остальных пар соседних строк с номерами 21 + 1 и 21, где I — 1,... ,Ь — 1, а так же для первой и последней строки. Обозначим через

Вп =

1 т т2 ... т4 т5 т6 т7

0 1 2т ... 4ш3 5 т4 6 ш5 7 т6

, вг = Юто3 15т4 21т5

10т2 20т3 35т4

0 ... 1 5т 15т2 З5т3

Р[

РІ

н

и X21 =

( 4 \

V «4 /

^■21+1

4 . где 1 = 0,1.

3В случаях I, И, IV и V обратная матрица существует при М ^ 7 ; М ^ 6 ; М ^ 4 и М ^ 3

соответственно.

И матрицу Во — В\Ах 1. Преобразованная система (7) расщепляется на системы для нахождения X21 и для нахождения Х21+1, I = 0,1,..., Ь — 1. Система для нахождения X2* имеет вид:

в2а21рь-1 \ -В2А21Рй

( —Е О О и -Е О

О и -Е

\ о о о .

то есть

.. 0 и \ / Х° \ /

.. 0 0 X2

.. 0 0 X4 =

.. и -Е V Х2Ь~2 ) \

= -ВгА^Р1, 1 = 0,1,..

= 0,1,..., Ь - 1:

-В2А21Р1

/ А0 Аг 0 0 ... О 0 \ / Х° \ ( Р° \

О 0 Ао Аг ... О О

\ О О О О ... Ао Аг /

X1 \ X21-1 )

то есть

АгХ21+1 = Р1 - А0Х21, I = О, 1.

(10)

(П)

Для непериодического сплайна система уравнений для определения коэффициентов имеет следующий вид:

Х° = Уо, ИХ21 - Х2^ =-ВгА^Р1, I = 0,1,..., Ь - 2, АгХ2Ш = Р1 - А0Х21, 1 = 0,1,...,Ь-1.

(12)

2. Существование и единственность 5-сплайна седьмой степени

Имеют место следующие теоремы.

Теорема 1. Пусть М ^ 5, тогда, для любой функции /(х), заданной на отрезке [а,Ь] своими значениями у к в точках х^ = а + кН, /г = ■■ , существует единственный непе-

риодический сплайн седьмой степени класса С2.

Доказательство. По формулам (12) последовательно находим Х°, X1,..., Х^-1. Тем самым все коэффициенты полиномов, составляющих сплайн, однозначно определены. □

Теорема 2. Пусть числа т и М ^ 5 таковы, что собственные числа матрицы и не равны корню степени Ь из единицы (Ь-число полиномов, составляющих сплайн). Тогда, для любой периодической функции /(х), заданной на отрезке [а, Ь] своими значениями ук в

точках Хк = а + кк, /г = ——, существует единственный периодический сплайн седьмой

степени класса С2.

Доказательство. Рассмотрим систему (10). Умножим первую строку системы на матрицу и и сложим со второй, полученную вторую строку умножаем на {7 и складываем с третьей и т.д. Поменяем знаки и получим следующую систему:

( Е 0 ... 0 -и

0 Е ... 0 ... -и2

0 0 ... Е -и1

V 0 0 ... 0 ... Е — UL

-1 5 А = т!3: :х0 wb2a21 і

\

X2

_y2(J-l)

\

( Do \ Di

D(i-1)

e-ul j \ x*L-v ) \ d(l_u J

i_1, I = 0,1,...,L-1.

По условию теоремы det(i? — UL) ф 0, следовательно

X2(L-1) = (E_ ul)-1DL-x, X21 = Di + ul+lX2(L~l\

а из системы (11) следует, что X2i+] = A2l{Pi — A\X21). Таким образом, все коэффициенты периодического сплайна найдены. □

Аналогичные утверждения справедливы для случаев I, ИДУ и V при соответствующих М.

Теорема 3. Пусть периодическая функция f(x) Е С8[а,Ь], и пусть выполнено условие \/{хк)~Ук\ ^ C7i8+e, е ^ 0. Пусть, кроме того, числа т, М, р, п таковы, что det(Ai) ф 0 и собственные значения матрицы U по модулю меньше единицы. Тогда периодический сплайн седьмой степени S7mM € С2[а,Ь] с узлами на равномерной сетке имеет дефект 5 и для х G [а, 6] справедливы следующие оценки:

ЙГ

|Г(х) “ й сг/*8"’

для г = 0,1,... ,7; при г = 3,... ,7, х ф в этом случае <р^(£/) = (£; + 0), где tp{x) =

’т,м(

f(X) - SIi.m(X)-

Доказательство производится аналогично доказательству теоремы о сходимости 3 сплайна в работах [1, 2, 3].

Аналогичные теоремы справедливы и для непериодического варианта, а также для случаев I, II, IV и V.

3. Устойчивость 5-сплайна седьмой степени

Для устойчивости й'-сплайпа необходимо, чтобы собственные числа матрицы 17 по модулю были меньше единицы (а если они еще и различны, то и достаточно). Собственные числа матрицы и определяются из уравнения

(М{и - АЕ) = 0. (13)

Для случая малых значений М (при 3 < М < 10) в результате расчета были получены значения собственных чисел матрицы II. В случае р = 0 матрица С/ состоит из одного числа. Показано, что при М = 7 и т = 1,... ,7 матрица II = 0. Некоторые наиболее интересные

полученные значения ш и М, при которых достигаются наименьшие значения тах | Л( и аппроксимация б'-сплайнами устойчива, представлены в таблице.

Как показано в случаях п = 3 и п = 5, для обеспечения этого условия устойчивости необходимо перекрывание. Это означает, что имеются такие элементы исходной таблицы значений функции, которые участвуют в определении коэффициентов не менее двух соседних полиномов, составляющих сплайн. Если перекрывание достаточно большое, то это в ряде случаев является и достаточным условием [1], [3]. На практике наиболее употребительными являются те сплайны, для построения которых используется небольшое число точек осреднения М.

Таблица

Собственные числа матрицы II

р м ш тах | А і | Р М т тах | А* Р М т тах | А*|

0 7 1...7 0 1 7 5 0,107 3 5 3 0,371

0 8 1 0,0006 1 7 6 0,173 3 6 1 0,499

0 8 2 0,0022 1 8 1 0,0931 3 6 2 0,495

0 8 3 0,0044 1 8 2 0,452 3 6 3 0,577

0 8 4 0,005 1 8 3 0,129 3 6 4 0,305

0 8 5 0,00435 1 8 4 0,0341 3 7 1 0,555

0 8 6 0,00218 1 8 5 0,0695 3 7 2 0,623

0 8 7 0,000622 2 5 1 0,570 3 7 3 0,305

0 9 1 0,00263 2 5 2 0,573 3 7 4 0,763

0 9 2 0,00774 2 6 1 0,272 3 7 5 0,655

0 9 3 0,0115 2 6 2 0,141 3 7 6 0,568

0 9 4 0,00749 2 6 3 0,233 4 4 1 0,690

0 9 5 0,00230 2 6 4 0,501 4 4 2 0,881

0 9 6 0,00807 2 7 1 0,324 4 5 1 0,715

0 9 7 0,00625 2 7 2 0,242 4 5 2 0,824

0 9 8 0,00226 2 7 3 0,248 4 6 1 0,756

1 6 1 0,167 2 7 4 0,0908 4 6 3 0,770

1 6 2 0,0667 2 7 5 0,321 4 7 1 0,787

1 6 3 0,0500 2 8 1 0,373 4 7 2 0,693

1 6 4 0,0667 2 8 2 0,387 4 7 3 0,790

1 6 5 0,167 2 8 3 0,131 4 7 4 0,817

1 7 1 0,165 2 8 4 0,271 4 8 1 0,812

1 7 2 0,0253 2 8 5 0,148 4 8 2 0,698

1 7 3 0,0570 3 5 1 0,428 4 8 4 0,765

1 7 4 0,0943 3 5 2 0,306 4 9 2 0,714

Авторы благодарят студента ВМК Кочнева Ю.К., который выполнял вычисление собственных чисел матрицы устойчивости II.

Литература

1. Силаев, Д.А. Приближение Б-сплайнами гладких функций / Д.А. Силаев, Г.И. Якушина // Труды семинара имени И.Г. Петровского. - М., 1984. - Вып. 10. - С. 197.

2. Полулокальные сглаживающие сплайны класса С1 / Д.А. Силаев, А.В. Амилющенко,

А.И. Лукьянов, Д.О. Коротаев // Труды семинара имени И.Г. Петровского. - М., 2007.

- Вып. 26. - С. 347 - 367.

3. Силаев, Д.А. Дважды непрерывно дифференцируемый полулокальный сглаживающий

сплайн / Д.А. Силаев // Вестн. Моск. ун-та. Сер. 1, Математика, механика, 2009. - №5.

- С. И - 19.

Силаев Дмитрий Алексеевич, кандидат физико-математических наук, доцент, кафедра «Общие проблемы управления», Механико-математический факультет, Московский государственный университет, [email protected].

Ингтем Женни Гастоновна, кафедра математической физики, факультет «Вычислительная математика и кибернетика», Московский государственный университет, [email protected].

Поступила в редакцию 24 мая 2010 г.