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

УДК 681.3.513

Е.А. Никулин

ПОЛИНОМИАЛЬНАЯ АППРОКСИМАЦИЯ ЛИНИЯМИ И ПОВЕРХНОСТЯМИ

Нижегородский государственный технический университет им. Р.Е. Алексеева

Тема работы: Построение аппроксимационных линий и поверхностей.

Цель работы: Разработка алгоритмов приближения полиномиальными линиями и поверхностями. Метод решения: Равномерный и взвешенный методы наименьших квадратов. Оригинальность: Использован блочный способ расчета коэффициентов полиномов. Выводы: Получены алгоритмы конструирования аппроксимационных полиномов.

Ключевые слова: полином, интерполяция, аппроксимация, матрица, функционал невязки.

Данная работа продолжает серию статей [1-2], посвященных блочному методу конструирования полиномиальных линий. На сей раз мы будем строить аппроксимационные линии, которые в меру своих изгибательных способностей, присущих степени полинома, проходят максимально близко к заданным узловым точкам, а через некоторые из них даже могут проходить точно. Полученные результаты будут распространены на поверхности.

Необходимость в полиномиальной аппроксимации возникает, по меньшей мере, в двух случаях, реализуемых понижением степени полинома m от ее интерполяционного значения n;

• для уменьшения межузловых осцилляций, присущих интерполяционным полиномам высоких степеней;

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

Вектор коэффициентов полиномиальной функции

p(t)=So +М+...+smtm = STm(t) Vtefo,tn], (1)

совпадающей в отсчетах параметра t0, t1, ..., tn с вершинами узловой сети u0, u, ..., un и записанной с помощью вектора базовых функций Tm(t)=[l t t2 ... tm]т в виде разложения по степеням параметра t, находится решением системы линейных алгебраических уравнений (СЛАУ)

SW=U (2)

с блочной матрицей Вандермонда W=[Tm(t0) Tm(t1) ... Tm(tn)], составленной из m+1 -й степени n+1 -го узлового отсчета параметра t, и блочным вектором узловых точек U=[u o ui .■■ un]. При m=n СЛАУ (2) с невырожденной (n+l)x(n+1)-матрицей W имеет интерполяционное решение

S=UW-1, (3)

обращающее в нуль квадратичный функционал невязки

n

3=2pfe)-2 =(SW-U)(SW-U)T . (4)

i=0

Для уменьшения межузловых осцилляций, присущих полиномам высоких степеней m > 7^8, можно задавать значение m меньше числа интервалов n ценой уменьшения изгибаемости линии и потери ею интерполяционных свойств. Система (2) с числом неизвестных m+1, меньшим числа уравнений n+1, становится переопределенной, а ее МНК-решение [3]

S=UH , где H=Wт (WWт )-1, (5)

доставляет критерию (4) минимальное значение

© Никулин Е.А., 2014.

3т1П =|-и|2 = и(1и+1 -HW)ит >0 .

Получаемая в результате аппроксимационная полиномиальная линия (1) не проходит точно ни через один узел. Это расплата за уменьшение межузловых осцилляций, которыми грешат интерполяционные полиномы высоких степеней при большом числе узлов. Иллюстрацией сказанного служат аппроксимационные полиномы степеней т=1^6, построенные на детерминированной сетке с числом интервалов п=7 (рис. 1, а, б).

При выборет =п решение (5) сводится к (3), сумма (4) достигает абсолютного минимума Зш1д =0, а полиномиальная аппроксимация превращается в интерполяцию — точное прохождение линии через все узлы с максимальными осцилляциями (рис. 1, в). Для преодоления плохой обусловленности матрицы Wт т)-1, при т=п теоретически равной W-1, следует заменить ее расчет прямым обращением матрицы

и

I ✓

■0U5

о

? U5

U

а)

б)

Рис. 1. Полиномиальные линии

в)

Часто координаты вершин узловой сети зашумлены случайными помехами. Допустим, у некоторой эталонной полиномиальной зависимости р э () степени тэ модель измерения

и = рэ (Ь V/ =0,П (6)

содержит вектор случайного смещения \ , с равномерно распределенными элементами гпй(2а)—ае(—а,а). В результате интерполяционная линия р(/) степени т=п, стремясь пройти через все точки и0 ^ип, приобретает хаотические изгибы (рис. 2,а).

Аппроксимационный полином степени т<п обладает замечательной способностью сглаживать осцилляции стохастического характера. По мере увеличения степени линия р(/) проходит все ближе и ближе к узлам и, (рис. 2, б), а значение функционала (4) уменьшается, достигая нуля при т=п, когда сглаживающий эффект исчезает. Наибольший спад зависимости 3(т) при разных параметрах разброса а стабильно наблюдается по достижении степени т=тэ (рис. 2,в). Эту закономерность можно использовать для идентификации неизвестной степени полинома р э (/).

3

А

<>äZö\ —1—1—

—►

123 m 56 78 m

а) б)

Рис. 2. Стохастическая МНК-аппроксимация

в)

4

4

0

Суммирование в (4) квадратичных невязок в равных пропорциях приводит к получению линии, примерно одинаково приближенной ко всем узловым точкам, независимо от ошибок их измерения. Однако во многих практических задачах узловые точки могут быть известны с разными погрешностями. Точные данные несут больше информации о неизвестной зависимости и вызывают больше доверия. Большие ошибки измерений приводят к чрезмерному «вилянию» интерполяционной линии, стремящейся пройти через каждый узел. Алгоритм аппроксимации должен давать решение, проходящее тем ближе к узловой точке, чем точнее она измерена.

Усовершенствуем метод (5) так, чтобы можно было регулировать близость линии (1) к вершинам узловой сети с учетом ошибок их измерения. Для этого введем в функционал (4) весовые коэффициенты ф.>0, входящие в диагональную матрицу Ф=diag^ ф2 ... фп] и обратно пропорциональные дисперсиям погрешностей измерения точек иг :

n

3=^фг. Ipfo. )—и. | =(SW-U^(SW-ü)T ^min . (7)

'=0 S При задании некоторого коэффициента ф^ = 0 вклад невязки |p(^ ) —и J в общую сумму исчезает, что делает линию p(t) полностью независимой от расположения точки uk. Напротив, при фк ^да расстояние между точкой р(^) и узлом ик должно стремиться к нулю.

Выполним минимизацию функционала (7) взвешенным методом наименьших квадратов (ВМНК). Из условия d3/dS=2SW®WT —2UФWT = Om+1 вытекает решение задачи — вектор полиномиальных коэффициентов

S=UH, (8)

где H=ФWT (WOWT )—1. Подставив его в (7), получим минимальное значение функционала невязки

3min = ü(ln+1 —HW)фиT . Очевидно, что функционал (4) и минимизирующее его решение (5) представляют собой частный случай взвешенной задачи (7) с единичной матрицей Ф.

Иллюстрацией работы описанного алгоритма взвешенной аппроксимации служат линии разных степеней на рис. 3, построенные на эталонных узлах рэ (t.) из рис. 2, смещенных

по (6) с разными максимальными отклонениями Go . Точки Uo и U3 сгенерированы при

малых отклонениях g = 0,01, точки выбросов Ui и u5 имеют большое значением g=1 , а

остальные узлы созданы при g = 0,1. Весовые коэффициенты функционала (7) задавались

—2

равными ф. =g. .

Рис. 3. Взвешенная МНК-аппроксимация

При любом выборе конечных весовых коэффициентов ф0 -^фи линия р(/), минимизирующая функционал (7), всегда будет аппроксимационной, поскольку для точного попадания в некоторый узел ик нужно задать недостижимое значение коэффициента ф^ . К то-

му же с ростом отношения фшах/фтщ ухудшается обусловленность обращаемой в (8) матрицы WФWт, что уменьшает вычислительную устойчивость ВМНК.

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

р'(/г-) или ускорения р''(/г-), а вокруг остальных узлов и^ дать возможность линии вести себя

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

Поставим задачу построения приближающего полинома р(/) степени т на сетке узлов {^, и,} V/=0, п, который обладал бы следующими свойствами (рис. 4):

• был закреплен в крайних узлах ио и ип , т. е. являлся интерполяционным решением системы из у=2 уравнений

S[T; ('о ) Т (п )]=[ио ип ] или SWИ = ии . (9)

Можно задать и другие V условий фиксации узловых состояний, скоростей и ускорений путем надлежащего формирования блочных матрицы еЯи вектора ии еЯV. Например, система из v=4 -х уравнений, фиксирующая концевые точки р(/0 )= и0, р(/п )= ип и векторы направлений р'(/0 )=У0, р'(/п )=Уп , имеет вид

S[Tm (/0) Т (п ) Т (/0) Т (хп )]=[и0 ип У0 Уп ];

• проходил в окрестности внутренних узлов и1 , т. е. являлся аппроксимационным решением системы из ц=п _1 приближенного уравнений:

S[T;(н) ... Т;(/п_1 )Ми ... ип_1 ] или SWa«иа. (10)

Число ц и состав аппроксимируемых узлов могут быть изменены путем формирования соответствующих блочных матрицы Wa еЯ (;+1)хц и вектора иа еЯц.

Р( ^р^^ р(^)

Un

U0C

Рис. 4. Приближающая линия

В случае иного, чем в (9) и (10), состава матрицы WH и Wa должны быть полноранговыми и не иметь одинаковых столбцов. Другими словами, условие ST {tj )=u;- должно входить лишь в одну из двух систем, а узел u должен быть либо интерполяционным, либо аппроксимационным, что является вполне естественным требованием. В зависимости от соотношений между числом неизвестных коэффициентов m+1 и размерностями v и ц обе системы могут быть различных типов определенности, поэтому выбор подходящей формулы вычисления вектора S не однозначен.

Алгоритм решения общей задачи условной аппроксимации — минимизации квадратичного функционала невязок

и-1

3 = 2 Р{^ )-Uj | ={SWa -Ua ){SWa-Ua )T ^ min (11)

j=1 S

в подпространстве решений системы 8^ = ии использует операцию псевдообращения матриц любых размерностей и рангов. Применительно к некоторой прямоугольной матрице ЛеЯпхт полного ранга г=шт{п,т} псевдообратная матрица А + еЯтхп вычисляется как

А + =

(АЛ')-'

Л1 (АЛ

-1

при п < т;

Л-1 при п=т; (лтл) Лт при п>т.

Псевдообращение матрицы неполного ранга выполняется методами ортогонального или сингулярного разложений. В частности, нулевая матрица Л=Опхт имеет нулевую псевдообратную матрицу Л + = Отхп .

Как показано в [3], решение линейной задачи наименьших квадратов с ограничениями-равенствами в принятых здесь обозначениях имеет вид

8 = ии^^(иа -ииХ^а Wи+ Wa )+ . (12)

Первое слагаемое 8и = ииWИ" представляет собой нормальное (минимальной длины или нормы) псевдорешение системы 8^ = ии, а второе слагаемое сдвигает вектор 8 и в точку минимума функционала (11).

Путем изъятия из общей задачи одной из систем уравнений получим два частных случая задачи полиномиального приближения:

• в отсутствие аппроксимационной СЛАУ SW ~ Па - задачу интерполяции с решением

8=ии ^ Wи )-1 Wит, которое при ии = и, т=п и Wи = W превращается в (3);

• в отсутствие интерполяционной СЛАУ 8^ = ии задачу аппроксимации с решением

8=иа ^^аГ (^^а WaT )-1,

переходящим в (5) при иа = и и Wa = W.

С помощью перестановок блоков выражение (12) вектора 8 через подмножества узлов и и из (9) и и а из (10) можно записать в типовом виде 8=ПН с регулярным вектором

узлов и=[и0 и1 ... ип]:

К - ^^ Wи+ Wa )+

(»а - Wи Wa )+

Условие единственности решения [3] - равенство ранга составной матрицы Wa ]

числу т+1 коэффициентов приближающего полинома, которое (а это условие существования решения) не должно быть меньшим числа интерполяционных уравнений:

у<т+1<у+^.

Для матриц Wи и Wa из (9) и (10) эти условия сводятся к неравенству 1<т<п . По изложенному алгоритму на рис. 5 построены семейства приближающих полиномов разных степеней т с концами, зафиксированными в крайних узлах. Утолщением выделены интерполяционные линии степеней т=п, обладающие наибольшими осцилляциями.

8 = [Пи Па ]

-^ -^^а )

" 1 ___0__ ]0_1п-11

и 077 0 п-1,\\ 1п-1

0 1 ! о1>И_1 ]

-^ ^ -^ )

ии

V

и

Рис. 5. Примеры приближающих линий

Распространим полученные достижения на трехмерный графический объект — поверхность, описываемую аналогичной (1) бипараметрической векторной функцией

(13)

^ Ш1П . 8

(14)

V ц

рМ=Е&/^ = Тт(08ТЦ(0 ^4о;п], тфо,Тт].

1=0 у=0

Скелетной основой построения поверхности служит матрично организованная узловая сеть и={игу } с индексами вершин г=0,п, у=0,т. Матрица полиномиальных коэффициентов 8, минимизирующая аналогичный (4) функционал невязки

п т , / \ 2

г=0 j=0

находится следующим образом:

• на узловых отсчетах параметров } и {т0,т17.,тт} формируются матрицы Вандермонда:

w=[^) ^) ... ^)]еЯ(^М^1,

П=[ТЦ(Т0) Тц(т1) ... Тц(тт)]еЯ(ц+1)х(т+1);

• при v<п и ц<т эти матрицы не квадратны и аналогично (5) получим аппроксимаци-онное МНК-решение

8=Нт Ш, где Н=Wт (WW т )-1, N=Пт (ППт )-1, (15)

доставляющее функционалу (14) минимум

Зшт = ^^^П-иХ W т8П-и)т );

• в предельном случае v=п и ц=т матрицы W и П становятся квадратными, аппрок-симационное решение (15) превращается в интерполяционное 8=(W-1)T иП-1, а функционал невязки достигает абсолютного минимума 3Ш1П =0.

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

Решение. Возьмем, например, квадратичную поверхность гиперболического параболо-

2 2

ида х - 4у - г = 0 с эталонной параметрической моделью

Рэ (¿,Т)=

12-Т2

*е[- 3,3], те[- 2,2].

(16)

Сформируем узловые отсчеты параметров = / - 3 V/ = 0,6, т/ = / - 2 V/ = 0,4 и сгенерируем вершины узловой сети с координатами

Ру = Рэ(^,Тj)+[гпё(2ст)-ст тпй(2ст)-ст тпй(2а)-а] при а = 0,3. (17)

Построим на одной и той же зашумленной узловой сети набор полиномиальных поверхностей (13) для всех сочетаний степеней полинома vе{l,2,6} и це{1,2,4} (рис. 6). Как и

5

I

Т

4

полагается, билинейная модель с V = ц = 1 представляет собой косую плоскость, поверхности с V = 1 имеют каркасные прямые в ? -направлении, с ц = 1 — в т -направлении, степени V = 2

и (1 = 2 создают параболические изгибы соответствующих каркасных линий.

у=6 3=6.7

3=1.3

3=29.5

|д=4

^щшшш - •

3=2.01

Рис. 6. Аппроксимационные поверхности

Анализируя форму поверхностей, замечаем, что аппроксимация со степенями V = ц = 2 наиболее гладко отражает общие тенденции расположения узловых точек (параболические выгнутости вверх вдоль оси х и вниз вдоль оси г ) и в то же время малочувствительна к их случайным смещениям. К тому же функционал невязки на этой паре {V, ц} претерпевает наиболее резкое уменьшение до 3« 2,3 в сравнении с соседними парами. Дальнейшее увеличение степеней аппроксимационного полинома вплоть до интерполяционных значений v = 6 и ц = 4 изгибает поверхность в локальных местах, приближая их к зашум-ленным узловым точкам и теряя общие закономерности их расположения.

По совокупности признаков заключаем, что эталонная поверхность описывается биквадратичной полиномиальной моделью (13) со степенями vэ =цэ = 2. Численным расчетом была получена матрица векторных коэффициентов этой модели

[_0.028 0.052 _0.082] [0.052 0.066 0.988] [0.053 _0.246 _0.023] S = [1.018 0.044 _0.028] [0.018 0.011 _0.018] [0.012 0.017 _0.012] [_0.002 0.26 0.002] [0.008 0.079 _0.008] [_0.003 0.003 0.003] тогда как эталонной модели (16) соответствует следующая матрица коэффициентов:

"[0 0 0] [0 0 1] [0 _0.25 0]" [1 0 0] [0 0 0] [0 0 0] [0 0.25 0] [0 0 0] [0 0 0]

В системах обработки данных вершины узловой сети и^ часто бывают измерены с разными погрешностями. Точно измеренные точки достоверно представляют неизвестную

Sэ =

поверхность, тогда как сильно зашумленные весьма от нее далеки. В результате суммирование в (14) квадратичных невязок с одинаковыми весами приводит к получению полиномиальной поверхности р(/,т), равномерно приближенной ко всем узловым точкам, независимо от погрешностей их измерения. Хотелось бы, однако, чтобы восстанавливаемая поверхность была тем ближе к узловой точке, чем точнее та измерена, а влияние на нее очевидных грубых выбросов было минимальным.

Усовершенствуем задачу (14) так, чтобы можно было регулировать близость поверхности (13) к вершинам узловой сети с учетом погрешностей их измерения. Для этого введем в квадратичный функционал невязки весовые коэффициенты фу >0, составляющие матрицу

Ф и обратно пропорциональные дисперсиям погрешностей измерения точек и- :

2

\ ( \ Рв-,т/hu

i=0 j=0

^ min.

S

(18)

При задании некоторого коэффициента фы = 0 вклад невязки |р(^ ,Т/)—в общую сумму аннулируется, что делает поверхность полностью независимой от расположения точки и ы .

Выполним минимизацию функционала (18) взвешенным МНК. Из-за произвольности весовых коэффициентов ф- представление взвешенной суммы квадратов невязок с помощью матриц Ф, 8, W и П отсутствует (при всех ф- =1 функционал (14) может быть записан как 1т но существует в поэлементной форме

3=Уф у ту - и у 2 Vi =0, п , у=0т, к=0, V, /=.

к ,/

Из условий дЗ/д$к1 = О Vk,/ составим систему из N=(у+1)(ц+1) линейных уравнений относительно такого же числа элементов матрицы 8:

X8 К^ф ^=Х у^к ^ и у Vк=, Х=.

к,Я.

'J

-v-

akl кЯ

'J

(19)

bt

Для решения этой СЛАУ развернем построчно прямоугольную матрицу 8 размерности ^+1)х(р,+1) в ^-вектор

S =[s00

s,

0ц ¡s10

j — |sv0

s

(20)

Используя принятые в (19) обозначения индексированных коэффициентов Ок/сЯ. и Ь к/: получим векторно-матричную СЛАУ #-го порядка:

SA=B,

где A = A 00 • •• A 0v II £ ak0к0 • • акцк0

_Av0 • •• A vv_ ak0кц • • акцкц

B=[b

00

Ь0ц |b10 ••• Ь1ц ¡•••jbv0

b

Решением этой системы, полностью определенной при ненулевой весовой матрице Ф, является вектор полиномиальных коэффициентов

5=ВА 1,

который по (20) сворачивается в матрицу 8, после чего по (13) строится аппроксимационная поверхность.

Иллюстрацией описанного алгоритма взвешенной аппроксимации служат поверхности на рис. 7, построенные на эталонной узловой сети из Примера, вершины которой смещены по (17) с разными максималь-

ными отклонениями равномерного распределения а- . Точки u00, u

06,

u23 и u46, обозна-

ченные символом О, сгенерированы при малых отклонениях а = 0,05, точка выброса и40 с большим значением а=5 обозначена □, а остальные точки и-- соответствуют а = 0,5 . Весо-

_2

вые коэффициенты функционала (18) задавались равными ф-- =а_ . Анализ формы поверхностей предоставляется читателю.

\'=1 ! \'=2 ! \'=6

у 3=270

Ц=1

U06 - о

'U, >■

^ X

¡1=2

3=196.4 .* о

3=38.6

3=15.7

3= 178.9

о

3=22.8

ц=4

Рис. 7. Взвешенная МНК-аппроксимация Библиографический список

1. Никулин, Е.А. Построение составных линий с различными свойствами сопряжения и произвольными условиями фиксации сегментов // Труды НГТУ. 2006. Т. 58. Вып. 11. С. 5-12.

2. Никулин, Е.А. Прямой параметрический синтез сплайновы х л ини й // Системы обработки информации и управления // Труды НГТУ. 2007. Т. 65. Вып. 14. С. 123-130.

3. Лоусон, Ч. Численное решение задач метода наименьших квадратов / Ч. Лоусон, Р. Хенсон. -М.: Наука, 1986.

Дата поступления в редакцию 22.04.2014

E.A. Nikulin

POLYNOMIAL APPROXIMATION BY LINES AND SURFACES

Nizhny Novgorod state technical university n.a. R.E. Alexeev

Subject: Construction of approximation lines and surfaces.

Purpose: Design of algorithms for polynomial approximation lines and surfaces.

Methodology: Uniform and weighted least squares methods.

Originality: Block method of calculating the coefficients of polynomials is used.

Findings: Algorithms of constructing approximation polynomials are obtained.

Key words: polynom, interpolation, approximation, matrix, the residual functional.