Учебное пособие_Численные методы -26 -дек_2014 (97-2003)

Конечные разности различных порядков удобно располагать в форме таблиц двух видов: горизонтальной таблицы разностей (табл. 4.1) или диагональной таблицы разностей (табл. 4.2).

Таблица 4.1	Таблица 4.2
Горизонтальная таблица разностей	Диагональная таблица разностей

	x		y		∆	∆2	∆3
				y		y	y
	x		y		∆	∆	∆
0		0		y0		2y0	3y0
	x		y		∆	∆	∆
1		1		y1		2y1	3y1
	x		y		∆	∆	∆
2		2		y2		2y2	3y2
	.		..		..	..	.
....		...		...		...	....

x	y	∆	∆2	∆3
		y	y	y
x	y
0	0
		∆
		y0
x	y		∆2
1	1		y0
		∆y		∆
		1		3y0
x	y		∆2
2	2		y1
		∆y
		2
x	y
3	3

4.3.Обобщенная степень

Вдальнейшем нам понадобится понятие обобщенной степени. Определение. Обобщенной n-степенью числа х называется

произведение п сомножителей, первый из которых равен х, а каждый следующий на h меньше предыдущего:

х[п]= х (х h) (х 2h) . . . [x (n 1)h],

(4.11)

где h некоторое фиксированное постоянное число.

Показатель обобщенной степени обычно записывается в квадратных скобках. Полагают x[0]=1. При h = 0 обобщенная степень (4.11) совпадает

с обычной х[п]= хп.
Вычислим			конечные разности для обобщенной степени,					полагая
∆x = h. Для первой разности имеем
x	[n]	(x h)	[n]	x	[n]	(x h)x...[ x (n 2)h] x(x	h)...[x (n 1)h]
								[n 1]		,
x(x h)...[x (n 2h)] {(x h) [x (n 1)h]} x(x h)...[x (n 2)h]nh nhx
т.е.
						x[n] nhx[n 1] .			(4.12)
Подсчитываем вторую разность
						101

x		(x		) (nhx	) nh (n 1)hx	[n 2]	nh		(n 1)x	[n 2]	.
2	[n]		[n]	[n 1]				2

Итак,

2	[n]	n(n 1)h	2	x	[n 2]
x

Методом математической индукции легко

.

доказать общую формулу

k

[n]

n(n 1)...[n (k 1)]h

k

x

[n k ]

,

x

где k 1,2,..., n. Очевидно, ∆kх[п]= 0 при k>n.

Из формулы (4.12) вытекает также простая формула конечного

суммирования.

Пусть x0 , x1, x2, ... равноотстоящие точки с шагом h:

xi+1 – xi=h(i = 0, 1, 2 ...) .

N 1

Рассмотрим сумму SN

xi[n] . Так как в силу формулы (4.12) имеем

i 0

[n 1]

N 1

n 1

n 1

x

x0

x[n]

, то

S N xi [n]

xN

.

h(n

1)

i 0

h(n 1)

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

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

Простейшая задача интерполирования заключается в следующем. На отрезке [а,b] заданы n+1 точки х0 , х1, ... , хn , которые называются узлами интерполяции, и значения некоторой функции f(х) в этих точках

f(x0)=y0, f(x1)=y1, ... ,f(xn)=yn. (4.13)

Требуется построить функцию F(х) (интерполирующая функция), принадлежащую известному классу и принимающую в узлах интерполяции те же значения, что и f(х), т. е. такую, что

F(x0)=y0, F(x1)=y1, ..., F(xn)=yn. (4.14)

Геометрически это обозначает, что нужно найти кривую у = F(х) некоторого определенного типа, проходящую через заданную систему

точек Mi(xi,yi) (i = 0, 1, 2, ...) (рис. 4.1).

В такой общей постановке задача может иметь бесчисленное множество решений или совсем их не иметь. Однако эта задача

102

становится однозначной, если вместо произвольной функции F(х) искать полином Рп(х) степени не выше п, удовлетворяющий условиям (4.14), т. е. такой, что

Pn(x0)=y0, Pn (x1)=y1, ... ,Pn(xn)=yn.

Полученную интерполяционную формулу y=F(x) обычно используют для приближенного вычисления значений данной функции f(х) для значений аргумента х, отличных от узлов интерполирования. Такая операция называется интерполированием функции f(х). При этом различают интерполирование в узком смысле, когда x x0 , xn , т. е.

значение х является промежуточным между х0		и хп, и
экстраполирование, когда	x x0 , xn . В дальнейшем под	термином

интерполирование мы будем понимать как первую, так и вторую операции.

4.5. Первая интерполяционная формула Ньютона

Пусть для функции y=f(x) заданы значения yi = f(xi) для равноотстоящих значений независимой переменной: хi = х0 + ih, i= (0, 1, 2,. . . , п), где h шаг интерполяции. Требуется подобрать полином Рn(х) степени не выше п, принимающий в точках xi значения

Pn(xi)=yi, i= (0, 1, 2, . . . ,

Условия (4.15) эквивалентны тому, что	Pn
	m

п).

x		y
	m
0			0

при

m

(4.15)

0,1,2,...n .

Следуя Ньютону, будем искать полином в виде

Pn x a0	a1 x x0 a2 x x0 x x1 ... an x x0 x x1 ... x xn 1				.	(16)
Пользуясь обобщенной степенью, выражение (4.15) запишем так:
Pn	x a0 a1 x x0	a2 x x0	... an x x0	.		(4.16')
	[1]	[2]	[n]
Наша задача состоит в определении коэффициентов ai (i= 0, 1,					2,	. . . , п)

полинома Рn(х). Полагая х = х0

в выражении (4.16'), получим Pn(x0)=y0=a0.

Чтобы найти коэффициент a1

,составим первую конечную разность

[1]

h 3a3 x

[2]

h ...

[n 1]

h .

Pn x a1h 2a2 x x0

x0

nan x x0

Полагая в последнем выражении х = х0,

получим

Pn x0 y0

a1h , откуда

a1

y

0

. Для определения коэффициента

a2 составим конечную разность

1!h

второго порядка

2 Pn x 2!h2a2 2 3h2a3 x x0 [1] ... n 1 nh2an x x0 [n 2] .

Положив х = x0 , получим 2 P

x

2 y

2!h2a

, откуда a

2 y0 .

n

0

0

2

2

2!h2

103

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

что

		y
a		i
			0
i		i!h	i

i

0,1,2,...,

n

, где положено

0! 1

и

0 y

y

.

Подставляя найденные значения коэффициентов ai в выражение

(4.16'), получим интерполяционный полином Ньютона

Pn x

			y				2					n
y				x x	[1]		y		x x	[2]	...	y
			0					0					0
	0		1!h	0			2!h	2	0			n!h	n

x

x
	[n]
0

.

(4.17)

Легко видеть, что полином (4.17) полностью удовлетворяет требованиям поставленной задачи. Действительно, во-первых, степень полинома Рп(х) не выше n, во-вторых, Pn x0 y0 и

P	x	k	y	0
n
y			k y
	0			0

y	0	x
h			k
k k 1
	2!

			y
x			2					x		x
						0
	0		2!h		2				k		0
							k 1 ...1
2		...		k
y
	0							k!

x	k	x
			1
y			1
k
		0

			y
				k		x
	...				0
			k!h		k
	k	y		y			,
			0			k

k	x	0	x	k	x	... x	k	x	k 1
					1

где k 1,2,...n .

Заметим, что при h→0 формула (4.17) превращается в полином

Тейлора для функции

того, очевидно,	lim x
	h 0

у. В самом

x		[n]	x x
	0
			0

			k
деле,		lim	h
		h 0
	. Отсюда
n

y		d	k	y
k	0			k
		dx			x x
					0

при h→0

y	k	x0	. Кроме

формула (4.17)

принимает вид полинома Тейлора:

P	x y(x	) y (x	) x x	...
n	0	0	0

y	(n)	(x	)
		0	x
	n!

x
	n
0

.

Для практического использования интерполяционную формулу Ньютона (4.17) обычно записывают в несколько преобразованном виде.

Для этого введем новую переменную q по формуле

q

xx0 h

, тогда

(x x

)

[i]

(x x

)

(x x

h)

[x x

(i 1)h]

0

0

0

...

0

q(q 1)(q 2)...(q i 1) (i 1,2,..., n).

h

i

h

h

h

где

Подставляя эти выражения в формулу (4.17), получим

P(x) y0

q y0

q(q 1)

2

...

q(q 1)...(q n 1)

n

,

(4.18)

2!

y0

n!

y0

q

x x

0

представляет

собой

число шагов, необходимых для

h

достижения точки x, исходя из точки х0. Это и есть окончательный вид

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

Формулу (4.18) выгодно использовать для интерполирования функции у=f(х) в окрестности начального значения x0 , где q мало по абсолютной величине.

Если в формуле (4.18) положить n=1, то получим формулу линейного интерполирования

P1 (x) y0 q y0 .

104

При n=2 будем иметь формулу параболического или квадратичного интерполирования

P2 (x)

y	0

q y	0

q(q 1) 2

y
2
	0

.

Если дана неограниченная таблица значений функции у, то число п в интерполяционной формуле (4.18) может быть любым. Практически в этом случае число п выбирают так, чтобы разность ∆nуi была постоянной с заданной степенью точности. За начальное значение х0

можно принимать любое табличное значение аргумента х.

Если таблица значений функции конечна, то число п ограничено, а именно: п не может быть больше числа значений функции у, уменьшенного на единицу.

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

4.6. Вторая интерполяционная формула Ньютона

Первая интерполяционная формула Ньютона практически неудобна для интерполирования функции вблизи конца таблицы. В этом случае обычно применяется вторая интерполяционная формула Ньютона.

Выводом этой формулы мы и займемся.

Пусть имеем систему значений функции yi y(xi ) (i 0,1,2,..., n) для

равноотстоящих значений аргумента

x	i

x	0

ih

.

Построим
P	x a
n

интерполирующий
0	a	x x	a	x x
	1	n	2	n

полином
x x	n 1	...

следующего вида:
a	x x	n	x x	n 1	... x x
n					1

или, используя обобщенную степень, получаем

Pn x a0 a1 x xn [1] a2 x xn 1 [2] a3 x xn 2 [3] ... an x x1 [n] .

Наша задача состоит в определении коэффициентов а0, а1 таким образом, чтобы были выполнены равенства Pn xi yi (i

(4.19) , a2, .... аn

0,1,2,..., n) .

Для этого необходимо и достаточно, чтобы

Pn xn i	yn i (i 0,1,..., n) .		(4.20)
i	i
Положим х = хn в формуле (4.19). Тогда будем иметь		Pn xn yn	a0 .

Следовательно, a0 yn .

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

Pn x a1 1h a2 2h x xn 1 [1] a3 3h x xn 2 [2] ... an nh x x1 [n 1] .

Отсюда, полагая х = хn-1 и учитывая соотношения (4.20), будем иметь

105

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

a
1

y	n 1
h

P	x	n 1	y	n 1
n

.

a1h

.

a	2

Аналогично составив вторую разность от

2

x a

2!h

2

a

3 2h

2

x

x

[1]

...

P

2

n 2

Полагая х

n

3

= хп-2,

находим

Pn xn 2 yn

2

2

y

2

.

n 2

2!h

2

Pn(х), получим

an n n 1 h				2	x x1	[n 2]	.
2	a2	2!h	2	и, таким образом,

Характер закономерности коэффициентов аi достаточно ясен. Применяя метод математической индукции, можно строго доказать, что

		y
a		i
			n 1
i		i!h		i

i

0,1,2,...,

n

.

(4.21)

Подставляя эти значения в формулу (4.19), будем иметь окончательно

P x y				y	n 1	x
n		n		1!h
Формула (4.22)
Ньютона.
Введем				более
тогда	x x			x x			h
	n 1				n
	h				h

		y				y
xn		2			x xn x xn 1 ...	n		x xn	... x x1 .	(4.22)
			n 2				0
		2!h		2		n!h	n

носит название второй интерполяционной формулы

удобную			запись формулы (4.22). Пусть				q	x x	n	,
								h
q 1	,	x x		n 2	q 2	и т.д.
			h

Подставив эти значения в формулу (4.22), получим:

(4.22’)

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

y=Pn(x).

Как первая, так и вторая интерполяционные формулы Ньютона могут быть использованы для экстраполирования функции, т. е. для нахождения значений функции y для значений аргументов x, лежащих вне пределов таблицы. Если x<x0 и x близко к x0, то выгодно применять первую интерполяционную формулу Ньютона, причем тогда

q	x x	0	0 .
	h

Если же x>xn и x близко к xn, то удобнее пользоваться второй

интерполяционной формулой Ньютона, причем

q x xn 0 . h

Таким образом, первая интерполяционная формула Ньютона обычно используется для интерполирования вперед и экстраполирования назад, а

106

вторая интерполяционная формула Ньютона, наоборот, для

интерполирования назад и экстраполирования вперед.

4.7. Интерполяционные формулы Гаусса

Пусть имеется 2n+1 равноотстоящих узлов

x

n

, x

(n 1)

,..., x

1

, x

, x

,..., x

n

0

1

где xi

xi 1 xi

h const

(i = - n, -(n-1),…,n - 1),

интерполирования

1 , xn ,

и для функции y=f(x) известны ее значения в этих узлах yi=f(xi)

(i=0, ±1,…, ±n).

Требуется построить полином P(x) степени не выше 2n такой, что

P(xi)=yi при i=0, ±1,…, ±n.

Из последнего условия вытекает, что для всех соответствующих

значений i и k
k P(x ) k y .		(4.23)
i	i

Будем искать этот полином в виде

P(x) a				0	a		(x x					0	) a		2	(x x					0	)(x x				) a		3	(x x		1		)(x x			0	)(x x		) a		4	(x x	1	)(x x	0	)(x x	)
						1																			1													1								1
(x x		2	) a			5	(x x						2	)(x x				1		)(x x				0	)(x x )(x x					2		) ...
																											1
... a	2n 2			(x x				(n 1)			)...(x x						1		)(x x				0	)(x x			)...(x x					n 1		) a	2n		(x x	(n 1)		)...(x x			1	)(x x	0	)(x x	)...
																										1																				1
...( x x		n 1			)(x x				n	).

(4.24)

Вводя обобщенные степени, получим

P(x) a			a (x x					)	[1]	a			(x x			)		[2]	a			(x x			)	[3]	a		(x x		)	[4]	...
		0					0				2				0									1				4		1
			1																		3
... a		(x x			)	[2n 1]				a				(x x						)	[2n]		.
	2n 1			(n 1)								2n					(n 1)

(4.25)

Применяя для коэффициентов ai (i=0,1,…,2n) тот же способ, что и при выводе интерполяционных формул Ньютона, и, учитывая формулу (4.23), последовательно находим:

							y						y						y						y						y								y
													2						3						4					2n 1								2n
a		y		,	a			0	,	a				1	,	a				1	,	a				2	,..., a					(n 1)		,	a					n		.
	0		0								2			2			3			3			4			4		2n 1					2n 1			2n					2n
					1		1!h						2!h						3!h						4!h					(2n 1)!h								(2n)!h

Далее, введя новую переменную

q	x x	0
	h

и сделав соответствующую

замену в формуле (4.25), получим первую интерполяционную формулу Гаусса

q(q 1)

2

(q 1)q(q 1)

3

P(x) y0

q y0

y 1

y 1

2!

3!

(q 2)(q 1)q(q 1)(q 2)

5

y 2

...

5!

(q n 1)...(q n 1)

2n 1

(q n 1)...(q n)

...

y (n 1)

(2n 1)!

(2n)!

(q 1)q(q 1)(q 2) 4 y 2 4!

(4.26)

2n y n ,

или короче

107

q

[2]

(q n 1)

[2n 1]

2n 1

P(x) y

q y

2

...

y

0

0

2!

y

1

(2n 1)!

(n 1)

где x = x0 + qh		и	q[m] = q(q - 1)…[q– (m-1)].
Первая	интерполяционная				формула Гаусса
разности y0	,	y 1 , y 1 ,		y 2 , y 2 ,		y 3 ,...
		2	3	4	5	6

(q n 1)	[2n]
(2n)!

содержит

	y n , (4.26`)
2n

центральные

Аналогично

Гаусса,	y
разности y 1 ,
	2

можно

	,	3		,
1		y	2

получить					вторую интерполяционную формулу
	содержащую							центральные
4		5		,	6		,...
y	2	, y	3		y	3

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

P(x) y

q y

(q 1)q

2

(q 1)q(q 1)

3

0

1

2!

y

1

3!

y

2

(q 2)(q 1)q(q 1)

4

......

(q n 1)...(q n 1)

2n 1

y

4!

y

2

(2n 1)!

n

(q n)(q n 1)...(q n 1)

2n

y

,

(2n)!

n

(4.27)

или в сокращенных обозначениях

(q 1)

[2]

(q 1)

[3]

(q 2)

[4]

P(x) y

q y

2

3

4

...

0

1

2!

y

1

3!

y

2

4!

y

2

(q n 1)

[2n 1]

2n 1

(q n)

[2n]

...

y

2n

y

,

(2n 1)!

n

(2n)!

n

(4.27`)

где x= x0 +qh.

4.8. Интерполяционная формула Лагранжа

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

интерполяционной формулой Лагранжа (рис. 4.2 а – рис. 4.2 б).

Пусть на отрезке [а,b] даны п+1 различных значений аргумента x0, x1, x2, .... xn и известны для функции у = f(x) соответствующие значения:

f(x0)=y0, f(x1)=y1,…f(xn)=yn.

Требуется построить полином Ln(x) степени не выше п, имеющий в заданных узлах x0, x1, .... хn те же значения, что и функция f(x), т. е. такой, что

Ln (xi)=yi(i = 0, 1, 2, ..., n).

108

а)

б)

4.2. Интерполяция методом Лагранжа (а – общая задача, б – частная задача)

Решим сначала частную задачу: построим полином pi(x) такой, что

pi(xj)=0 при i≠j и pi(xi )= 1 (рис. 4.2 б).

Другими словами, эти условия можно записать следующим образом:

			1, если	j i,
pi	(x j	) δij
				j i,
			0, если

(4.28)

где δij символ Кронекера.

Так как искомый полином обращается в нуль в п точках х0 , х1, ... , хi-1 , хi+1, ... , хn, то он имеет вид

pi (x) Сi (x x0 )(x x1 )...(x xi 1 )(x xi 1 )...(x xn ) ,

(4.29)

где Ci постоянный коэффициент. Полагая x = xi в формуле (4.29) и учитывая, что рi(xi)= 1, получим

	Сi (xi x0 )(xi				x1 )...(xi xi 1 )(xi xi 1 )...(xi			xn ) 1.
	Отсюда находим
		Ci			1				.
				(xi x0 )(xi x1 )...(xi xi 1 )(xi xi 1 )...(xi xn )
	Подставив это значение в формулу (4.29), будем иметь
	pi (x)	(x x0 )(x x1 )...(x xi 1 )(x xi 1 )...(x xn )					.	(4.30)
		(xi	x0 )(xi		x1 )...(xi xi 1 )(xi xi 1 )...(xi xn )

Теперь перейдем к решению полинома Ln(x), удовлетворяющего Этот полином имеет вид

общей задачи, т.е. к отысканию указанным выше условиям Ln(x)= yi.

n		n
	(x)
L
		i 0

p	(x) y	i
i

.

(4.31)

В самом деле, во-первых, очевидно, степень построенного полинома Ln(x) не выше п и, во-вторых, в силу условия (4.28) имеем

109

содержит явно yi, что

получим

полинома

(4.35’)

(4.33)

(4.34)

n
Ln (x) pi (x j ) yi	p j (x j ) y j	y j ( j 0,1,..., n) .
i 0

Подставив в формулу (4.31) значение pi(x) из (4.30), получим

n

(x x

)(x x )...(x

x

)(x x

)...(x x

)

Ln (x) yi

0

i 1

n

.

1

i 1

(x

x

)(x

x )...(x

x

)(x

x

)...(x

x

)

i 0

0

i

i

n

i

i

1

i 1

i 1

i

(4.32)

Это и есть интерполяционная формула

Докажем единственность полинома

Предположим противное. Пусть

Ln (x)

Лагранжа.

Лагранжа.

полином, отличный от Ln(x),

степени не выше п, и такой, что Ln	(xi ) yi . (i = 0, 1, 2, ...). Тогда полином
Qn Ln (xi ) Ln (x) , степень которого,	очевидно, не выше n, обращается в

нуль в п+1 точках х0 , х1, ... , хn, то есть Qn≡0. Следовательно, Ln (xi ) Ln (x) . Отсюда, в частности, следует, что если узлы интерполирования

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

Формуле (4.32) Лагранжа можно придать более сжатый вид. Для этого введем обозначение

n 1 (x) (x x0 )(x x1 )

...(x xn ) .

Дифференцируя по х это произведение, получим:

		n
	(x)		(x x	0	)(x x )...( x x	j 1	)(x x	j 1	)...(x
n 1					1
		j 0

x	n	)

.

Полагая х =хi (i = 0, 1, 2, ...), будем иметь

n

(x)

(x

x

)(x

x )...( x x

j 1

)(x

i

x

j 1

)...(x

n 1

i

0

i

1

i

i

j 0

Внося выражения (4.33) и (4.34) в формулу (4.32),

n

yi

Ln (x) n 1 (x)

.

n 1 (xi )(x

xi )

i 0

x	n	)

.

Следует отметить, что формула Лагранжа бывает иногда важно.

Рассмотрим два частных случая интерполяционного Лагранжа.

При п=1 мы имеем две точки, и формула Лагранжа представляет в этом случае уравнение прямой у = L1 (x), проходящей через две заданные точки:

y	x b	y			x a	y ,
	a b		0		b a 1

где а, b абсциссы этих точек.

110