УП_Вычисл_матем_Кузина-Кошев

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

Вкачестве базисных рассмотрим полиномы степени n, удовлетворяющие следующему условию:

		1,	если x x j ,		j k;
		lkn (x)		если x x j ,	j k.
		0,
Такому условию удовлетворяют полиномы вида:
lkn (x)		(x x0 ) (x x1 ) ... (x xk 1 ) (x xk 1 ) ... (x xn )
							(39)
	(xk x0 ) (xk x1 ) ... (xk			xk 1 ) (xk xk 1 ) ... (xk xn )
d0 d1		x d2 x2 ... dn xn .

Тогдав качестве интерполяционного полинома можно взять сумму вида

n	n
Ln (x) yk lkn (x) . Имеем	Ln (xi ) yk lkn (xi ) yi . Таким образом,
k 0	k 0

значения функции Ln (x) в узлах интерполяции xi совпадают с заданными значениями yi .

Полином

	n	n	(x xi )
	Ln (x) yk			(40)
			(xk xi )
	k 0	i 0
		i k

называется интерполяционным полиномом Лагранжа.

Оценка погрешности интерполяции полиномом Лагранжа

Если известно, что интерполируемая функция (x) имеет непрерывные частные производные до (n 1) -го порядка включительно, то разность между точным значением функции и её полиномиальным приближением

(x) Ln (x) в любой точке х из интервала интерполирования

x [a,b]

можно оценить по формуле (x) Ln (x)

(n 1)

( )

n

(x x j ) , где

[a,b].

(n 1)!

j 0

Из этой формулы вытекает, что

(n 1)

(x)

n

Rn (x)

(x) Ln (x)

max

x x j

.

(41)

x [a,b]

(n 1)!

j 0

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

Обобщением понятия производной является понятие разделенной разности. Разделенная разность нулевого порядка совпадает со значением

91

самой функции y(xi ) . Разделенной разностью первого порядка называется

величина y(x , x

)

y(xi ) y(x j )

.

j

i

xi x j

Аналогично,

величина

y(xi , x j , xk )

y(xi , x j ) y(x j , xk )

называется

xi xk

разделенной

разностью

второго

порядка,

а

y(x0 , x1 ,..., xk )

y(x0 , x1 ,..., xk 1 ) y(x1 , x2 ,..., xk )

–

разделенной

разностью

x0 xk

порядка k.

y(xi ) , i 1, , n

Таблицу разделенных

разностей

функции

можно

представить в виде следующей схемы:

y(x1 )

y(x1 , x2 )

y(x2 )

y(x1 , x2 , x3 )

y(x2 , x3 )

y(x1 , x2 , , xn )

y(x3 )

y(xn 1 , xn )

y(xn )

Предположим, что функция y представляет собой полином y Pn (x) порядка n, тогда разделенные разности также являются полиномами:

P(x, x0 ) P(x) P(x0 ) – полином порядка ( n 1), x x0

P(x, x0	, x1 )		P(x, x0 ) P(x0 , x1 )				– полином порядка ( n 2 ),
				x x1
P(x, x0	,..., xk	)		P(x, x0 ,..., xk 1 ) P(x0 ,..., xk )			– полиномпорядка(n k 1 ).
				x xk
						можно записать полином Pn (x) в
Используя приведенные разности,
следующей форме:

P(x) P(x0 ) P(x, x0 )(x x0 )

P(x0 ) (x x0 )P(x0 , x1) (x x0 )(x x1)P(x, x0 , x1)

P(x0 ) (x x0 ) P(x0 , x1 ) (x x0 )(x x1 ) P(x0 , x1 , x2 ) ...

(x x0 )(x x1 )...(x xn 1 ) P(x0 , x1 ,..., xn )

n

P(x0 ) (x x0 ) (x x1 ) ... (x xk 1 ) P(x0 , x1 ,..., xk ).

k 1

92

Если такой полином Pn (x) взять в качестве интерполяционного полинома для функции y, то, так как в узлах интерполяции Pn (xi ) yi , совпа-

дут и все разделенные разности, вычисленные в этих точках для функций Pn (x) и y. Тогда интерполяционный полином можно записать в виде:

	n
Pn (x) y0	(x x0 ) (x x1 ) ... (x xk 1 ) y(x0 , x1 ,..., xk ) .	(42)

k1

Всилу способа построения Pn (xi ) yi , а Pn (x) представляет собой

интерполяционный полином, называемый интерполяционным полиномом Ньютона.

Если узлы интерполяции на интервале [a, b] расставлены равномерно, т.е. шаг интерполяции h xi xi 1 для всех i 0,1, ..., n , то интерполяционный полином Ньютона можно записать в следующем виде:

	P (x) y				0		y0
	n						h 1!
		n y0		(x x0 )
		hn n!

(x x	0	)	2 y0		(x x	0	) (x x ) ...
			h2 2!				1
							(43)

... (x xn 1 ).

Здесь k yi k 1 yi 1 k 1 yi – неразделенные разности, которые вычис-

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

Таким образом, для построения интерполяционного полинома Ньютона необходимо предварительно найти (n 1) неразделенных разностей.

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

P(x) a0 a1 x a2 x2 ... an xn

к виду

P(x) a0 x (a1 x (a2 ... x (an 1 x an ))) ,

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

Пример 2 9 .

По данным табл. 1 вычислить приближенное значение функции f(x) в точке x 21,85 средствами MS Excel с использованием интерполяционного

полинома Ньютона.

				Таблица 1
i	0	1	2		3
xi	–0,7	5,5	11,7		17,9
yi	25,25	20,03	44,93		10,61

93

Решение.

Интерполяционный полином Ньютона (43) будет

P (x) y	0		y0 (x x			)		2 y0		(x x			)(x
4			h 1!		0			h2 2!				0
			3 y0	(x x			)(x x )(x x				2	).
			h3 3!			0		1

иметь вид:

x1 )

(44)

Исходные данные из табл. 1 введем в ячейки A1:E3 (рис. 24). В ячейку G2 введемшаг, вI2 – точкуx, вкоторойтребуетсявычислитьприближенное значение функции f(x).

Рис. 24. Вычисление значения таблично заданной функции в MS Excel

94

Вячейку B4 занесем формулу для вычисления конечной разности первого порядка и скопируем ее во все ячейки, где будут вычисляться конечные разности второго и третьего порядков.

В7-ю строку введем формулы для вычисления промежуточных коэффициентов в соответствии с формулой (44), затем просуммируем их в ячейке E8 и вычислим в ячейке E9 значение функции в заданной точке. На рис. 24, а показан лист Excel с формулами (для этого задать: Формулы – Зависимости формул – Показать формулы) и числовыми значениями – рис. 24, б.

Построение полинома Ньютона и графика с использованием MathCAD для этого примера показано в лабораторном практикуме (см. пример 15 [21]):

P (x) 25,25 5,22			(x 0,7)	30,12	(x 0,7)(x 5,5)
4	6,2	1!		6,22 2!

89,34 (x 0,7)(x 5,5)(x 11,7) 6,23 3!

0,062x3 1,423x2 5,99x 20,338.

ВычисленноевMathCAD значение функции f (21,85) 83,151 несколь-

ко отличается от полученного в Excel, что говорит о погрешности вычислений.

6.3. Принципы двухмерной интерполяции

Предположим, что функция f (x, y) зависит от двух переменных и известны ее значения при некоторых значениях аргументов {xi , y j }. Предположим, что нам необходимо вычислить значение функции в

межузловой точке		f (xr , ys ) ,	причем r i ,	s j для любых i и	j . В этом
случае для	всех	значений	y j строятся	интерполяционные	полиномы
P(x, y j ) , т.е.	P(x, y1 ), P(x, y2 ), ..., P(x, ym ) .			Далее для последовательности

точек P(xr , y1 ), P(xr , y2 ), ..., P(xr , ym )			строится интерполяционный полином
PP xr , y , и затем вычисляется требуемое значение					f (xr , ys ) PP xr , ys .
	Пример 3 0
	Вычислить значение функции двух переменных z					= f(x, y) в точке
z = (0,25; 2,5), если функция задана табл. 2.							Таблица 2
	yi/xi	0,1		0,2			0,3
	1	2		3			4
	2	3		4			5
	3	4		5			6

95

Решение.

Построим интерполяционные полиномы для всех значений yi по формуле (44), используя значения неразделенных разностей из табл. 3.

hx 0,1 ; принимаем x0					0,2.						Таблица 3
xi	0,1	0,2	0,3		xi	0,1	0,2	0,3		xi		0,1	0,2	0,3
z(xi, y1)	2	3	4		z(xi, y2)	3	4	5		z(xi, y3)		4	5	6
y1	1	1			y2	1	1			y3		1	1
2y1	0				2y2	0				2y3		0

P

(x, y ) 2

1

(x 0,1) 10x 1;

P (0,25; 3) 3,5;

1

0,1 1!

P

(x, y2 ) 3

1

(x 0,1) 10x 2;

P (0,25; 4) 4,5;

0,1 1!

P

(x, y3 ) 4

1

(x 0,1) 10x 3;

P (0,25; 5) 5,5.

0,1 1!

Для

последовательности точек P (0,25, y1 );

P (0,25, y2 ); P (0,25, y3 )

строим

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

полином P (0,25, y1 ) , используя

значения

неразделенных разностей из табл. 4.

hy 1; принимаем y0

4.

Таблица 4

yi

1

2

3

P(x)

3,5

4,5

5,5

P(x)

1

1

2 P(x)

0

P (0,25, y) 3,5

1

(y 1) 2,5 y

1 1!

Вычисляем требуемое значение P (0,25; 2,5) 5.

6.4.Кубическая сплайн-интерполяция

Внекоторых задачах требуется интерполировать не только саму функцию, но также ее первую и вторую производные. Для решения таких задач используется сплайновая интерполяция. Сплайн – это функция, которая на каждом межузловом интервале совпадает с некоторым полиномом, своим для каждого интервала. Полиномы смежных интервалов склеиваются так, чтобы функция и ее производная были непрерывны.

Вслучае кубической сплайн-интерполяции на каждом из интервалов

[xi 1 , xi ] , где i 1, 2, ..., n , строится интерполирующий кубический полином

96

(полином третьего порядка) таким образом, чтобы в точках xi выполнялись следующие условия:

1)P(xi ) yi ;

2)P (xi 0) P (xi 0) ;

3)P (xi 0) P (xi 0) ,

где P (xi 0) и P (xi 0) – соответственно первые и вторые производные функции P(xi ) слева и справа.

Кубический полином, определенный на интервале [xi 1 , xi ] , имеет вид

P(x) a

i

b (x x

i 1

) c

i

(x x

i 1

)2

d

i

(x x

i 1

)3

,

(45)

i

где x [xi 1 , xi ] . Для всего интервала будет соответственно n кубических полиномов, отличающихся коэффициентами ai , bi , ci , di .

Требование выполнения условий 1, 2, 3 приводит к следующим линейным уравнениям, составленным для каждого из n интервалов:

P(xi 1 ) ai ;		i 1, 2, ..., n;									(46)
P(x ) a b h c h 2 d					h 3	y	;	h x x		;
				i
i	i	i i i i			i	i		i	i i 1
	2 ci (x xi 1 ) 3		di (x xi 1 )					2
bi									P (x);
2 ci 6 di (x xi 1 ) P (x).

На интервале [xi , xi 1 ]

P(x) ai 1 bi 1 (x xi ) ci 1 (x xi )2 di 1 (x xi )3 .

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

bi 1 2			ci 1 (x xi ) 3 di 1	(x xi )2 P (x);	(47)
	ci 1
2		6 di 1 (x xi ) P (x).

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

Приравнивая производные для P(х) первого и второго порядка в конце i-го и в начале (i 1) -го интервалов, получим набор уравнений:

b

2 c

h

3 d

h 2 b

2 c

h

3

d

h 2

;

i 1, 2, ..., n 1.

i

i

i

i

i

i 1

i 1

i 1

i 1

i 1

2 ci 6 di

hi 2 ci 1 6 di 1 hi 1;

Эти соотношения дают

(2n 2)

уравнения,

всего

получим (4n 2)

уравнения.

97

4) В точках x0 и xn считается, что f (x0 ) f (xn ) 0 . Это дает еще два уравнения:

с	0;	(48)
1
сn 3 dn hn 0.

Таким образом, имеем 4n линейных уравнений с 4n неизвестными: ai, bi, ci, di, i = 1, …, n. Решив эту систему, найдем совокупность из n кубических полиномов, значения которых приближенно равны значениям интерполируемой функции в соответствующих точках.

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

Замечание. Чаще всего узлы сплайновой интерполяции располагают равномерно, т.е. hi const , хотя это необязательно.

	Пример 3 1 .
	Построить кубический		сплайн для		функции,	заданной таблично
(табл. 5):							Таблица 5
	i		0		1		2
	xi		1		2		3
	yi		1		5		8

Определить значения функции в точках x 1,25 и x 2,5. Решение.

Весь интервал делим на подынтервалы: [1, 3] = [1, 2] [2, 3], на каждом из которых строим полином 3-го порядка.

На интервале [1,2], в соответствии с формулами (45, 47), получаем:

y a1 b1 (x 1) c1 (x 1)2 d1 (x 1)3 ,

y b1 2c1 (x 1) 3d1 (x 1)2 , y 2c1 6d1 (x 1).

На интервале [2,3] имеем:

y a2 b2 (x 2) c2 (x 2)2 d2 (x 2)3 ,

y b2 2c2 (x 2) 3d2 (x 2)2 , y 2c2 6d2 (x 2).

Составим систему уравнений. Для каждого интервала по формулам (46) получим уравнения, соответствующие прохождению полинома через

98

начальную и конечную точки. Здесь узлы интерполирования располагаются равномерно, h 1. Для интервала [1,2] получаем:

1 a1;

5 a1 b1h c1h2 d1h3 a1 b1 c1 d1,

для интервала [2,3]:

5 a2 ;

8 a2 b2h c2h2 d2h3 a2 b2 c2 d2 .

Приравниваем в точке x 2 , являющейся концом 1-го и началом 2-го интервалов, соответствущие производные первого и второго порядков:

b1 2c1 3d1 b2 ;

2c1 6d1h 2с2 .

По формулам (48) получаем еще 2 уравнения: c1 0; c2 3d2 0.

Полученную систему из 8-ми уравнений с 8-ю неизвестными можно решить любым известным методом. Приведем решение в системе MathCAD:

a1 0

b1 0

c1 0

d1 0

a2 0

b2 0

c2 0

d2 0

Given

a1

1

a1 b1 c1

d1

5

a2

5

a2 b2 c2

d2

8

b1 2c1

3d1

b2

1

2c1 6d1

2c2

4.25

c1

0

0

2c2 6d2

0

0.25

Find(a1 b1 c1

d1 a2 b2 c2 d2)

5

3.5

0.75

0.25

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

y 1 4,25(x 1) 0,25(x 1)3 0,25x3 0,75x2 3,5x 3

y 5 3,5(x 2) 0,75(x 2)2 0,25(x 2)3 0,25x3 2,25x2 9,5x 7 .

Определим значение функции в точке x 1,25:

y(1,25) 0,25 1,253 0,75 1,252 3,5 1,25 3 2,059

и в точке x 2,5:

y(2,5) 0,25 2,53 2,25 2,52 9,5 2,5 7 6,594 .

Графики кубических сплайнов в системе MathCAD показаны на рис. 25.

99

y1(x)

y2(x)

y1(x) 1 4.25(x 1) 0.25(x 1)3

y2(x) 5 3.5(x 2) 0.75(x 2)2 0.25(x 2)3

8

6

4

2

0

2

0.5	1	1.5	2	2.5	3	3.5
			x

Рис. 25. Графики кубических сплайнов в системе MathCAD

Контрольные вопросы и задания

1.Какие точки называются узлами интерполирования?

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

3.Что такое интерполянта?

4.В чем заключается задача обратного интерполирования?

5.В каком методе полиномиальной интерполяции используются конечные разности?

6.Какие виды интерполирования, кроме полиномиального, Вам известны?

7.Что такое кубическая сплайн-интерполяция?

8.Можно ли по экспериментальным данным получить аналитическую зависимость?

9.В каких случаях интерполяционные полиномы Ньютона и Лагранжа совпадают? Приведите пример.

100