45

q(q −1), f (x0 ; x1; x2 ; x3 )(x − x0 )(x − x1 )(x − x2 )=

Параметрами программы являются вектора исходных данных x и y . Компоненты вектора z равны: zi = xi − xT , n - число точек табличного задания функции y = f (x). Для получения матрицы по схеме Эйткена необходимо определить лишь вектор z и саму матри-

цу M :

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

Если

интерполируемая

функция

задана

на

таблице

с постоянным шагом

h (xi = x0 + ih, i = 0,1,..., n), то можно использовать связь между конечными и разделенными

разностями:

f (x

; x

;...; x

)=

k yi

. В этом случае многочлен Ньютона можно записать не-

i

i+1

i+k

hk k!

сколько в ином виде:

Пусть

q =

x − x0

,

h = x − x

0

= const, P (x)= y

0

+ f (x

0

; x

)(x − x

0

)+

1

n

1

h

+ f (x0 ; x1; x2 )(x − x0 )(x − x1 )+ f (x0 ; x1; x2 ; x3 )(x − x0 )(x − x1 )(x − x2 )+...

Преобразуем разделенные

разности

в

конечные:

f (x

0

; x

)(x − x

0

)= y0 q,

f (x

0

; x

; x

2

)(x − x

0

)(x − x )=

1

1!

1

1

=

2 y

q

(x − x

0

)+

(x

0

− x )

=

2 y

0

1

0

2!

h

2!

f (x

; x

)=

y0

,

тогда

0

1

1!h

2 y

0

(x − x

0

) (x − x )

=

1

2!

h

h

=

3 y0 q(q −1)

x − x0 + x0 − x2

=

3 y0 q(q −1)(q − 2) и так далее.

3!

h

3!

Тогда многочлен Ньютона можно переписать в следующем виде:

P

(x)= y

0

+

y0

q +

2 y0

q(q −1)+

3 y0

q(q −1)(q −2)+... +

n y0

q(q −1)(q −2)...(q −n +1). (2.10.1)

n

1!

2!

3!

n!

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

P (x)= y

n

+

yn−1

q +

2 yn−2

q(q +1)+

3 yn−3 q(q +1)(q + 2)+...

n

1!

2!

3!

(2.10.2)

n y0

q(q +1)(q + 2)...(q + n −1).

+

n!

x

sin t

dt в точках x(1) = 0.1 и x(2) = 0.7 , используя форму-

Пример. Вычислить

six = ∫

t

0

лы (2.10.1) и (2.10.2)

при

n = 3, если

x

six

y

2 y

3 y

4 y

5 y

0.0 0.00000

		19956
0.2	0.19956	-266
		19690

0.40.39646 -523

19167

0.60.58813 -770

		18397
0.8	0.77210	-999
		17398

-257
	10
-247	8

18 -229

1.0

0.94608

Здесь x0 = 0.0,

h = 0.2, q =

x − x0

. Для

x

(1)

используем формулу (2.10.1), так как значение

h

0.1 − 0.0

1

x(1) = 0.1 расположено в начале таблицы, тогда q(1) =

=

.

0.2

2

si(0.1)≈ si(0.0)+

y0

q(1) +

2 y0

q(1)(q(1) −1)+

3 y0

q(1)(q(1) −1)(q(1) − 2)=

1!

2!

3!

1

1

1

1

(−

0.00266)+

1

1

1

3

(− 0.00257)=

= 0.0 +

0.19956 +

−

−

−

2

2

2

2

6

2

2

2

= 0.0 + 0.09978 + 0.0003325 − 0.0003212 = 0.10270.

Для x(2) = 0.7 лучше использовать формулу (2.10.2), так как

x(2) = 0.7 расположено

ближе к нижнему краю таблицы. Тогда xn =1.0, h = 0.2,

q

(2)

=

x(2) − xn−1

=

0.7 −0.8

= −0.5.

h

0.2

46

si(0.7)≈ si(0.8)+	y	3	q(2) +	2 y	2	q(2)(q(2) +1)+	3 y	q(2)(q(2) +1)(q(2) + 2)=
							1
	1!			2!			3!

=0.77210 + (− 0.5) 0.18397 + (− 0.5) (0.5) 12 (− 0.00770)+ (− 0.5) (0.5) 1.5 16 ×

×(− 0.00247)= 0.77210 − 0.091985 + 0.0009625 + 0.0001544 = 0.68123.

2.11.Лабораторная работа № 3. Интерполирование и экстраполирование данных. Интерполяционный многочлен Ньютона

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

Интерполяционные многочлены Ньютона нужных степеней строятся по формулам

(2.7.1) или (2.10.1) и (2.10.2).												y = f (x)
	Пример. Найти приближенное значение функции													при данном значении ар-
гумента xT		с помощью интерполяционного многочлена Ньютона, если функция задана в не-
равноотстоящих узлах таблицы и xT									= 0.552 :
	x			0.35		0.41			0.47		0.51		0.56			0.64
	y			2.73951		2.30080			1.96464		1.78776		1.59502			1.34310

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

f (x ; x

2

)=

f (x2 )− f (x1 )

,

f (x

; x

2

; x

3

)=

f (x2 ; x3 )− f (x1; x2 )

,

1

x2

− x1

1

x3

− x1

... f (x1; x2 ; x3 ; x4 ; x5 ; x6 )= f (x2 ; x3 ; x4 ; x5 ; x6 )− f (x1; x2 ; x3 ; x4 ; x5 ). x6 − x1

Запишем получаемые конечные разности в виде таблицы, однако, в целях экономии места, не будем писать разности между строк xi и xi+1 . Тогда получим

xi	f (xi )	f (xi ; xi+1 )	f (xi ;...xi+2 )	f (xi ;...xi+3 )	f (xi ;...xi+4 )	f (xi ;...xi+5 )
0.35	2.73951	-7.311833	14.243057	-15.227431	-102.232694	825.705728
0.41	2.30080	-5.602667	11.806667	-36.696296	137.221967
0.47	1.96464	-4.422000	6.302222	-5.135244
0.51	1.78776	-3.854800	5.429231
0.56	1.59502	-3.149000
0.64	1.34310

Далее по формуле (2.7.1) вычисляем

P5 (x)= ∑6 f (x1 ; x2 ;...; xk )ωk (x)= f (x1 )+ f (x1 ; x2 )(x − x1 )+... + f (x1 ; x2 ;...; x6 )(x − x1 )...(x − x5 )=

k =1

47

=2.739510 − 7.311833(x − 0.35)+14.243057(x − 0.35)(x − 0.41)−15.227431(x − 0.35)×

×(x − 0.41)(x − 0.47)−102.232694(x − 0.35)(x − 0.41)(x − 0.47)(x − 0.51)+ 825.705728 ×

×(x − 0.35)(x − 0.41) (x − 0.47) (x − 0.51)(x − 0.56)= 825.705728x5 − 2001.355769x4 +

+1898.621006x3 −870.466826x2 +184.903812x −11.051930.

Подставляя теперь значение xT = 0.552						вместо аргумента x в найденную формулу,
вычислим приближенное значение функции P5 (0.552)=1.62450.
Задание № 1. Построить интерполяционный многочлен Ньютона по неравноотстоя-
щей сетке узлов и найти приближенное значение интерполируемой функции y = f (x)													при
значении аргумента xT .
I.
	№ варианта		xТ			x				y
	1		0.154			0.05				4.4817
	6		0.257			0.19			4.9530
	11		0.291			0.21				5.4739
	16		0.335			0.27			6.0496
	21		0.367			0.32				6.6859
	26		0.412			0.34			7.3891
						0.39				8.1662
						0.45			9.0250
II.
	№ варианта		xТ			x				y
	2		0.015			0.01			9.9182
	7		0.195			0.11				9.5194
	12		0.137			0.16			9.1365
	17		0.245			0.23				8.8769
	22		0.385			0.28			8.4164
	27		0.478			0.39				8.0779
						0.46			7.7530
						0.50				7.4412
III.
	№ варианта		xТ			x				y
	3		1.330			1.30			4.7556
	8		1.455			1.45				5.3533
	13		1.813			1.65			6.4552
	18		2.290			1.90				7.5618
	23		2.750			2.40			8.6734
	28		3.040			2.55				9.7904
						2.80			10.9131
						3.20				12.0419

IV.

48

№ варианта		xТ			x				y
4		3.905			3.50				33.1154
9		4.110		4.55				34.8133
14		5.995		5.60				36.5982
19		7.213		6.20				38.4747
24		9.115			7.75				40.4473
29		10.050			8.80				42.5211
					9.45				44.7012
				10.95				46.9931

V.

№ варианта		xТ			x				y
5		1.090			1.01				12.6183
10		1.250			1.08				12.7644
15		1.275		1.11				12.9122
20		1.316		1.21				13.0617
25		1.488		1.26				13.2130
30		1.343			1.33				13.3660
					1.46				13.5207
				1.51				13.8357

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

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

мента xT1 = 0.4675 и xT2 = 0.5410.

x

0.45

0.46

0.47

0.48

0.49

0.50

y	20.1946	19.6133	18.9425	18.1746	17.3010	16.3123
x	0.51	0.52	0.53	0.54	0.55	0.56
y	15.1984	13.9484	12.5508	10.9937	9.2647	7.3510

Составим таблицу конечных разностей. Ограничимся при этом разностями четвертого порядка, так как они практически постоянны. Для вычисления значения функции при xT1 = 0.4675 воспользуемся формулой Ньютона для интерполяции вперед (2.10.1). При этом

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

x0

= 0.46, q =

x − x0

=

0.4675 − 0.46

= 0.75

и

f (0.4675)≈ P(0.4675)= f (x0 )+

y0

q +

0.01

1!

h

+

2 y0

q(q −1)+

3 y0

q(q −1)(q − 2)+

4 y0

q(q −1)(q − 2)(q −3)=19.6133 −

0.6708

0.75 −

2!

3!

4!

1!

−0.09712! 0.75 (−0.25)− 0.00863! 0.75 (−0.25) (−1.25)− 0.00084! 0.75 (−0.25) (−1.25) (− 2.25)=

=19.6133 −0.3354 + 0.0091−0.0067 + 0.0000 =19.2803

сточностью ε =10−4.

x		y		y		2 y		3 y		4 y
0.45		20.1946

49

0.46	19.6133
0.47	18.9425

0.4818.1746

0.49	17.3010
0.50	16.3123
0.51	15.1984

0.5213.9484

0.53

12.5508

0.5410.9937

0.559.2647

0.567.3510

-0.5813
		-0.0895
-0.6708				-0.0076
		-0.0971				-0.0010
-0.7679				-0.0086
		-0.1057				-0.0008
-0.8736				-0.0094
		-0.1151				-0.0007
-0.9887				-0.0101
		-0.1252				-0.0008
-1.1139				-0.0109
		-0.1361				-0.0006
-1.2500				-0.0115
		-0.1476				-0.0004
-1.3976				-0.0119
		-0.1595				-0.0005
-1.5571				-0.0124
		-0.1719				-0.0004
-1.7290				-0.0128
		-0.1847
-1.9137

	Аналогично, для xT	= 0.5410, xn = 0.55, q =	x − xn	=	0.5410 − 0.55	= −0.9. Воспользу-
			h		0.01
	2

емся формулой (2.10.2) для интерполяции назад:

	f (0.5410)≈ P(0.5410)= f (xn )+	yn−1	q +	2 yn−2	q(q +1)+	3 yn−3	q(q +1)(q + 2)+
		1!		2!		3!

+4 4y!n−4 q(q +1)(q + 2)(q + 3)= 9.2647 − 1.72901! (− 0.9)− 0.17192! (− 0.9) − 0.01243! ×

×(− 0.9) 0.1 1.1− 0.00054! (− 0.9) 0.1 1.1 2.1 = 9.2647 +1.5561 + 0.0077 + 0.0002 +

+0.0000 =10.8287

стой же точностью, что и в предыдущей формуле.

Задание № 2. Вычислить приближенное значение функции y = f (x) по интерполяционной формуле Ньютона для интерполяции вперед или назад.

I.

№ варианта		xТ			x				y
1		4.413			4.35				16.3597
6		6.085		4.60				17.7334
11		5.218		4.85				18.7686
16		6.480		5.10				20.0334
21		4.390			5.35				22.2846
26		6.020			5.60				23.5973
					5.85				25.0811
					6.10				26.5278
				6.35				28.3944
				6.60				29.9902

50

II.

№ варианта		xТ		x		y
2		1.210		1.15		66.1659
7		4.680		1.55		63.9989
12		2.000		1.95		61.9658
17		3.870		2.35		60.0551
22		2.411		2.75		58.2558
27		3.275		3.15		56.5583
				3.55		54.6807
				3.95		52.7220
				4.35		50.5229
				4.75		48.1091

III.

	№ варианта		xТ			x				y
	3		0.105		0.10				1.2618
	8		0.121			0.16				1.2764
	13		0.625			0.22				1.2912
	18		0.438			0.28				1.3062
	23		0.534		0.34				1.3213
	28		0.229		0.40				1.3366
					0.46				1.3521
					0.52				1.3677
						0.58				1.3836
IV.						0.64				1.3995
	№ варианта		xТ			x				y
	4		1.871			1.350				4.2556
	9		1.369			1.355				4.3533
	14		1.382			1.360				4.4552
	19		1.394		1.365				4.5618
	24		1.354		1.370				4.6734
	29		1.362		1.375				4.7904
					1.380				4.9131
						1.385				5.0419
						1.390				5.1778
						1.395				5.3202
V.
	№ варианта		xТ			x				y
	5		1.285		1.20				10.6044
	10		1.211		1.21				11.3276
	15		1.272		1.22				11.9671
	20		1.228		1.23				12.5324
	25		1.263			1.24				13.0328
	30		1.237			1.25				13.4776
						1.26				13.8759
					1.27				14.2367
					1.28				14.5688
					1.29				14.8809

51

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

	ORIGIN := 1
	0.35		2.73951
	0.41			2.30080
	0.47			1.96464
x :=			y :=
	0.51		1.78776
	0.56			1.59502
	0.64			1.34310

M 1 := x M 2 := y i := 1...5

M i,3 :=

(M i+1,2

− M i,2 )

xi+1

− xi

i := 1...4

M i,4

:=

(M i+1,3

− M i,3 )

xi+2

− xi

i :=1...3

Mi,5

:=

(Mi+1,4

− Mi,4 )

xi+3

− xi

i :=1...2

M i,6

:=

(M i+1,5

− M i,5 )

xi+4

− xi

M i,7

:=

(M 2,6

− M1,6 )

x6

− x1

Если теперь набрать M = , то на экран выведется точно такая же таблица, что и полученная нами в первом примере этой лабораторной работы.

Получить аналитический вид многочлена можно с помощью меню Symbolics (Символьные выражения) – Simplify (Упростить), набирая по формуле (2.7.1) числовые значения коэффициентов из матрицы M :

P1(z):= 2.73951+ (−7.311833) (z −0.35)

5.298652 − 7.311833 z

P2(z):= 2.73951 + (− 7.311833) (z − 0.35)+14.243056 (z − 0.35) (z − 0.41)

7.342530 −18.136556 z +14.243056 z 2

P3(z):= 2.73951 + (− 7.311833) (z − 0.35)+14.243056 (z − 0.35) (z − 0.41)...

+ (−15.227431) (z − 0.35) (z − 0.41) (z − 0.47)

8.369544 − 25.760930 z + 32.972795 z 2 −15.227431 z3

В Mathcad длинный оператор, не помещающийся в одну строку на странице, можно перенести в следующую строку, нажав [Ctrl][Enter]. Это проделано в предыдущем операторе. Аналогично

P4(z):= 2.73951 + (− 7.311833) (z − 0.35)+14.243056 (z − 0.35) (z − 0.41)...

+(−15.227431) (z − 0.35) (z − 0.41) (z − 0.47)...

+(−102.232684) (z − 0.35) (z − 0.41) (z − 0.47) (z − 0.51)

4.853055 + 7.239956 z − 82.345571 z 2 +162.657283 z3 −102.232594 z 4

52

P5(z):= 2.73951 + (− 7.311833) (z − 0.35)+14.243056 (z − 0.35) (z − 0.41)...

+(−15.227431) (z − 0.35) (z − 0.41) (z − 0.47)...

+(−102.232684) (z − 0.35) (z − 0.41) (z − 0.47) (z − 0.51)...

+825.705728 (z − 0.35) (z − 0.41) (z − 0.47) (z − 0.51) (z − 0.56)

−11.051930 +184.903812 z −870.466826 z 2 +1898.621006 z3 ...

+(- 2001.355769) z 4 + 825.705728 z5

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

xt := 0.552 xx1 := P1(xt) xx2 := P2(xt) xx3 := P3(xt) xx4 := P4(xt) xx5 := P5(xt)

1.262520

⎜1.671068

xx = 1.635251

1.6251511.624499

Эту же задачу можно решить средствами программирования пакета Mathcad. В среде Mathcad потребность в программировании ощущается не только для расширения и совершенствования базового набора математических инструментов, но и в учебных целях. Пакет содержит одиннадцать функций экстраполяции и интерполяции, в частности линейную интерполяцию (linterp) и интерполяцию сплайнами (cspline, pspline, lspline). Полиномиальная интерполяция осуществляется функциями polint и polin2, однако они не являются встроенными и становятся таковыми лишь после подгрузки соответствующего электронного учебника.

Составим поэтому программу, вычисляющую интерполяционный многочлен Ньютона по разделенным разностям. Эту задачу удобно разбить на три части: а) составление таблицы разделенных разностей; б) вычисление коэффициентов функции ω n (x) в формуле (2.7.1) и

в) вычисление самого многочлена Pn (x).

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

на:

Здесь x и y - вектора аргумента и функции, n - число узлов заданной сетки. После обращения к данной подпрограмме

MM := Newt(x, y,6)

53

можно еще раз напечатать полученную матрицу, набрав MM = . Эта таблица будет опять точно такой же, как и в первом примере.

Нахождение коэффициентов функции ω n (x), включая раскрытие скобок и приведе-

ние подобных членов, - задача значительно более трудная. Без сомнения, можно предложить несколько вариантов этой подпрограммы, которые будут длиннее или короче друг друга в зависимости от выбранного алгоритма. Один из вариантов может быть, например, таким:

Подпрограмма omega(x, n) выдает матрицу an×(n+1), в которой по строкам записаны коэффициенты всех функций ω n (x) от первой степени ω 1(x) до n - степени ω n (x), причем коэффициенты расположены начиная со старшей степени аргумента. Последняя подпрограмма PolinNewt реализует формулу (2.7.1) и выдает вектор коэффициентов интерполяци-

онного многочлена Ньютона степени m (1 ≤ m ≤ n −1).

Можно проверить правильность предыдущих вычислений, набрав

		825.705728
		− 2.001360 10
				3
a1: = PolinNewt(x, y,6,5)		1.898621 10	3
	a1 =
		−870.466939		.
		184.903845
		−11.051934

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

54

Структура формулы (2.10.1) позволяет при			вычислении произведений вида
q (q −1)(q − 2)...(q − n +1)		использовать уже готовую	функцию ω n (x). Действительно,

q = x −hx1 = 1h (x − x1 )= 1h ω1 (x), q(q −1)= x −hx1 x −hx2 = h12 ω2 (x) и так далее.

Подпрограмма, реализующая формулу Ньютона для интерполяции вперед, может быть, например, такой:

При написании этой подпрограммы использовались операторы break – выхода из текущего цикла и операции отношения >,<. Операции отношения берутся из панели равенств и отношений (Evaluation and Boolean Palette), которую открывает кнопка .

Эта панель имеет вид

Щелчком по кнопкам этой панели в рабочий документ вставляются соответствующие знаки отношений и т.п.

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

Введем с клавиатуры

	0.01			9.9182
	0.06			9.5193
						M 2 : = Tabl(x, y,6)
x :=	0.11		y :=	9.1365
	0.16			8.7691
				8.4164
	0.21
	0.26			8.0779

0.01		9.9182		− 0.3989		0.0161	− 7 10−4	1.7763568 10−15	2 10	−4
	0.06	9.5193		− 0.3828		0.0154	− 7 10−4	− 2 10−4	0
	0.11	9.1365		− 0.3674		0.0147	−5 10−4	0	0
M 2 =		8.7691		− 0.3527		0.0142	0	0	0
0.16
	0.21	8.4164		− 0.3385		0	0	0	0
	0.26	8.0779		0		0	0	0	0
ε : =10−6 a2 : = Polin Newt Const(x, y,6,5,ε)
		− 0.933333
				3.388
		a2 =	−8.211147
			9.999974

56