К задаче полиномиального интерполирования с кратными узлами Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2010. Вып. 3

УДК 519.65

А. Ю. Утешев, Г. Ш. Тамасян

К ЗАДАЧЕ ПОЛИНОМИАЛЬНОГО ИНТЕРПОЛИРОВАНИЯ С КРАТНЫМИ УЗЛАМИ*)

Задача. Построить полином у = Н(х), имеющий заданные значения и значения своих производных в узлах интерполяции {Х1,. ,.,х8}с€:

Эта проблема возникла в начале XIX в. как обобщение интерполяционной формулы Лагранжа в задаче разложения рациональной дроби на простейшие. В дальнейшем потребность в ее решении возникала у Эрмита в задаче оценки величины определенного интеграла по заданным значениям подынтегральной функции и ее производных на концах интервала интегрирования [1, 2]. Именно Эрмит поставил задачу поиска полинома Н(х), минимально возможной степени удовлетворяющего таблице значений (1):

и показал единственность такого полинома. Хотя в большинстве современных публикаций [3-6] данный полином называется интерполяционным полиномом Эрмита, но различные формы его записи вызвали к жизни названия «интерполяционный полином Серре» [7] и «интерполяционный полином Лагранжа-Сильвестра» [8]. В XX в. интерес к решению задачи стимулировался теорией матриц - при вычислениях аналитических функций от матриц, имеющих кратные собственные значения.

К сожалению, известное в настоящее время аналитическое представление полинома Эрмита весьма громоздко [4-6]:

Утешев Алексей Юрьевич — профессор кафедры математического моделирования энергетических систем факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 31. Научные направления: символьные (аналитические) алгоритмы для систем полиномиальных уравнений и неравенств, вычислительная геометрия. E-mail: [email protected].

Тамасян Григорий Шаликович — доцент кафедры математической теории моделирования систем управления факультета прикладной математики—процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 17. Научные направления: недиффе-ренцируемая оптимизация, негладкий анализ, вариационное исчисление, теория управления, вычислительная геометрия. E-mail: [email protected].

*) Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 09-01-00360).

© А. Ю. Утешев, Г. Ш. Тамасян, 2010

F (xi ),...,F(ni-i) (xi), F (x2),..., F(n2-i) (x2),

F(xs),..., F(n"-i) (xs).

(1)

deg H < ni + П2 + ...+ ns - 1,

(2)

1=3

и, как правило, используется для практических целей только при небольших значениях индексов: либо при в = 2, либо при и\ = ... = п8 = 2.

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

Ш(х) = П(х - х3 )Щ .

3 = 1

Положим

шз(х) = 7-^Г ПРИ 3 = 1,3.

(х — Ху )пз

Утверждение 1. Существует набор полиномов {и1(х),... ,и8(х)} С С[х], удовлетворяющих тождеству

и!(х)Ш1(х)+и2(х)Ш2(х) + ... + и8(х)Ш8 (х) = 1. (4)

При условии (\egUj < щ {] = 1, в) такой набор полиномов будет единственным; в этом случае полином иу (х) совпадает е отрезком ряда Тейлора разложения дроби 1

. . по степеням х — ж,-. Шз (х) 3

Доказательство. Воспользуемся теоремой о разложении правильной рациональной дроби на сумму простейших [9, 10]: в качестве искомых полиномов можно взять числители разложения дроби

1 _ М1(ж) П2{х) и8{х)

Ш (х) (х — х1)П1 (х — х2 )П2 (х — х8)

на сумму простейших. Из тождества (4) выразим и1(х):

1 и2(х)Ш2(х) +... + и8(х)Ш8(х)

ив

и1(х) =

Ш1(х) Ш1(х)

Все полиномы Ш2(х),... ,Ш8(х) имеют множителем (х — х1 )П1, поэтому при разложении правой части последнего тождества в ряд Тейлора по степеням (х — х1) члены

степеней, меньших п\, образуются только при разложении дроби -——. Таким обра-

Ш^х)

зом, при ограничении deg и1(х) < П1 полином и1 однозначно определяется отрезком разложения функции ^ в ряд Тейлора по степеням {х — х{)ш, в этом разложении

рассматриваются члены степеней до (п1 — 1)-й включительно.

Для остальных полиномов иу(х) доказательство проводится аналогично. I

Х = Хо

Обозначим при j = 1, s

fiix) = F(Xj) + x - Xj) + - Xjf + ...

F(nj -i)(T ■)

+ -*>>•'-'■ (5)

Утверждение 2. Полином

G(x) = fi(x)ui(x)Wi(x) + f2(x)u2(x)W2(x) + ... + fs(x)us (x)Ws (x) (6)

является обобщенным интерполяционным полиномом.

Доказательство. Покажем, что величины полинома G(x) и его производных в точке xi определяются первой строчкой таблицы значений (1). Из (4) выразим ui(x)Wi(x) и подставим в (6):

G(x) = fi(x) + [f2(x) — fi(x)]u2 (x)W2 (x) + ... + [fs(x) — fi(x)]us(x)Ws(x). (7) Все полиномы W2(x),..., Ws(x) имеют множителем (x — xi)ni, поэтому

G(xi) = fi(xi), G'(xi) = f (xi),..., G(ni-i)(xi) = f[ni-i)(xi),

и по способу определения полинома f i (x) полученные значения совпадают с табличными.

Доказательство для остальных узлов интерполяции проводится по аналогии. ■ Замечание! Имеет место следующая оценка степени обобщенного интерполяционного полинома:

deg G^ni + ... + ns — 1 + max щ. — 1,

k=l,s

т. е. его степень, как правило, выше степени (2) полинома Эрмита.

Получим теперь явное представление для полиномов Uj{x) (j = l,s) из тождества (4).

Утверждение 3. Имеет место представление uj (x) в виде

jy jJ 0<ki + ... + ks-kj <nj-i i=i v j l

£=j

Здесь CM = N!/[M!(N — M)!] обозначает биномиальный коэффициент.

Доказательство. Воспользуемся второй частью утверждения 1. Найдем, к примеру, явное представление для ui(x). С этой целью разложим дробь

1 _ 1

Wi{x) (х - х2)П2 . .. (х - xs)n°

в ряд Тейлора по степеням x — xi, ограничившись членами степеней, меньшими ni.

Для каждой дроби 1/(х — ж^)™3 (j = 2, s) выписываем формулу биномиального разложения

1 = 1 у^ се f xt-xV

(х - xj)ni (xi - xj)ni j-* гч+£-1 _xjj

и перемножаем получившиеся произведения. Тогда 78

8 П1-1 , N к1

1 1 ТТ тг гк* I Х1~Х ) ,

11 I _ _ „ I + • • • ■

(х) ^1(Х1) VХ1 - х,

а поскольку deg и1 < П1, получаем

и1(г)=й4л 5: (8)

0<к2 + ... + к„<П1-1 ,=2 у 1 '

Из (8) в случае произвольного ] имеем

и- (х)

= шТТл при = I] (9)

к=0

и»= 53 (Ю)

к1 +... + ка-к; =к£=1v , 01

Очевидно, что для всех ] = $

и-о = 1 • (11)

■

Следующее утверждение дает альтернативное представление для коэффициентов и0к.

Утверждение 4. Обозначим

в

е=з

для j = 1, в, то = — 1. Тогда справедлива такая рекурсивная формула вычисления и-к при к ^ 1:

1к

т=1

в ней, согласно (11), считаем и-0 = 1.

Доказательство проиллюстрируем для случая ] = 1. На основании определения (9) полинома и1(х) получаем его в результате разложения дроби Ш1 (х1)/^.(х) по степеням х — Х1; в этом разложении оставляем члены степеней до (п1 — 1)-й включительно. Тогда в разложении произведения ^.(х)и1(х)/^1 (Х1) по степеням х — Х1 коэффициенты мономов (х — Х1), (х — Х1)2,...,(х — Х1 )П1 -1 должны обращаться в нуль при свободном члене, равном 1. Разложим полином Ш1(х)/Ш1(х1) по возрастающим степеням х — х1:

(х)

„г , \ = 1 + Ш1{х - X!) + и2{х - Х1)2 + ... . (14)

1^1 (х1)

Получаем систему линейных уравнений

и

ип + ш1и10 = 0,

и 12 + + ^2^10 = 0, (15)

и1к + Ш1и1,к~1 + ...+ ^к ию = 0.

Выразим теперь коэффициенты ш1,...,шк через х1,...,ха. Имеем

(х) = (х — Х1 — (х2 — Х1))П2 . .. (х — Х1 — (х8 — Х1))Пв ,

и, с одной стороны, искомые коэффициенты получаются по формулам Виета. С другой стороны, формулы Виета не очень удобны в виду того обстоятельства, что интересуют коэффициенты при возрастающих степенях х — х\. Потому воспользуемся другим подходом для их определения - через суммы Ньютона [8, 9]. Если ввести новую переменную X = х — х1, то корнями полинома

(X — (х2 — х1 ))п2 ...(X — (хя — х1)Уп

будут числа х2 — х1,...,х8 — х1 в соответствующих кратностях. Выражения

П2(х2 — х1)т + ... + п8(х8 — х1)т

при натуральных т называются суммами Ньютона этого полинома. Имеются рекурсивные формулы (Ньютона), связывающие суммы Ньютона с коэффициентами полинома. В смысле такого определения, становится очевидным содержание введенного в условии теоремы обозначения (12): выражение

п2 пч

+ ...+

(х2 — х1)т (х8 — х1 )т

представляет собой сумму Ньютона для полинома

у_ 1 Г...1У-

(х2 — х1^у ' " V (х8 — х1)/ Для получения выражений коэффициентов из (14) представим ^1(х)/^1 (х1) в виде W1(x) 1

П2

у__/у__при У = _

Wl(xl) Уп2 +'"+па У х2 — х1 у \ х8 — х1 у х — х1

и для полинома по У из правой части последнего тождества выписываем формулы Ньютона

#11 + Ш = 0,

Б12 + ^1#11 + 2^2 = 0, (1б)

#1к + + ... + ^к-1#11 + к^к = 0.

Итак, системы (15) и (16) позволяют вывести зависимость коэффициентов и^ от величин (12) - осталось только избавиться от вспомогательных величин шц. Сделаем это

п

следующим приемом: дополним уравнения системы (16) последним уравнением системы (15) и перепишем получившиеся соотношения в матричном виде

/ 5-11 1 0 0 .. 0 0 1 \ 0

5*12 5ц 2 0 .. 0 0 0

^13 512 5ц 3 .. 0 0 Ш2 0

51к 5^-1 51,к-2 51 к-3 . . 5и к ик- 1 0

\ И1к Ц ,к-1 И1,к-2 И1,к-3 . . и и и 10 ) V ик / 0

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

5ц 1 0 0 . .. 0 0

512 5ц 2 0 .. 0 0

513 512 5ц 3. .. 0 0

51 к 51,к-1 51,к-2 51,к-3 . . 5и к

Ц1к И1к-1 И1к-2 И1к-з . .. и11 и 10

0 .

(17)

(к+1)х(к+1)

Разложением определителя по последней строке можно получить линейную зависимость и1к от и1}к-1, и1,к-2, ,...,и10. Однако, для того чтобы показать ее эквивалентность соотношению (13), придется произвести некоторые дополнительные изыскания.

Будем доказывать справедливость представления (13) индукцией по к. Проверка справедливости этой формулы для к =1 и к = 2 очевидна. Предположим, что эта формула справедлива при всех значениях индекса от 1 до к — 1. Преобразуем определитель (17), вычтем из его последней строки первую, умноженную на Ц^к^/к, вторую, умноженную на И ,к-2/к, и т. д., к-тую, умноженную на Ию/к. Элементы последней строки получившегося определителя оказываются равными

и^к — (511^1,й-1 + 5*12^1,к-2 + + • • • + ¿ПяДю),

к 1 1

—^—и 1^-1 — —{Бци 1^-2 + ¿>12 к-3 + • • • + к 2 1

—^—— -^(Зци^к-з + • • • + й^-аСЫ

^п-^ьО.

На основании индукционного предположения все выражения, кроме первого, обращаются в нуль. Раскладывая полученный определитель по последней строке, установим эквивалентность равенства (17) равенству

к!

U\k — ^(¿>11^1,fc-i + Sí2Uítk-2 + SisUi^-з + ... + Slku10)

0,

откуда и следует (13). ■

Формула (13) гораздо более удобна для вычислений в сравнении с формулой (10), которая содержит большое количество громоздких слагаемых с биномиальными коэффициентами.

Пример 1. Имеем

из1 = ^'Ь = Т^5?! + из3 = б^1 + +

Полученное в утверждении 2 выражение для обобщенного интерполяционного полинома О(х) можно использовать для выведения представления интерполяционного полинома Эрмита Н(х). Действительно, как отмечалось выше, Н(х) является остатком от деления полинома О(х) на Ш(х). Применяя представление (6), получаем, что этот остаток равен сумме остатков от деления полиномов /у (х)иу (x)Wj (х) на Ш(х), каждый из которых, в свою очередь, представйм в виде г у (х)Шу (х), где г у (х) означает остаток от деления /у (х)иу (х) на (х — ху )п^, т. е.

Н (х) = т1(х)Ш1(х) + ... + т3(х)Ш3 (х).

Имея явные представления (5) для полинома /у (х) и (9) для полинома иу (х), приходим к следующему результату.

Утверждение 5. Интерполяционный полином Эрмита вычисляется по формуле

у=1 ■" -1

Ее можно переписать в следующем виде:

j=l Jy ■" k=0 1=0

j=l k=0 1=0 Jy Это представление эквивалентно формуле (3), поскольку

1 I dk 1 | Ujk

1 (х — хе)п I Ш (ху)'

„ J

Однако для приложений более удобен другой вариант формулы (18). Введем в рассмотрение следующую операцию. Для представленного по формуле Тейлора полинома /(х) = а0 + а1(х — х*) + а2(х — х*)2 + . .. + ап(х — х*)п обозначим

\п—к

/1 ](х) = ао + а\(х - х*) + а2(х - х*) + ... + ап-к(х - х*)

при к = 1,п и будем считать f^(x) = f(x).

Следствие 1. Интерполяционный полином Эрмита можно представить в виде

(х) р«х)

5 = 1

к=0

к\

П (х)(х - х 5 )к .

(19)

Пример 2. Легко заметить, что при П1 = = ... = п8 = 1 полином Н(х) становится интерполяционным полиномом в форме Лагранжа. При П1 = П2 = ... = п8 = 2 получаем

Н (х) = ]Г

/

Е

5 = 1

Р(х3 ) + (х - х3 ) (

Р '(хд ) + 2Р (хд )]Г

1

х£ - х^

£=з

Р (хз)

1

1 + 2(х - хЛ) ]Г — + ^{х^х - хЛ)

\

"1 х£ - хз

Пример 3. Найти интерполяционный полином Эрмита по следующей таблице:

х3 -1 0 1 2

16 7 8 217

-1 -4 1375

6 -44

-126

Для этого примера имеем: в = 4,П1 = 1,П2 = 3,пз = 4,П4 = 2. Для нахождения полиномов и (х) по формуле (9) потребуются выражения для их коэффициентов П^к из формул (13). Приведем подробные вычисления коэффициентов полинома Пз(х):

#3,

1 3 2 3

х1 - х3 х2 - х3 х4 - х3 -2 -1 1 2

П1 п2 П4

#3,2 = #3,3 =

П1

+

П2

+

П4

21

(х1 - х3)2 (х2 - х3)2 (х4 - х3)2 4'

П1 П2 П4 9

+

+ ■

(х1 - х3)3 (х2 - х3)3 (х4 - х3)3

Тогда

П31

31

1; изА = 5з,1 = ^3,2 = ^(#3,2 + ^Г3,1#3д)

1 39

^3,3 = ^(#3,3 + #3,2^3,1 + #3,1^3,2) = — "д"

15

Т'

и, следовательно,

15,

39

и3(х) = 1 - -{х - 1) + —{х - 1)^ - —{х -

1

Аналогично находим остальные полиномы:

43 35

и Ах) = 1, и2(х) = 1 + 4ж + — ж2, £/4(ж) = 1 - — (ж - 2).

4 6

Обобщенный интерполяционный полином получаем по формуле (6) при Ш1(х) = х3(х — 1)4(х — 2)2, Ш2 (х) = (х + 1) (х — 1)4(х — 2)2, Ш3(х) = (х + 1)х3(х — 2)2, Ш4(х) = (х + 1)х3(х — 1)4, /1(х) = 16,/2(х) = 7 — х + 3х2, /3(х) = 11 — 23х + 41 х2 — 21 х3, /4(х) = —2533+ 1375х. Из (6) имеем

ч 819 12 3507 10522 1П 24709 9 69505 8 27775 7

П(х) = -ж--ж Ч--ж--ж--ж Ч--ж —

к ' 16 8 9 144 16 4

37493 6 244079 5 168767 4 45065 3 2

24~Х 48 Х + 36 Ж 36+ Ж ~Х+ ' и полином Эрмита равен остатку от деления полученного полинома на

Ш (х) = (х + 1)х3(х — 1)4(х — 2)2.

Проиллюстрируем сейчас алгоритм непосредственного построения полинома Н(х) -т. е. без промежуточного нахождения О(х). В самом деле, для полинома Эрмита имеется представление (19):

W3(x)

F "(0)

U2(x)F( 0) + (1 +4ж)^'(0)ж+ -^-х1

+

+

W3(1)

3 15

U3(x)F(l) + (1 - -(ж - 1) + —(ж - 1)2)^'(1)(ж - 1) +

+

2 7 2 4 ' 6

= 2x9 - 3x8 - 4x5 + 5x4 - x3 + 3x2 - ж + 7.

Еще одно упрощение представления полинома Эрмита основано на формуле (7): она позволяет сэкономить на вычислении полинома ui(x). Если в этой формуле в качестве первого узла интерполяции взять тот, что соответствует min п^, и переразложить fi(x)

k=l,s

по степеням ж — Xj:

П1 — 1

\к

f1(x) akj (Ж - Xj )к

k=0

для j G 2, s, то аналогом формулы (19) будет

здесь полагаем а ¡у = 0 для всех к = п\, щ — 1 при п^ > п\. 84

Пример 4. При s = 2 и ni = П2 = n имеем

*<*> = +ш^тг ё № -«) <* - -

к=0 4 у т=0 4 '

здесь через ао,..., ап-1 обозначены коэффициенты разложения полинома /.(х) по степеням X — Х2.

Пример 5. В примере 3 имеем П1 = 1, /.(х) = 16 ив использовавшихся обозначениях находим

F" (0)

U2(x)(F(0) - 16) + (1 +4x)F'(0)x +-—х^

+

+ [Ы*){Р{1) - 16) + (1 - |(х - 1) + Ц{х - l?)F'{l){x - 1) +

+ (1 - |(х - - I)2 + ^(х - I)3] +

+ №)№) - 16) + F'(2)(х - 2)] =

= 2x9 - 3x8 - 4x5 + 5x4 - x3 + 3x2 - x + 7.

Литература

1. Hermite Ch. Sur la formule d'interpolation de Lagrange // J. f. d. reine u. angew. Mathematik. 1878. Vol. 84. P. 70-79 (см. также: Hermite Oh. "Oeuvres" Gauthier-Villars. Paris, 1912. Vol. 3. P. 432-443).

2. Крылов В. И. Приближенное вычисление интегралов. М.: Наука, 1967. 500 с.

3. Mejiering E. A chronology of interpolation: from ancient astronomy to modern signal and image processing // Proc. of the IEEE. 2002. Vol. 90, N 3. P. 319-342.

4. Калиткин Н. Н. Численные методы. М.: Наука, 1978. 512 с.

5. Березин И. O., Жидков Н. П. Методы вычислений: в 2 т. М.: Наука, Гл. ред. физ-мат. лит-ры, 1959. Т. 1. 464 c.

6. Самарский А. А., Гулин А. В. Численные методы: учеб. пособие для вузов. М.: Наука, Гл. ред. физ-мат. лит-ры, 1989. 432 с.

7. Оерре И. А. Курс дифференциального и интегрального исчислений. СПб.; М.: Изд-во Вольф, 1883. Т. 1. 573 c.

8. Гантмахер Ф. Р. Теория матриц. М.: Физматлит, 2004. 560 с.

9. Утешев А. Ю., Калинина Е. А. Лекции по высшей алгебре. Ч. II: учеб. пособие. СПб.: Соло, 2007. 279 c.

10. Фаддеев Д. К. Лекции по алгебре: учеб. пособие для вузов. М.: Наука, Гл. ред. физ.-мат. лит-ры, 1984. 416 с.

Статья рекомендована к печати проф. В. Ф. Демьяновым. Статья принята к печати 1 апреля 2010 г.