239287

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РФ ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ «ВОРОНЕЖСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ»

ИНТЕРПОЛЯЦИЯ АЛГЕБРАИЧЕСКИМИ МНОГОЧЛЕНАМИ.

СПЛАЙН-ИНТЕРПОЛЯЦИЯ

Учебно-методическое пособие для лекционных занятий в вузах

Составители: В.П. Трофимов, А.П. Карпова, М.Н. Небольсина

Издательско-полиграфический центр Воронежского государственного университета

2012

1

Утверждено научно-методическим советом математического факультета 30 ноября 2011 г., протокол № 0500-09

Рецензент канд. физ.-мат. наук, доц. С.А. Шабров

Учебно-методическое пособие подготовлено на кафедре математического моделирования математического факультета Воронежского государственного университета

Рекомендуется для студентов 4-го курса дневного отделения и 5-го курса вечернего отделения математического факультета.

Для направлений: 010100 – Математика, 010200 – Математика. Прикладная математика; специальности 010101 – Математика

2

ТЕОРИЯ ПРИБЛИЖЕНИЯ ФУНКЦИЙ ОДНОЙ ВЕЩЕСТВЕННОЙ ПЕРЕМЕННОЙ

1. Интерполяция алгебраическими многочленами

1. 1. Постановка задачи интерполяции

Пусть для функции f : X → R, X R известны ее значения в (n + 1)-й

точках xi X, i =0,…, n. Запишем эти значения функции f									в табл. 1.1
									Таблица 1.1
x		x 0		x 1			…		x n
f ( x )		f ( x 0 )		f ( x 1 )			…		f ( x n )
Далее будем считать, что выполнено условие
		a ≤ x 0 < x1 <			< x n ≤ b .
Задача приближенного вычисления для заданной табл. 1.1 значения
функции f ( x )		при x ≠ xi ,	i = 0,…, n		называется задачей интерполя-

ции (распространения внутрь).

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

Pn ( x ; x 0 , x 1 , … ,			x n ;	f ) = Pn ( x ; f ),	(1.1)
принимающийвточках x 0 ,	x1 ,	, x n	тежезначения, чтоифункция f		:
f ( x i )	= Pn	( x i ;	f ),	i = 0 , 1, … , n .	(1.2)

Интерполяционным многочленом (интерполянтой) для табл. 1.1

называется многочлен (1.1) степени не выше				n , удовлетворяющий усло-
вию (1.2). Точки x 0 , x1 ,	, x n называются узлами интерполяции.
Вычислениезначения	f ( x ) при x	≠	x i	, i = 0 ,… , n поформуле
f ( x ) ≈ Pn ( x ;			f )	(1.3)
называется интерполяцией функции f		с помощью алгебраического мно-
гочлена.	x [a ; b ],	то вычисление f ( x ) с помощью
Замечание 1.1. Если

(1.3) называют экстраполяцией.

Замечание 1.2. Существуют различные формы записи интерполяционного многочлена.

Теорема 1.1. Для табл. 1.1 интерполяционный многочлен существует и единственен.

Доказательство. Запишем интерполяционный многочлен в виде

Pn ( x ; x 0 , x1 ,… , x n ; f ) = a 0 + a1 x +

+ a n x n .

3

Из условия (1.2) получим систему линейных алгебраических уравнений относительнокоэффициентовинтерполяционногомногочлена a 0 , a1 , , a n

a0a0

a0

+ a x	0		+ a	x2	+ + a		n		xn =		f (x		),
1			2	0					0			0
+ a x			+ a x2		+ + a			xn =			f (x ),
1	1		2	1			n		1			1	(1.4)
+ a x		n	+a	x2	+ + a	n		xn		=	f (x	n	).
1			2	n					n

Так как a ≤ x 0 < x1 < < x n ≤ b , то определитель системы (1.4) (определитель Вандермонда) отличен от нуля:

1 x0 x02 …x0n

=1 x1 x12 …x1n =

1xn xn2 …xnn

=(xn − xn−1 )…(xn − x0 )(xn−1 − xn−2 )…(xn−1 − x0 )…(x1 − x0 ) ≠ 0.

Следовательно, для любой табл. 1.1 можно построить единственный интерполяционный многочлен (его коэффициенты определяются единственным решением системы (1.4).

Теорема доказана.

Важное замечание 1.1. Наиболее общим способом решения задачи интерполяции для табл. 1.1 является линейное интерполирование: интерполянт разыскивается в виде обобщенного полинома – линейной комбинации за-

данной системы линейно независимых базисных функций {ϕ j ( x)}nj =0 :

Φ n ( x ; a 0 , a 1 , … , a n ; f ) = Φ n ( x ; f ) = ∑n

a j ϕ j ( x ) ,

j = 0

где коэффициенты a 0 , a 1 ,

, a n

определяются из условия интерполяции

f ( x i ) = Φ n

( x ;

f ) = ∑n

a j f ( x j ),

j = 0 ,1, … , n .

j = 0

В этом случае для определения коэффициентов

a 0 , a 1 , , a n по-

лучаем аналог системы (1.4)

a

ϕ

0

(x

0

) + a ϕ (x

0

) + + a

n

ϕ

n

(x

0

) = f

(x

0

),

0

1 1

a0ϕ0 (x1 ) + a1ϕ1 (x1 ) + + anϕn (x1 ) = f (x1 ),

a

ϕ

0

(x

n

) + a ϕ (x

n

) + + a

ϕ

n

(x

n

) = f

(x

n

).

0

1 1

n

4

В силу линейной независимости базисных функций {ϕ j ( x)}nj =0 оп-

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

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

ϕ j ( x) = x j	,	j = 0 ,1, … , n .
Для периодической функции	f (t )		с периодом Τ интерполянт разы-
скивается по системе базисных тригонометрических функций
ϕ j (t ) = a j cos		jπt	+ b j sin	jπt	.
		Τ
				Τ

Такая интерполяция называется тригонометрической. В этом случае интерполянт – тригонометрический полином степени N = 2 n + 1 имеет вид

	n			jπt			jπt
	Φ N (t; f ) = a 0 + ∑	a j	cos		+ b j	sin		.
				Τ			Τ
	j =1

Предложение 1.1. Пусть Ρ ( n ) – пространство многочленов степени не выше n . Если f Ρ ( n ) , то для любой табл. 1.1

Pn ( x ; x 0 , x1 ,… , x n ; f ) ≡ f ( x ).

Задание. Докажите предложение 1.1.

Замечание 1.3. Решение системы (1.4) является достаточно сложной вычислительной задачей. Поэтому способом построения интерполяционного многочлена, установленным при доказательстве теоремы 1.1, на практике обычно не пользуются. Более удобный способ построения интерполяционного многочлена был предложен Лагранжем.

1.2. Интерполяционный многочлен Лагранжа

Для каждого узла интерполяции						xk , k = 0,… , n найдем много-
член lk( n ) ( x ) степени n ,	равный нулю во всех узлах, кроме k -го, в кото-
ром он равен 1:		0 ,	i	≠ k
l k( n ) ( x i ) =					,	i , k = 0 ,1… , n .
		1,	i	= k

Систему многочленов l k( n ) ( x ) принято называть базисом Лагранжа.

Нетрудно построить полиномы l k( n ) ( x ). Действительно, зная корни

полинома, имеем:

1) многочлен (x − x0 )(x − x1)…(x − xk−1)(x − xk+1 )…(x − xn ) равен нулю во всех узлах интерполяции, кроме x k ;

5

2) многочлен

(x − x0 )(x − x1 )…(x − xk −1 )(x − xk +1 )…(x − xn )

(xk − x0 )(xk − x1 )…(xk − xk −1 )(xk − xk +1 )…(xk − xn )

равен 1 при x = x k и равен нулю во всех остальных узлах.

Итак, для любого k = 0,… , n

l k( n ) ( x ) =

(x − x0 )(x − x1 )…(x − xk −1 )(x − xk +1 )…(x − xn )

.

(xk − x0 )(xk − x1 )…(xk − xk −1 )(xk − xk +1 )…(xk − xn )

Введем многочлен ω(x) (n + 1)-й степени,

построенный по узлам

интерполяции:

ω(x) = (x − x0 )(x − x1 )…(x − xn−1 )(x − xn ) .

Заметим, что производная многочлена ω(x) в точке x k

′

ω(xk ) = (xk − x0 )(xk − x1 )…(xk − xk−1 )(xk − xk+1 )…(xk − xn ) .

Теперь многочлен l k( n )

( x ) можно записать в виде

l k( n ) ( x )

=

ω(x)

.

(1.5)

(x

− xk ) ω′(xk )

По построению многочлен f ( x k ) l k( n )

( x ) имеет степень n ,

при-

нимает в узле x k

значение

f ( x k

) и равен нулю во всех остальных узлах

интерполяции.

Следовательно, многочлен

(n)

ω(x)

n

n

Ln (x; f ) = ∑f (xk )

lk (x) =∑f (xk )

(1.6)

′

k=0

k=0

(x − xk ) ω(xk )

является интерполяционным многочленом для табл. 1.1 (имеет степень не выше n и Ln (xi ; f ) = f (xi ) , i = 0 ,… , n ).

Формулу (1.6) называют интерполяционной формулой Лагранжа, а полином Ln (x; f ) – интерполяционным многочленом Лагранжа.

Замечание 1.4. Число арифметических операций, необходимых для

вычисления по формуле (1.6), имеет порядок O (n 2 ).

Пример. Найдем интерполяционный многочлен Лагранжа для n = 1 . В этом случае формула (1.6) примет вид

L (x ; f ) = f (x		)	x −x1		+ f (x )	x −x0		=	f (x0 )(x −x1 ) − f (x1 )(x −x0 )		.
1	0		x	−x	1	x	−x		x −x
		0		1		1	0	0		1
Графиком функции L1 (x; f )						является прямая, проходящая через точ-
ки (x0 ; f (x0 )) и	(x1; f (x1 )). Такая полиномиальная интерполяция называ-

ется линейной полиномиальной (не путать с линейным интерполированием (см. важное замечание 1.1)).

6

Задание. НайдитеинтерполяционныемногочленыЛагранжадля n = 2, 3. Замечание 1.5. Поскольку интерполяционный многочлен (1.6) линейно зависит от значений функции f (xi ) , то интерполяционный многочлен для

суммы функций равен сумме интерполяционных многочленов слагаемых.

1.3. Погрешность интерполяции

Погрешностью интерполяции называется разность

rn ( x ; f ) = f ( x ) − Pn ( x ; f ).	(1.7)
Очевидно, что в узлах интерполяции x 0 , x 1 ,	, x n

rn ( x i ; f ) = f ( x i ) − Pn ( x i ; f ) = 0 .

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

Замечание 1.6. Из предложения 1.1 следует, что погрешность интер-

	поляции r ( x i ; f ) ≡ 0			для любой функции f					Ρ ( n ) , где Ρ ( n )		– про-
	странство многочленов степени не выше n .										n + 1
	Найдем погрешность интерполяции для многочлена степени
	( f Ρ( n +1) ,	deg f = n +1).
	В этом	случае	rn ( x; f ) = f ( x) − Pn ( x; f ) есть многочлен степени
	n + 1 и узлы интерполяции x k ,					k	= 0 , … , n		являются его корнями.
	Следовательно,
	rn (x; f ) = c ω(x) = c (x − x0 )(x − x1 )…(x − xn−1 )(x − xn ) ,										(1.8)
	где c = const .
	Продифференцировав по x это равенство n + 1 раз, получим
	rn(n+1) (x; f ) = f (n+1) (x) − Pn(n+1) (x; f ) = f (n+1) (x) = c (n +1)!,
	так как P (x; f ) – многочлен степени n , то P(n+1)								(x; f ) ≡ 0.
	n							n
	Отсюда найдем c =				f (n+1) (x)	. Таким образом, для f Ρ( n+1) погреш-
					(n +1)!
	ность интерполяции имеет вид							f ( n +1 ) ( x )
	rn ( x ; f )		= c ω ( x )				=			ω ( x ).	(1.9)
								( n + 1)!
	Однако для произвольной функции, заданной только табл. 1.1, ничего
	конкретного сказать о погрешности интерполяции нельзя.
	Если функция f C ( n +1) [a ; b						] ( C ( n + 1 ) [a ; b ] – пространство функ-
	ций, n + 1 раз непрерывно дифференцируемых на отрезке [a ; b ]),										то для

погрешности интерполяции можно получить формулу, аналогичную (1.8).

Теорема 1.2. Если f C ( n +1) [a; b ], то для любого x [a ; b ] погрешность интерполяции определяется формулой

7

	rn ( x ; f ) =	f		( n + 1 ) ( ξ )	ω ( x ),	(1.10)
			( n + 1)!
где ξ – некоторая точка отрезка [a ; b ] ( ξ = ξ( x ) [a ; b ]).
Доказательство. Будем разыскивать погрешность интерполяции в ви-
де (1.8), положив с = с( x ),
	rn ( x ; f ) = c ( x ) ω( x ).
Зафиксируем произвольное x [a ; b ], x ≠ xi , i = 0,…, n и рас-
смотримвспомогательную функцию ϕ отпеременной z:
	ϕ ( z ) = rn ( z ; f ) − c ( x ) ω ( z ) =
	= f ( z ) − Pn ( z ; f ) − c ( x ) ω ( z ).
Очевидно, что ϕ C ( n +1 ) [a ; b ] и обращается в нуль в n + 2 точках
отрезка	[a ; b ]: z = x , x 0 , x 1 ,			, x n . По теореме Ролля		функция
ϕ ′ (производная от функцииϕ по				z ) обращается в нуль по крайней мере в

n + 1 точках отрезка [a ; b ], при этом ϕ′ C ( n ) [a ; b ], функцияϕ ′′ равна

нулю по крайней мере в n точках этого отрезка, ϕ ′′ C ( n −1 ) [a ; b ] и так далее.

Таким образом, ϕ( n +1) (z) (ϕ ( n +1) C [a;b]) обращается в нуль по

крайней мере в одной точке ξ [a ; b ] и ξ = ξ( x ). Учитывая, что для любого x

Pn( n +1 ) ( x ; f ) ≡ 0 и ω ( n + 1 ) ( x ) ≡ ( n + 1)! ,

получаем

ϕ ( n + 1 ) ( ξ ) = f ( n + 1 ) ( ξ ) − c ( x ) ( n + 1)! = 0 .

Следовательно,

	c ( x ) =			f	( n +1 ) ( ξ )
и				( n + 1)!
			f ( n + 1 ) ( ξ )
rn	( x ; f ) =						ω ( x ),
			( n		+ 1)!

где ξ = ξ ( x ) [a; b ].

Теорема доказана.

Важное замечание 1.2. Из формулы (1.10) следует, что погрешность интерполяции зависит от выбора узлов интерполяции x 0 , x 1 , , x n и

гладкости функции f .

8

Замечание 1.7.

Из доказательства

теоремы 1.1 получаем, что

ξ = ξ(x) удовлетворяет условию

min{x, x0 , x1,…, xk−1}< ξ < max{x, x0. , x1,…, xk−1}.

Следствие 1.2.

Пусть

f C ( n +1) [a; b ]

и

M n+1 = max

f ( n+1) ( x)

.

Тогда

M n +1

a≤x≤b

rn ( x ; f )

≤

ω ( x )

для любого x [a ; b ], (1.11)

( n + 1)!

M n +1

max

rn ( x ; f )

≤

max

ω ( x )

,

(1.12)

( n + 1)!

a ≤ x ≤ b

a ≤ x ≤ b

max

rn ( x ;

f )

≤

M n +1

(b − a ) n +1 .

(1.13)

( n + 1)!

a ≤ x ≤ b

Важный пример (погрешность линейной полиномиальной интер-

поляции).

Пусть

n = 1

и

f C ( 2) [a;b].

Обозначим

h = x1 − x0,

q =

x − x0

, 0 ≤ q

≤ 1. В этом случае ω(x)

= (x −x )(x −x ) и интерпо-

h

0

1

ляционный многочлен может быть записан в виде

L1 ( x ;

f ) = (1 − q ) f ( x 0 ) + q f ( x1 ).

(1.14)

Наибольшее значение

ω(x)

на отрезке [x 0 ; x1 ] достигается в точке

τ = ( x0 + x1 ) 2 и

max

ω( x)

=

ω(τ)

=

1

( x1

− x0 )( x0

− x1 )

= h 2 .

x0 ≤x ≤ x1

4

4

Отсюда по формуле (1.12) получаем максимальную оценку погреш-

ности линейной интерполяции

h 2

max

r1 ( x ; f )

≤

max

f ′′( x )

.

(1.15)

8

x 0 ≤ x ≤ x 1

x 0

≤ x ≤ x 1

Если

f C ( 2 ) [a; b ], то оценка (1.15) не имеет места.

Найдем оценку погрешности интерполяции при минимальных требо-

ваниях к гладкости функции

f .

[a ; b ], то есть

Пусть

f удовлетворяет условию Липшица на отрезке

для любых ξ, η [a ; b ] выполнено неравенство

f (ξ) − f (η)

≤ α

ξ − η

,

где α = const ≥ 0 . Тогда максимальная оценка погрешности линейной интерполяции имеет вид

max		r1 ( x ; f )		≤ α	h	.	(1.16)
					2
x 0 ≤ x ≤ x 1

Действительно, используя формулу (1.14), имеем

9

f ( x ) − L1 ( x ; f ) = f ( x ) − (1 − q ) f ( x 0 ) − qf ( x1 ) = = (1 − q ) f ( x 0 ) + qf ( x1 ) − (1 − q ) f ( x ) − qf ( x ) ≤

≤ (1 − q ) f ( x 0 ) − f ( x ) + q f ( x1 ) − f ( x ) .

Поскольку f ( x1 ) =

f ( x 0

+ h ) и f ( x ) = f ( x 0 + qh ), имеем

f ( x1 ) − f ( x )

=

( f ( x 0 + h ) + f ( x 0 + qh )

≤ α (1 − q ) h ,

f ( x 0 ) − f ( x )

≤ α qh .

Так как (1 − q )q ≤ 1

4

при 0 ≤ q ≤ 1 , то

≤ 2 (1 − q ) q α h ≤ α h

f ( x ) − L1

( x ; f )

2

.

Оценка (1.16) доказана.

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

x 0 , а правый – за x1 , и	f ( x ) вычисляют по формуле линейной интерполя-
ции (1.14). Погрешность интерполяции для f C ( 2 ) [a; b ] оценивается по
формуле (1.15).
Замечание 1.8. Для класса n + 1 раз непрерывно дифференцируе-
мых функций константу	M n+1 = max		f ( n+1) ( x)		в (1.11)–(1.13) улучшить
	a≤x≤b

(уменьшить) нельзя, так как для случая f Ρ ( n +1 ) неравенство (1.11) пре-

вращается в равенство. Что касается величины max ω ( x ) , то она суще-

a ≤ x ≤b

ственно зависит от выбора узлов интерполяции и, следовательно, может быть уменьшена при специальном выборе узлов.

1.4. Интерполяционный многочлен Ньютона

Удобным представлением интерполяционного многочлена для практических вычислений является его запись в виде интерполяционного мно-

гочлена Ньютона.

Введем в рассмотрение разделенные разности функции f для заданной табл. 1.1.

Значения функции f в узлах xi , i = 0,…, n будем называть раз-

деленными разностями нулевого порядка.

Для любой пары узлов x i , x j величины

10