book1989

Л е м м а 1. Пусть для квадратурной суммы Ih справедливо раз­ ложение (25) и сетка измельчается по правилу hk=qhh, k = 0, 1, ...

. . . , m. Тогда для величин Ihl_v определенных согласно (29), спра­ ведливы разложения

iliL=/ +*№+&Ю1+... +ь!№+о {hlT),

/ = 1. 2, ... , m, k = 1, 2, . . . , m — / + 1,

(30)

где коэффициенты ур не зависят от сетки.

Д о к а з а т е л ь с т в о . Проведем его индукцией по /. При /= 1

равенство (30) выполняется с bP = a{ согласно (25). Предположим, что равенство (30) выполняется при / = / и докажем, что оно вы­ полняется при / = /+1.

Имеем

m

t i l l	=	! J r	S		+ 0 (/г“-		Ь1)-
			i = l				(31)
			rn
/!£ =		/ +	^ ь К		+ оць™ )
			i = l
и, следовательно,		m
йЧ -	= s w (h					° ( C r 1).
		i = l
Далее, подставляя полученную разность в							(29) при }= 1 и учи­
тывая (31), получим
пг				тп
= / + 2		+ — Цу- s ^					+ :° ( С п .
i=i			1~■<?	;=/
т. е.				а/	а/
			(;)ьа£
/ « ^ = 7 +		2		ч	—ч	+	0(/гДГ1)-
				1— 4
	i = l

Отсюда получаем

i=/+i

т. е. равенство (30) выполняется с / = /+ 1, причем

ь г ь р ,

t = 7 + 1,

m.

1 — 4

Лемма 1 доказана.

171

Из леммы следует, что суммы совпадают с интегралом I

kOh '

с точностью до величин О (И. 1),т. е. порядок точности повышается по сравнению с исходной формулой в aja.! раз.

Изложенный метод повышения точности называется методом экстраполяции Ричардсона. Его можно применять не только к квадратурным формулам, но и к самым различным сеточным функ­ циям, если только для них существуют асимптотические разложе­ ния по степеням h. Подробное изложение метода экстраполяции по Ричардсону содержится в книге [23]. Применительно к формуле трапеций данный метод называется методом Ромберга. Существу­ ют стандартные программы вычисления определенных интегралов методом Ромберга. Пример, приведенный в конце п. 4, является частным случаем метода (29), когда /„ — квадратурная сумма, со­ ответствующая методу трапеций, m= 1, <7= 0,5.

Отметим еще, что для формулы трапеций

	I ж	h =	N	0,5 ifс + ft-1) ft,	fi = / ( X i ) ,
			2
			i =l
разложение (25) имеет вид
ь
h = $ f ( x ) d x	+ ^	if' (b) -		Г (a)) - ™ (Г (b) - Г (a)) +
a
+ - ^ - ( f <>	) - / > »		+	. . . + ^ /1 !'</м	( б ) - Л - и М ) + 0 ( 'О .
30240					(32)

Здесь коэффициенты c2r совпадают с коэффициентами разложения функции

	G(h)	h	eh+ 1
		2	eh— 1

в ряд Тейлора:

G (h) = 1 + c2ftz + cji4-]-... + c2r/i2r+ ...

Доказательство формулы (32), называемой формулой Эйлера, можно найти, например, в [2, с. 165].

§2. Квадратурные формулы интерполяционного типа

1.Вывод формул. Будем рассматривать формулы приближенно­ го вычисления интегралов

$Р (x)f(x)dx,	(1)
а

где р(х)>0 — заданная интегрируемая функция (так называемая весовая функция) и f(x) — достаточно гладкая функция. Рассма-

172

триваемые далее формулы имеют вид

	ь	П
	$ р (x)f(x)dxzz 2		Ckf(Xk),	(2)
где	[а, Ь] и ск— числа, k = 0,	1, ...,	п.
	В отличие от предыдущего	параграфа, не будем разбивать от­

резок [а, Ь] на частичные отрезки, а получим квадратурные форму­ лы путем замены f(x) интерполяционным многочленом сразу на всем отрезке [а, Ь]. Полученные таким образом формулы называ­ ются квадратурными формулами интерполяционного типа. Как пра­ вило, точность этих формул возрастает с увеличением числа узлов интерполирования. Рассмотренные в § 1 формулы прямоугольни­ ков, трапеций и Симпсона являются частными случаями квадра­ турных формул интерполяционного типа, когда п = 0, 1, 2, р(х)==1.

Получим выражения для коэффициентов квадратурных формул интерполяционного типа. Пусть на отрезке [а, Ь] заданы узлы ин­ терполирования хк, k = 0, 1 Предполагается, что среди этих узлов нет совпадающих, в остальном они могут быть расположены

как угодно на [а, Ь].		интерполяционным мно­
Заменяя в интеграле (1) функцию f(x)
гочленом Лагранжа
U (х) = 2	<0 (х)	f(Xk),	(3)

(X — X k) (O' (Xk)

6 = 0

где

П

“ (*) = Г[ (х — xfi, = Ц (хк — х,), i=o

получим приближенную формулу (2), где

	_	Г__PJ	■dx, k = 0, 1, . . п.	(4)
Ск		р (х) (О (*)
		J (х- Xk) (O' (Хк)

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

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

1

J / (*) dx	(5)
о

и выберем в качестве узлов точки х о = 0 , дц=0,5, лг2 = 1 .

Квадратурная формула интерполяционного типа, построенная по заданным узлам, совпадает с формулой Симпсона

1

j / (х) dx ж J (/ (0) + 4f (0,5) + /(1)).

(6 )

О

173

Заменим теперь в (5) функцию }(х) многочленом ф(х) наилучшего средне­ квадратичного приближения первой степени. Согласно (18) из § 5 гл. 3 этот многочлен имеет вид

Ф (*) = (/ (1) - f (0)) (* — 0.5) + — (/ (0) + / (0.5) + /(1)).

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

1

j / W dx я J (/ (0) + / (0,5) + / (1)).

(7)

о

не совпадающей с (6 ).

2. Оценка погрешности. Получим выражение для погрешности квадратурной формулы интерполяционного типа. Представим функ­ цию f(x) в виде

f(x)= L n(x) +rn+i(x),

где L„(x) — интерполяционный многочлен для f(x), построенный по узлам х0, хи .. . , х„ и гп+,(х) — погрешность интерполирования. Тог­ да получим

ь	ь	ь
^ Р М / М dx =	^ р (х) Ln (х) dx +	j р (х) гПт1 (х) dx =
а	а	а
		п	Ь

	= 2	(*k) + j р w ' n+i
	k=o	а
(4)	Таким образом, погрешность ф„+1 квадратурной формулы ( 2) ,
	равна
	ь	(8)
	Ф« + 1 — ^ Р (%) Гп+ 1 (•£) dx,
	а
где гп+1(х) — погрешность интерполирования.
	Вспоминая выражение для погрешности интерполирования (см.
(3)	из § 2 гл. 3)

Сп+1 (*) — f(n+l) (£(*)) (я+ О'

получаем

I

Фп+1 --

(Я+1)1

J р {х) ш (х)

(£(*)) dx.

(9)

Отсюда приходим к следующей оценке погрешности квадратур­ ной формулы интерполяционного типа:

Мг

(л+ 1)1 J

а

где Мп+1— шах | / (п+1) (х) |. Из формулы (10) видно, что справед-

хе[а,Ь]

либо следующее утверждение.

174

Квадратурная формула интерполяционного типа, построенная по п + 1 узлу х„, х1г . . . , хп, является точной для любого многочлена степени п, т. е. если f(x) — многочлен степени п и ск— коэффициен­ ты, вычисленные согласно (4), то имеет место точное равенство

J РМ / (*) dx = J] ckf (л*).		(11)
a	k=o
Справедливо и обратное утверждение.
Т е о р е м а 1. Если квадратурная формула

	j	Р (х) / (х)			dkf (**)	О 2)
	a			k=o
точна для любого многочлена степени п, то она является квадра­
турной формулой интерполяционного типа.
Д о к а з а т е л ь с т в о .			Достаточно		показать, что dh = ch, где с*
определены согласно		(4),	k — 0,	1, . . . , п. Рассмотрим		многочлены
	Фа (х)	<о(*)		k = 0, 1, ... , п,
		— x k) о ' ( xk)
	( *
имеющие степень п, и вычислим интегралы
			ь
		h	= j р (х) Фа(х ) dx.
			а
По условию теоремы справедливы точные равенства
			/ а — 2	d.t(pk(xt).
Поскольку			1 = 0
				0,	k ^ t ,
		Фа (XI) =
				1,	k = l,
получаем Ih = dk, k = 0, 1, ..., n.
С другой стороны, согласно (4) имеем
		ь		<Д(-к)
	I* =				Ck.
			(х — xk) ш '
					( xk)
Таким образом, dk=ch, k = 0,				1........ п, что и требовалось.
3.	Симметричные формулы. Для некоторых квадратурных фор­
мул оценка погрешности (10) является грубой, так как она не учи­
тывает симметрии формул.
Рассм отрим , например, ф орм улу Симпсона
	1
	J / (х) dx « -i (/ ( - 1) + 4/ (0) + / (1))					(13)
						173

для функции f( x ) = x \ В данном случае имеем п = 2, /(п+1) (х) =6, поэтому согласно (9) погрешность можно представить в виде

1

Фз = $ “ (*) dx> -I

где (ц(л) =х(х2—1). Благодаря симметричному расположению узлов имеем ф3= 0. В то же время правая часть неравенства (10), равная

1

| \u>{x)\dx = —, “1

отлична от нуля. Таким образом, оценка (10) не является точной для формулы Симпсона.

Квадратурная формула (2) называется симметричной, если

1)п четно;

2)узлы расположены симметрично относительно середины от­

резка	[а, Ь], т. е.
	а4- b	а + Ь	,	А ,	(14)
	---------Xk = X n-k -------— ,		« =	0, 1, . . ., «/2;
3)	Ск = сп. к,	k = 0,	1, . . . ,	п/2.	(15)

Свойство (15) коэффициентов квадратурной формулы опреде­ ляется не только симметричным расположением узлов, но и сим­ метрией весовой функции р(х). Говорят, что р(х) — четная функция

относительно середины отрезка [а,	Ь], если
для всех x e[0 , (b—а) 12].			»,б >
		и п четно, то соот­
Л е м м а 1. Если ск определены согласно (4)
ношения (15) являются следствием условий (14),		(16).
Д о к а з а т е л ь с т в о . Покажем сначала, что
“ (*) = П	(х — xk)		(17)
k=o
является нечетной функцией относительно точки
хп/г=(а+Ь)12.
Имеем
п/2-1
10 (х) = (х хп/2) Г]	(х— Хк) {X — Xn_k),
к=о
откуда, учитывая условия (14), получим
П/2—1
(0(х) = (Х — Хп_2) Т] [(х— хп/2)2—(хк—хп/2)'2].			(18)
k=Q
176

Следовательно, при любом t

rt/2-l

“ (дгп/а + ^) = < Я ^а — (xk—*п/г)2] = — ы(хл/2 — 0.

	Л—О
т. е. функция м(л)	нечетна относительно точки х п/г.
Из формулы	(18) следует также, что <о'(аг) — четная функция

точки Хп/2, поэтому

W'(JC„_A) = CO'(Xa).

Рассмотрим теперь разность

ь

ck— cn_k = \ р М ш(т) ц (1)Ф,

а

где

относительно

(19)

( 20)

		1		1
		Р(х) = (х — xk)m' (хк)	(х -	хп_к) со' (xn_k)
Учитывая (19)		видим, что
, ,	1	х п - x n-k	1	x k ~ x n-k
Р (X) =	“"7 '( * * ) ( X — x k ) (Х — Хп - к )		( Xk )	Н Х — Х п / 2У2— (Х п / 2 ~	ХьУ2}
откуда следует четность ц(дг) относительно точки д:п/2.
Таким	образом, нодынтергальная функция в (20) нечетна относительно сере­
дины отрезка [а, b\, и, следовательно, интеграл				равен нулю. Лемма 1	доказана.
Покажем теперь, что наличие симметрии повышает точность
квадратурных формул. Справедлива
Т е о р е м а		2. Пусть р(х) — четная функция относительно точки

(а+Ь)12 и пусть выполнены условия (14), где п — четное число. Тогда, если квадратурная формула интерполяционного типа (2) точна для любого многочлена степени п, то она точна и для любого многочлена степени п+ 1.

Д о к а з а т е л ь с т в о . Достаточно показать, что формула точи»

для многочлена
f (х) = (х—хпПу	v n / 2 ‘ =0,5 (а+6).
Поскольку
^ р (х) (х — xn/2)n+l d x = 0

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

	In =	2 ckf Ы =	0.
		А—0
Представим /„ в виде суммы i f + In \			где
	л/2-1
l f	= 2	ck(xk- x n/2r \
	k=Q
I f	=	2
	fe=n/2+l

177

Из условий (14) получим

		Д	2 , =	2	Ck ( x n„ -	x n- k)n+1
				к=п/2Ы
	1(п1=	2 Cn~i (Хп/* — х^м = — 2				Cn-k (Хк~~Хп^		1-
Последнее равенство справедливо в силу того, что п — четное число.
Таким образом, получаем
				Л / 2 - 1
		in= / ? ’ +	i	f = 2	fa -	fa -	*«//*■
				fc=o
Согласно лемме 1 имеем сЛ=с„_А, &=0,						1,	и/2—1,	т. е. /„= 0 ,
что и завершает доказательство теоремы 2.
4.	Формулы Ньютона— Котеса. Численная устойчивость ква­
дратурных		формул.	Формулами Ньютона — Котеса					называются

квадратурные формулы интерполяционного типа, построенные на равномерной сетке, когда хк—xk-i = h, k = 1, 2, . . . , п.

Различают два типа формул Ньютона — Котеса: формулы замк­ нутого типа и формулы открытого типа. В формулах замкнутого типа х0= а и хп = Ь, а в формулах открытого типа хотя бы один из узлов х0 или хп не совпадает с соответствующей граничной точкой отрезка [а, Ь]. Для простоты изложения рассмотрим лишь случай формул замкнутого типа, когда xh = a+ kh, k = Q, 1........п, hn = b—а.

В случае равномерной сетки можно упростить выражения для коэффициентов квадратурных формул. В формуле (4) сделаем за­ мену x = a+ th, O^Lt^in. Простые выкладки, которые мы опускаем, показывают, что в результате замены формула (4) примет вид

сц = ф — a) b f\

где
Ы =	f — fp (а + /А)' (l ~	11 •"	~ ">- Д	■	(21)
	k\ (п— k)\ \' п Jо	t — k	'\
Отметим,	что формулы Ньютона — Котеса с		редко		ис­

пользуются из-за их численной неустойчивости, приводящей к рез­ кому возрастанию вычислительной погрешности. Причиной такой неустойчивости является то, что коэффициенты формул Ньютона — Котеса при больших п имеют различные знаки, а именно при п ^ 5; 10, р(х) = 1 существуют как положительные, так и отрицатель­ ные коэффициенты.

Остановимся подробнее на значении знакопостоянства коэффи­ циентов для устойчивости вычислений.

Рассмотрим квадратурную сумму
/„ = 2 <*/(**)•	(22)
k=0

178

Предположим, что значения функции f(x) вычисляются с некото­ рой погрешностью, т. е. вместо точного значения получаем прибли­ женное значение f (xh) = f(xh) + 8 h. Тогда вместо /„ получим сумму

/л '•= ^ Ск(/ (-*■'<) Ф" 8*) —	+ б/„,
к=о
где
би = 2	(23>
й=о
Поскольку квадратурная формула (2),	(4) точна для /(*)== I,
имеем

пЬ

2	Ck = ^ р (х) dx,
k=Q	а
т. е. при р(х)>0 сумма
	= М	(24)
	к=0

ограничена числом Л4>0, не зависящим от п.

Предположим теперь, что все коэффициенты ch неотрицательны.

Тогда из (23), (24)	получим оценку
| 8 / „ | < 2		|с*| | б * [ = 2	с*|б*|^(гпах (6й1)М,	(25>
^		,	о<ck<n
k=0		k~0

которая означает, что при больших п погрешность в вычислении квадратурной суммы (22) имеет тот же порядок, что и погрешность в вычислении функции. В этом случае говорят, что сумма (22) вы­ числяется устойчиво.

Если коэффициенты ск имеют различные знаки, то может ока­ заться, что сумма

2 1 * 1

к=о

не является равномерно ограниченной по п и, следовательно, по­ грешность в вычислении /„ неограниченно возрастает с ростом п. В этом случае вычисления по формуле (22) будут неустойчивы и пользоваться такой формулой при больших п нельзя.

Таким образом, если необходимо сосчитать интеграл (1) более точно, то имеются две возможности. Во-первых, можно разбить весь отрезок \а, Ь] на несколько частичных отрезков и на каждом из частичных Отрезковприменить формулу Ньютона -— Котеса с не­ большим числом узлов. Полученные таким образом формулы назы­ ваются составными квадратурными формулами. Они часто приме­ няются на практике, хотя и не являются достаточно экономичными, поскольку требуют многократного вычисления значений функции

179

,f(x). Во-вторых, можно попытаться выбрать узлы квадратурной формулы так, чтобы полученная формула имела большую точность, чем формула Ньютона — Котеса с тем же числом узлов. В следую­ щем параграфе рассматривается один из методов, основанных на

.выборе узлов квадратурной формулы, а именно метод Гаусса. Он приводит к квадратурным формулам с положительными коэффи­ циентами при любых п и существенно более точным, нежели форму­ лы Ньютона — Котеса.

§3. Метод Гаусса вычисления определенных интегралов

Т.Постановка задачи. В предыдущем параграфе предполага­ лось, что узлы квадратурных формул заданы заранее. Было пока­ зано, что если использовать п узлов интерполяции, то получим ква­ дратурные формулы, точные для алгебраических многочленов сте­ пени п—1. Оказывается, что за счет выбора узлов можно получить ■квадратурные формулы, которые будут точными и для многочленов

степени выше п—1. Рассмотрим следующую задачу: построить ква­ дратурную формулу

Ьп

$Р (x )f(x )d x ^ '2 } ckf(xk),		(I)
a	k=i

которая при заданном п была бы точна для алгебраического много­ члена возможно большей степени. Обратим внимание, что здесь в отличие от § 2 для удобства изложения нумерация узлов начинает­ ся с k= 1.

В настоящем параграфе будет показано, что такие квадратур­ ные формулы существуют. Они называются квадратурными форму­ лами наивысшей алгебраической степени точности или формулами Гаусса. Эти формулы точны для любого алгебраического многочле­ на степени 2п—1.

Итак, потребуем, чтобы квадратурная формула (1) была точна для любого алгебраического многочлена степени ш. Это эквива­ лентно требованию, чтобы формула была точна для функций f(x) = = ха, а = 0, ш. Отсюда получаем условия

пЬ

^	CkX%= \.р{х) х?йх, сс = 0, 1, . .., т,	(2)
k= i	а

которые представляют собой нелинейную систему т+ I уравнений -относительно 2п неизвестных

Ci, С2, . . . , С п ,

Х 1 , Х 2, . . . , Х п .

Для того чтобы число уравнений равнялось числу неизвестных,

.надо потребовать т = 2п—1. В дальнейшем будет показано, что си­ стема (2) при т = 2п—I имеет единственное решение. Однако сна­ чала рассмотрим несколько частных случаев, когда решение систе­ мы (2) можно найти непосредственно.

180