book1989

явных двухшаговых итерационных методов:
В- (**+1 — **) + (1 —ak+1) (хк — xk-1)	Ас* = /, k = 1 , 2 , ... (26)

Здесь а„+1, тЬч— итерационные параметры, которые будут опреде­ лены далее. Для начала счета необходимо задать два начальных приближения х0 и х,. Начальное приближение х„ будем задавать произвольно, а вектор х, вычислять по одношаговой формуле, ко­

торая получается из (26) при k = 0, a ^ l ,	т. е. по формуле
+ Лх0 = / .	(27)
Ti

Если параметры ak+i, тЛ+1 найдены, то новое приближение хк+1

выражается через д	е э предыдущих значения хк и xh- t по формуле
x k+ t ~	a k+ i X k -f- (1	ctfc+i)x k- i Tft+ia k + i w k ,		(28)
где wh = B~lrh, гк-=Ахк—f.
5. Минимизация	погрешности. Перейдем к		вопросу о выборе
итерационных параметров в		методе (26). Для	погрешности	zk=
= хк—х получим уравнения
zk+i = ah+i(E—Tk+lB -,A)zh-lr ( l —ah+i)zh-i,			&=1, 2, ...,
	z,= (£—■xiB~1A)za.

Введем, как и ранее, вспомогательную функцию vk = A'/2zk, для которой |1щ|| = ||24||л. Функция vk удовлетворяет уравнениям

Щ+1 =	<Xk+i {Е — xk+1C) vk + (1 — a.k+i) vk-u k — \,2, . . .,	(29)
Vi =	(E — TjC) v0,	(30)

где C= A,/2B~lA,/2. Будем считать, что А и В — симметричные по­ ложительно определенные матрицы, удовлетворяющие неравен­ ствам

YiSs^ А ^ « ( 2В,	> 7*> 0.	(31)
Тогда С* = С>0, причем		(32)
"fiE^C^-y^E.
Исключая последовательно из уравнений (29),		(30) векторы
vu v2, , щ_1; найдем, что		(33)
vk = Ph(C)v0,
где Рк{С) — многочлен степени k	от оператора С,	удовлетворяю­

щий условию ДДО) =Е.

Поставим задачу выбрать итерационные параметры тк, ак так, чтобы при любом /1= 1, 2, ... была бы минимальной ||у„|| = ||z„IL. Обратим внимание на отличие от постановки задачи, возникающей при построении оптимального чебышевского набора итерационных параметров (см. § 6). Там при фиксированном п требовалось най­ ти параметры, минимизирующие ||2„||А. Теперь же требуется боль­ шее— минимизировать ||zJ!Aпри каждом п.

121

Параметр TI находится из условия минимума ||щ|| при заданном векторе у0. Так же, как и в методе скорейшего спуска, получаем

	т ~ ? Г(Соь, ув)		(34)
	1	ПС«оР

Отметим, что при таком выборе tj выполняется равенство (Cvlt у„) = = 0, т. е. векторы и, и и0 ортогональны в смысле скалярного про­ изведения

(и, v)c~ (Си, v).

Далее, рассмотрим погрешность vk=Pk(C)v0,

возникающую на k-й итерации, и запишем многочлен Рк(С) в виде

Pk(C)=E + 2 арС\	(35)
1=1

где ар — числовые коэффициенты, определяемые параметрами а,-,

Ti, t = l, 2, ... , k. Тогда

k

Vk = V0+ 2 aik)Clv0, * = 1 , 2 , ... (36) i=l

Найдем условия, которым должны удовлетворять коэффициен­ ты а!"), минимизирующие ||п„||2. Из (36) получим

I VnIP =

2

f lW

С'и») + 2

i

^

( v °> C'vo) + I

г .

(37)

i,/=i

/=1

т. е. Цу„||2 является многочленом второй степени по переменным

... ,а р .

Приравнивая

к нулю

частные производные

<3||ц„||г/

/дар,

j= 1, 2 , ..., п, приходим к системе уравнений

2

а‘п) (СЧ , СЧ ) + (C’v0, Va) = О,

(38)

L= 1

решение которой а{р,

/ = 1, 2 , ..., я,

и обращает в минимум ||цп||2.

6.

Выбор итерационных параметров в методе сопряженных гра­

диентов. Целью дальнейших рассуждений является нахождение па­

раметров ак,

т„, k = l,

2, ..., п, для

которых выполнены

условия

(38). Заметим прежде всего, что (38)

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

(

С

Ч

vn),

= 0,

/=

1, 2, ..., п.

(39)

Л е м м а

1. Условия (39) эквивалентны условиям

(Cvjt Vn)= 0,

/ = 0,

1,

...,

л—1.

(40)

Д о к а з а т е л ь с т в о .

Согласно (36) имеем

Си/ = Cv0-f 2 a\nCinv0,

122

поэтому
(Cvh vn) = (Cv0, vn) -f 2	a\n (C‘+1v0, vn) =
i=l
	= (Cv0, V n) -f 2 a<i-x (c 4 . Vn)• (41)
		1=2
Пусть выполнены условия (39). Тогда,		если /+ 1 ^ п (т. е.
^ п—■1 ), то
(Cv0, vn) = 0,	(C2va, vn) =0, ...,	(Ci+lv0, vn) = 0.
Поэтому из (41) при j ^ n —1 получим
	(Cvj, vn) = 0.

Итак, условия (40) следуют из (39). Покажем, что верно и об­ ратное, т. е. из (40) следует (39). Доказательство проведем индук­

цией по числу /. Условие (40)				при /'= 0 совпадает с условием			(39)
при	/=1. Предположим,		что	условия (39)	выполнены	для	/ =
= 1 , 2 , ..., k, и покажем, что				они выполнены и для / = fe+ 1 , где
k ^ n —1.
Из (40) при j = k получим, учитывая (36),
О =	(Су*, vn) = ( Cv0+	J	aVC^X, vn
		i= l
				= {Cvo, Vn) +	(c 4> Vn)-		(42)
					1=2
По	предположению	индукции		условия (39)	выполнены	при	/ =
= 1, 2, .... А. Поэтому из (42) получим

(с*+х, vn) = о.

Поскольку ар ф 0 (так как Рк(С) — многочлен степени k), отсюда получаем (Ch+lv„, у„)=0, т . е. условия (39) выполнены и при /' =

= А + 1.

Лемма 1 доказана. Она потребуется для построения оптималь­ ных итерационных параметров в методе (26).

Заметим, что число п в лемме 1 предполагалось фиксирован­ ным, в то время как при постановке задачи оптимизации мы тре­ бовали, чтобы ||ип|| была минимальной при любом п= 1, 2, ... По­ этому оптимальные параметры надо отыскивать не из условий (40) при фиксированном п, а из условий

(Cvu vn)= 0, п= 1,2, ... , / = 0, 1, ..., п—1.

(43)

Если такие параметры будут найдены, то это будет означать, что

построена ортогональная	(в смысле		скалярного произведения
{и, v)c=(Cti, v )) система	векторов	v0,	vu .... vh, ...	Поскольку
пространство решений системы (1 )		имеет размерность		т, постро­
				123

енная ортогональная система будет содержать не более т векто­ ров. Это означает, что начиная с некоторого п (п^.т) погрешности vn обратятся в нуль, т. е. метод сойдется за конечное число ите­

раций.

Перейдем к построению итерационных параметров, для которых выполнены условия (43). Параметры а, и т, найдены согласно (34):

Ctj — 1 , Tj — (Суо, tip)	(44)
II Ci>„||2

Пусть параметры щ, т2, ■• •, т*, а,, а2, ■■■, ак уже выбраны опти­ мальным образом. Тогда согласно (43) выполняются условия

(Cvj,Vi)= 0, i= 1, 2, . . . , А, / = 0, 1, ..., i—1.

(45)

Построим оптимальные параметры т*+1, аЛ+). Согласно лемме 1

при n=k+ 1 должны выполняться условия
{Cvh t'fi+1)=0,		j = 0,		1, ..., k.		(46)
Часть из этих условий, а именно условия (46) при / = О, I, . . . ,						k—2,
следует из (45) . Действительно, согласно (29) имеем
(vk+1, Cvj) = а*+1 (vk, Cv,) — сснЛп (Cvk, Cv,) -f (1 — ct*+1) {vk-u						Cv,).
Из (45) при i= k и i = k—1 получим, соответственно,
(Cv;, vk) = Q,	/ =	0,	1,	....	k — 1,
(Cv,,vk_l) = 0,	/ =	0,	1,	. . .,	k — 2.
Поэтому						(47)
(v/l+l, CVj)=—al!+iTl4-1(Cvk, Cv,)

при / = 0, 1, .. ., k—2 .

Покажем, что для этих же значений / выполняются равенства (Cvh, Cvj)=0. Запишем уравнение (29) при k = j:

i>+1= «и-1{Е—т,+,С) иj+ (1—ocj+,) Uj.,

и найдем отсюда

CVj = — —	V j --------- -------- [Vjrl — (1 — К/+1)
т/+i	ai+iTi+i

Умножая предыдущее соотношение скалярно па Cvh и учитывая симметричность матрицы С, получим

(Cvh Cvk) =

—- —— (vit Cvk) --------	!-----		[(иу-+1, Cvk) — (1 — ot/+i) (f/-i, Cvk)\ =
T/+i	a/+iT/+i
= — — (Си/,	yA) -----------!------[(C u/n, Vk)				(1 a, rl) (Cw/.J, Уд.)].
T/+1		K/T1T/+1
Из (45) при i= k имеем
(Сиу, uft) =		0,	/ = 0, 1	, . . .,	k — 1 ,
{Cvj+1, vk) =		0,	j = 0, 1	, . . .,	k — 2,
{Cvhu	vk) =	0,	/ = 1 , 2	, . ..,	k.
124

Следовательно			(Сщ, Cvh) =0 при / = 0, 1, ..., k—2, и согласно			(47)
имеем	(vh+i, Cvj) = 0, / —0, 1, ..., k—2.
Итак, из всех условий (46) остаются лишь два:
			(Со*-,, t’kki) =		0,	(48)
				(Си*, и*+1) =	0.	(49)
Подставляя в (48) значение vh+l из (29), получим
0	= а*+1 (vk, Cvk-x) — rJ-kцТ*н (Со*, Си*-,) + (1 — a*+1) (о*_ь Co*-,).
Согласно		(45)	при i= k, j = k—1 имеем (Сщ_ь vh)=0, так что пре­
дыдущее уравнение принимает вид
		1	Tfe-f.i(Со*,	Cyfi-i) -r(l	a*+1) (Сщ_ь I'^-j) =0.	(50)
Далее,		подставляя в		(49) значение щ+1 из (29), получим
	0 =	a*+1 (vk, Cvk) — а*+1т*+, (Со*, Со*) -f- (1 —■а*+1) (о*-,, Со*).

Последнее слагаемое в этом тождестве равно нулю, так как соглас­ но (45) при i = k, j = k—1 имеем

{vh- u	Cvk) = (Сщ_,,			vk) =0.
Таким образом, приходим к тождеству
aA+i[(Cu/,, vk)		Tft+i (Co*,		Co*) ] =0,
из которого находим значение параметра хЛ+1:
	Ti+l	(Cvk■vk)			(51)
		II Со*||*
Обратимся теперь к уравнению (50)				и исключим из него выра­
жение (Со*, Со*-,). Из уравнения (29) получим
Со*_, = — о*_, — ------[vk— (1 — a*) и*_2],					(52)
следовательно,		klk
(Со*, Со*_,) = — (Со*, vk-!)-----[(Си*, vk) — (1 — a*) (Cvk, vk-2)].
4	ak4
Согласно (45) имеем
(Co*, vk-i) =		(o*, Cvk-J = 0,
(Co*, o*_2) =		(vk, Ськ-2) = 0.
Поэтому
(Си*, Co*_,) =		-------- (Co*, vk).
			anxk
Подставляя это выражение в (50), получим
afe+iTfe*i	(с^> vt)		1	— a*+1 = 0.
akxk	(Cvk-i- vk-i)
					125

Отсюда приходим к рекуррентной формуле для параметров оц+1:

TA+l	1	(Cvk, vk)
—	'«A	k = 1, 2, . . ., otj = 1. (53)
Ч		(Cvk-V Vk~i)

Формулами (51), (53) задаются выражения для итерационных параметров в методе сопряженных градиентов. Скалярные произ­ ведения, входящие в эти выражения, вычисляются в процессе ите­ раций. Учитывая, что

vk= Av,zk,	С = AV'B~1AV‘, zk= xk	— x, Azk — Axk — f = r k,B-1rk= w k,
получим	Cvk= A l/2wk, (CVb, vh) =	(wk, rk), (Cvh, Cvh) = (Awh, w„).

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

:	к . rk)		6 = 0, 1 , ....			(54)
	(Awk, wk)
®*+i — 1		4+1 1	(®A- 4 )	k = l , 2	,	ax= 1 . (55)

4ak (Awk~v wk~i)

7.Оценка погрешности в методе сопряженных градиентов. Выше отмечалось, что в методе сопряженных градиентов точное ре­ шение системы уравнений ( 1) получается за конечное число итера­ ций, равное порядку системы. Если порядок системы велик, то мо­ жет оказаться полезной и оценка погрешности. Эта оценка не хуже, чем в одношаговом итерационном методе с чебышевским набором параметров. Действительно, из выражения для погрешности (33) получаем

llvJ = llxn-xllA^HPn(C) Hllxa-xllA.

Поскольку Рп(С) — многочлен степени п от оператора С, удо­ влетворяющий условию Р„(0)=.Е, выполняется оценка

2р"	1 - К б
\Рп(С)ШТп{С)1	PI — 1+ / I
1 + рГ

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

["Ь	Та]. 7\.(0) = 1-
	Таким образом, для погрешности метода сопряженных градиен­
тов справедлива оценка
	\\хп—* IIA sS< 7 * 1 1 - 4 — я |1а .
где	g„= 2p"/(l +pjn).
126

Г Л А В А 3

ИНТЕРПОЛИРОВАНИЕ И ПРИБЛИЖЕНИЕ ФУНКЦИЙ

В настоящей главе излагаются вычислительные аспекты некото­ рых задач теории приближения функций. Задача интерполирования состоит в том, чтобы по значениям функции f(x) в нескольких точ­ ках отрезка восстановить ее значения в остальных точках этого отрезка. Разумеется, такая задача допускает сколь угодно много решений. Задача интерполирования возникает, например, в том случае, когда известны результаты измерения yh= f{ x k) некоторой физической величины f(x) в точках хк, &= 0, 1, ..., п, и требуется определить ее значения в других точках. Интерполирование используется также при сгущении таблиц, когда вычисление зна­ чений f\x) трудоемко. Иногда возникает необходимость прибли­ женной замены или аппроксимации данной функции другими функциями, которые легче вычислить. В частности, рассматривает­ ся задача о наилучшем приближении в нормированном простран­ стве Я, когда заданную функцию f e / / требуется заменить линей­ ной комбинацией ср заданных элементов из Я так, чтобы отклоне­ ние ||/—ср|| было минимальным. Результаты и методы теории ин­ терполирования и приближения функций нашли широкое примене­ ние в численном анализе, например при выводе формул численного дифференцирования и интегрирования, при построении сеточных аналогов задач математической физики.

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

1.Интерполяционная формула Лангранжа. Пусть на отрезке

ai^Lx^Lb заданы точки xh, £= 0, 1 , ..., п (узлы интерполирования), в которых известны значения функции f(x). Задача интерполиро­ вания алгебраическими многочленами состоит в том, чтобы по­

строить многочлен

Ьп(х)==а0+ а ^ + . .. + апхп

(1)

степени п, значения которого в заданных точках xk, k = 0, 1, ..., п, совпадают со значениями функции f(x) в этих точках.

Для любой непрерывной функции f(x) сформулированная зада­ ча имеет единственное решение. Действительно, для отыскания ко­ эффициентов а„, «1, ■■•, ап получаем систему линейных уравнений

aQ+ a^i + а2х} + . . . + апх?= / (*,-), i = 0, 1 , .. ., п,

(2)

определитель которой (определитель Вандермонда [12, с. 33]) от­ личен от нуля, если среди точек хи i = 0, 1, ..., п, нет совпадающих.

Многочлен Ln{x), удовлетворяющий условиям
Ln{Xi)=f(Xi), i=0,	п,	(3)
называется интерполяционным многочленом	для	функции f(x),
построенным по узлам {*;}?.
		127

Решение системы (2) можно записать различным образом. Наи­ более употребительна запись интерполяционного многочлена в форме Лагранжа и в форме Ньютона.

Интерполяционная формула Лагранжа позволяет представить многочлен Ln(x) в виде линейной комбинации

П

и (*) = 2 с><м / (**)	(4)
k=Q

значений функции f(x) в узлах интерполирования.

Найдем явное выражение для коэффициентов ск(х). Из условий

интерполирования (3) получаем
П
2 ck ( X i ) f {Xu) = f ( X i ) ,	i = 0, 1, . . ., n.
£=0

Эти соотношения будут выполнены, если на функции ск{х) нало­ жить условия

Ck (Xi) =		0,		k,
		1 .		■k,	i = 0, 1 ,	п,
которые означают, что каждая				из	функций	ck{x), k = 0, 1 , ..., п,
имеет не менее п нулей на			[а,	Ь].	Поскольку L„(x) — многочлен
степени п, коэффициенты ск(х)				естественно		искать также	в виде
многочленов степени п, а именно в виде
ск(х) =Кк(х—х0) (х—х1)...(х —хк- 1) (х—хк+1)...(х —хп).
Из условия с„(хь) = 1 находим
—{Хи	Х(^ (х'и Хк)		. . . (Хи		■Xu—i)(Xu	Xuui) • ■- (Хи	Хп).
Таким образом, коэффициенты ск{х) интерполяционного много­
члена (4) находятся по формулам
		Ck (х) =		П
				— --------- .			(5)
				П (xk ~ xi)
				/Ф&
Часто коэффициенты ск(х) записывают в другом виде. Введем
многочлен oi(x)	степени п+ 1 :

со(х) = (х—х0) (x—xl)...(x —xk- i) (х—хк) (х—хк+1)...(х —хп) (6)

и вычислим его производную в точке хк:

со'(хк) = (хк—х0)... {хк—хк^ ) (xh—xh+i) ... (хк—хп) .

Тогда получим, что

<о (х)

Ck ( х ) = (X — хк) со' (хк)

128

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

Ln (х) = 2	f(Xk)	(7)
&=0 (X — Xk ) ( O ' ( x k )
или, более подробно,
£ « (* )= 2 -------------	/"Ы-	(8)
ft=o4 (Xk~~Х!^
i¥*k

2. Интерполяционная формула Ньютона. Эта формула позволя­ ет выразить интерполяционный многочлен Ln(x) через значение /(х) в одном из узлов и через разделенные разности функции f(x), построенные по узлам х0, х„ . .., хп. Она является разностным ана­ логом формулы Тейлора

f(x) = f (х0) + ( х - х 0)Г (Х0) +

Г (Х„) + . ..

Сначала приведем необходимые сведения о разделенных раз­ ностях. Пусть в узлах хк^[а, b), 6 = 0, 1, ..., п, известны значения функции f(x). Предположим, что среди точек хй, 6 = 0 , 1, ..., п, нет совпадающих. Разделенными разностями первого порядка назы­ ваются отношения

	f(Xi, х,) =	f (*/) — f (*,)	i, / = 0	, 1	, . . .,	n, 1ф'].
		Xj — x;

Будем рассматривать разделенные разности, составленные по соседним узлам, т. е. выражения f(x0, х,), /(х,, х2), ..., /(хп_ь хп).

По этим разделенным разностям первого порядка можно по­ строить разделенные разности второго порядка:

f(xlt xd — f(xn, X,)

! (*о> хи х2) =

Х о — Х 0

f (*i, хг, х3) = f ( Х о , х3) — f(xx, Х п )

Х 3— X j

f {хп_х, хп) — f (*„_2, xn_J

/ (Xfi—2, Xfi-1»хп) =

Аналогично определяются разделенные разности более высоко­ го порядка. Например, если известны разности 6-го порядка

f	X j + 1, . . . , Xj+fc) ,	/ (Xj+i,
то разделенная разность (6 + 1)-го порядка определяется как
/ (X/, Х/+1, . . . ,	Х/+А, Х/+£+1) =
	__ /	f {Xj, х)+1......... Х/+&)
		xi+k^ — Xj

5 А . А . С а м а р с к и й , А . В . Г у л и н	129

При	вычислении	разделенных разностей принято записывать
их в виде таблицы
Хо	f (Хо)
	f (*0.		*l)
X 1	f (Xl)		f (xa, xlt x2)
	f	(Xl,	X2)
	f (x%)		/ (*„, x u . . . Xn)

f (xn-2’ xn-1 >xn)

•f ( xn-V xn)

хп f(Xn)

Разделенная разность k-ro порядка следующим образом выра­ жается через значения функции f(x) в узлах:

1+к	/ (*,■)	(0)
/ (xiI ■*■/+.ii • • ■I xi+k) = 2 7^	•

l~' П (Xi — xfr

1фi

L=j

Эту формулу можно доказать методом индукции. Нам потребуется частный случай формулы (9):

/(*„.*!.	*n) = ^		----- =
		П	(JC£— JCf)
		/I=0,
				f (Xi)		( 10)
		(Xi — *o) (*i — Xi)		(X; —	•
					x {_ j ) (XC— * l4 l) • ■• ( x c — Xn )
Интерполяционным многочленом					Ньютона называется много­
член
Рп (X) =		(* — хо)/ (*0. xi) +		{X — XQ) {X — *i) f (XQ, xu x2) + . . .
=	/ (*о) +
	...	+ (X — x0) (x — Xj) ...			(x — x„_!) / (*0, Xi, .. ., xn).	(II)

Покажем, что многочлен Рп(х) совпадает с многочленом Ла­ гранжа (8). Рассмотрим наряду с Ln(x) многочлены L„(x) =} (ха),

130