Вычислительная математика лекции

3.Интерполирование функций.

3.1.Постановка задачи. Алгебраическое интерполирование.

Интерполирование функций – один из вариантов аппроксимации, иными словами замены исходной функции другой,

близкой функцией, удобной для проведения расчетов.

Аппроксимируемая функция задается в табличной форме (сеточная функция) {xi , y(xi ) yi}in 0 ; x [a,b]. В качестве аппроксимирующей функции часто берут обобщенный полином

n (x) an n (x) an 1 n 1(x) ... a0 0 (x),

где {i (x)}in 0 заданный набор базисных функций, ai i=0,…,n

подлежащие расчету коэффициенты.

Интерполирование использует критерий близости в виде

				n (xi ) yi ;			i=0,n.
Для определения n+1 неизвестных ai								i=0, n нужно решить
систему n+1 линейных алгебраических уравнений с n+1
неизвестными Pa=y,
n (x0 )				0 (x0 )				an		y0
			(x )			(x )		a		y
P		n	1		0	1		, a= n 1		, y= 1	.
	n
			(xn )	0 (xn )				a0		yn
Введем векторы φj=(φj(x0) , φj(x1) ,…,φj(xn))T,									j=0,n. Будем

говорить, что система функций линейно независима на узлах {xi }in 0 ,

n
если i j	0 тогда и только тогда, когда γi=0, i=0,n, j=0, n. В
i 0

этом случае решение системы существует и единственно.

21

Интерполирование называется алгебраическим, если

i (x) xi , i=0,n. Такая система базисных функций линейно независима на несовпадающих узлах. Отсюда в частности следует,

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

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

Пусть заданы узлы x0,…,xn, среди которых нет совпадающих, и

значения y0,…,yn функции y(x) в этих узлах . Существует один и только один алгебраический многочлен P m(x) степени m n,

принимающий в заданных узлах xk значения yk=y(xk).

Из теоремы в частности следует, что если y(x)=Pm(x) и m n, то интерполяционный полином Pn(x), построенный на n+1 узлах совпадает с Pm(x).

В дальнейшем ограничимся анализом алгебраической интерполяции.

3.2. Формулы алгебраической интерполяции.

Ниже приводятся формулы, позволяющие записать алгебраический интерполяционный полином без составления и решения системы алгебраических уравнений. Полагаем, что задана сеточная функция с n+1 узлами {xk , y(xk ) yk }kn 0 .

3.2.1 Интерполяционный полином в форме Лагранжа.

n

Ln(x)= wk (x) yk .

k 0

22

Полином wk(x) будем искать в следующей форме

			n					1, если x=x		k
wk (x) ck			(x xi ). Потребуем wk
							(x)				, i k, i=0, n.
			i 0,i k						0, если x=x
										i
		Константа ck определяется из условия wk(xk)=1.Откуда
ck			1	.
		n
		(xk xi )
		i 0,i k
					n		xi
		Таким образом, wk			(x)	x		. Нетрудно убедится, что
						x	x
					i 0,i k k		i

wk(xi)=0, когда i≠k, i 0, n.

Окончательный вид интерполяционного полинома в форме Лагранжа

	n	n	x	xi
	Ln (x) yk				.
			x
	k 0	i 0,i k		x
			k	i

Ln(x) – интерполяционный полином, построенный на n+1 узлах.

Действительно, Ln(x) есть полином степени не выше n, так как wk(x)

– полином n – ой степени, а, кроме того, Ln(xk)=yk ( k 0, n ).

Эквивалентный вариант записи формулы Лагранжа:

	n	(x)		n
	L(x) yk			, (x) (x xi ).
		(x x )	'(x )
	k 0			i 0
		k	k

Упражнения.

1.Записать интерполяционный полином a)L1(x); b)L2(x).

2.Дана таблица функции y=ln(x)

i	0	1	2	3	4
x	1.0	1.1	1.2	1.3	1.4
y	0.000000	0.095310	0.182322	0.262364	0.336472

23

Вычислить приближенное значение ln(1.23) для а) линейной интерполяции;

б)квадратичной интерполяции.

Оценить погрешности.

3.Построить L3(x) для табличной функции

x	0	2	3	4
y	2	0	1	6

3.2.2. Интерполяционный полином в форме Ньютона.

Nn (x) y(x0 ) (x x0 ) y(x0 , x1) (x x0 ) x x1 y(x0 , x1, x2 ) ... (x x0 )...(x xn 1 ) y(x0 , x1,..., xn )

n

y(x0 , x1,..., xk ) k (x);

k 0

k 1

0 (x) 1; k (x) (x xi ).

i 0

Докажем, что Nn(x) есть интерполяционный полином. Действительно это полином не выше n – ой степени. Кроме того,

для любого xi (i 0, n )

Nn (xi ) f (x0 ) (xi x0 ) y(x0 , x1) (xi x0 ) xi x1 y(x0 , x1, x2 ) ...

(xi x0 )...(xi xi 1) y(x0 , x1,..., xi ) y(xi )

Последнее равенство следует из формулы, связывающей значения сеточной функции с разделенными разностями.

В отличие от формы Лагранжа при добавлении узла в сеточной функции форма Ньютона не требует пересчета всех её слагаемых:

Nn+1(x)=Nn(x)+y(x0,x1,…,xn,xn+1)ωn+1(x).

В силу теоремы о существовании и единственности алгебраического интерполяционного полинома при любом способе вычисления получают идентичные полиномы.

24

3.3.Остаточный член интерполяционного полинома. Оценка погрешности интерполирования.

Назовем остаточным членом интерполяционного полинома функцию

rn(x)=y(x)-Pn(x).

Если для y(x) известны только значения в узлах, то никаких суждений о r(x) сделать невозможно. Дополнительную информацию можно получить, если предположить y(x) C[na,b1].

Рассмотрим расширенную систему узлов x0,x1,…xn,xn+1=x (x [a,b]). Для вычисления y(x) применим формулу, связывающей значения сеточной функции с разделенными разностями

y(x)=y0+(x-x0)y(x0,x1)+…+(x-x0)…(x-xn-1)y(x0,…xn)+(x-x0)…(x- xn)y(x0,…,x)=

=Nn(x)+(x-x0)…(x-xn)y(x0,…,x).

Отсюда rn(x)=y(x)-Nn(x)=(x-x0)…(x-xn)y(x0,…,xn,x).

Учитывая, что y(x , x ,..., x , x)							y(n 1)	( )	,	[a,b], получаем
			0		1	n	(n 1)!
r (x)		y(n 1)	( )	,	[a,b].
	n 1	(n
n			1)!

Задача оценки погрешности состоит в отыскании max | rn (x) | .
						x [a,b]
max | rn	(x) | max | n 1	(x) |	Mn 1		, где Mn 1	max | y	(n 1)	(x) | .
			(n 1)!
x [a,b]	x [a,b]					x [a,b]

3.4. Минимизация погрешности интерполирования.

Пусть число узлов (n+1) фиксировано и о функции y(x)

известно лишь, что она имеет (n+1) ограниченных производных.

Тогда единственным средством уменьшить погрешность является

25

выбор узлов сетки. Этот выбор должен удовлетворять решению

		n
характерной задачи на минимакс min	max | (x xi ) | . Из всех
x0 ,x1 ,...,xn a x b		i 0

возможных стандартных полиномов ωn+1(x) путем выбора корней нужно подобрать такой, чтобы модуль его максимального значения на заданном интервале был минимальным. Эта задача, известная как задача о построении полинома, наименее уклоняющегося от нуля,

была решена (по другому поводу) П.Л. Чебышевым, который доказал, что такими свойствами обладают нормированные полиномы, известные под его именем.

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

(n+1) на интервале [a,b].

Свойства полиномов Чебышева.

1. Полином Чебышева n – ой степени Tn(t)=cos(n arccos(t)) (t [-

1,1]) можно рассчитать по рекурсивным формулам T0(t)=1; T1(t)=t; Tn(t)=2tn-1(t)-Tn-2(t); n=2,3,… .

2. На отрезке [-1,1] полином Tn(t) имеет n вещественных корней

tk cos (2k 1) ; k=0, n 1. 2n

3. Коэффициент при tn равен 2n-1.

4. При n≥0 max | Tn (t) | 1.

1 t 1

Если n≥1, то этот максимум достигается ровно в n+1 точках,

которые находятся по формуле tm		m	, m=0, n. При этом
	cos
		n

Tn(tm)=(-1)m, т. е. максимумы и минимумы многочленов Чебышева чередуются.

Предположим, интерполяционный полином строится на отрезке [-1,1] . Тогда для минимизации погрешности нужно

26

выбрать узлы, совпадающие с корнями полинома Чебышева. То есть

	(t)	1	T	(t). Тогда max | r (t) |	M n 1	. Произвольный
n 1		2n	n 1	1 t 1 n	(n 1)!2n

отрезок [a,b] приводится к стандартному [-1,1] заменой переменных

a b

b a

2

a b

x

t; t=

x

.

При этом оказывается

2

2

b

a

2

n

n

b a

a b

b a n 1

n

n 1 (x) (x xi )

t

xi

t ti

2

2

i 0

i 0

2

i 0

2

a b

ti

=

xi

;

b a

2

x

a b

b a

t .

i

2

2

i

Для минимизации погрешности значения ti должны совпадать

корнями Tn+1 полинома Чебышева : t

cos

(2i 1)

;

k=0,n.

i

2n

Таким образом, выбор чебышевских узлов приводит к следующему результату

	(t) |	M	n 1		b a n 1
max | rn							.
		(n 1)!2		n
a x b						2

Можно показать, что выигрыш по сравнению с равномерным

4 n

выбором узлов составляет раз.

e

3.5. Сходимость интерполяционного метода.

Будем говорить, что при заданной стратегии выбора узлов

			0 .
метод интерполяции сходится, если Pn= max	y(x) Pn	(x)
x	a,b	n
27

Из конструктивной теории функций известно, что для заданной степени n существует полином Qn(x) наилучшего приближения такой, что

min max		y(x) Qn	(x)		En ( y)	. Здесь En(y) наименьшая
Qn x a,b
достижимая погрешность. Алгоритм построения такого полинома
настолько сложен, что его практическое использование
нецелесообразно.

Известна асимптотика En(y): En(y)=o(1/nr), где r наибольший порядок ограниченной производной y(x). В этом случае En(y)→0 при n→∞.

Интересно выяснить, насколько близка погрешность интерполяции к En(y).

Теорема. Пусть Pn(x) интерполяционный полином для y(x) на интервале [a,b] с чебышевской сеткой. Тогда Pn (1+ln n) En(y).

Таким образом, погрешность интерполяционного полинома с чебышевской сеткой близка En(y).

Для интерполяционного полинома с равномерной сеткой аналогичная оценка будет иметь величину порядка 2n.

При n→∞ в случае использования чебышевской сетки Pn→0

даже, если r=1; для равномерной сетки в аналогичных условиях

Pn→∞.

Рассмотренная выше ситуация обобщается следующими теоремами:

Теорема1. Какова бы не была стратегия выбора узлов,

найдется непрерывная на интервале [a,b] функция такая, что метод интерполяции расходится, то есть

Pn(y)→∞ при n→∞.

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

28

Теорема2. Пусть в качестве узлов интерполяции выбирают чебышевские узлы. Тогда для любой функции, обладающей ограниченными производными на [a,b], метод интерполяции сходится.

Ниже приведен пример функции с несуществующей первой производной внутри заданного интервала.

3.6. Устойчивость интерполяционного полинома относительно погрешности вычисления y(x).

Пусть yi вычисляются с погрешностью y*i=yi+δi, причем | δi| δ.

Используем лагранжеву форму интерполяционного полинома

n				n
Ln(x)= yi wi	. Тогда	L*n	max | L*n (x) Ln (x) | (	max | wi (x) |) .
i 0			x [a,b]	x [a,b]i 0

n
Величина Λn=max | wi (x) | называется константой Лебега Она не
x [a,b] i 0
зависит от [a,b], но только от расположения узлов Для
чебышевской сетки	Λn=(2/π) ln(n+1)+1 Для
	29

	2n 1
равномерной сетки Λn>	(2n 1) n (n≥4 Из формул следует, что

интерполяционный полином с чебышевскими узлами малочувствителен к вычислительным ошибкам.

3.7. Интерполяция сплайнами.

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

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

реализующий одну из модификаций указанной идеи.

Пусть отрезок [a,b] разбит точками a=x0 < x1 < … xn =b на частичные отрезки [xi-1 ,xi ] i=1,2,…,n. Сплайном степени m

называется функция Sm (x), обладающая следующими свойствами:

1. Sm (x) непрерывна на [a,b] со всеми своими производными до некоторого порядка p.

2.На каждом частичном отрезке [xi-1 ,xi ] функция Sm (x)

совпадает с некоторым алгебраическим полиномом Pmi (x) степени m.

Разность m-p называется дефектом сплайна. Наибольшее распространение получили сплайны 3-ей степени (кубические) с

дефектом 1 или 2.

Сплайн называется интерполяционным, если Sm (xi) =yi для i=0,1,…,n.

30