book1989

Тогда система (13) решается в явном виде,
Ck=(f,(fk)H,	k = 0 , l, . . . ,n ,	( 16)
а погрешность приближения определяется формулой
II/ — ф||« =	Ш н - 2 &	(17)
	6=0

Числа ск, определенные согласно (16), называются коэффици­

ентами	Фурье элемента f^ H	по ортонормированной системе
( Ф а } а = о ,	а обобщенный многочлен
		П
	ф =	2 ^
		k--=a

называется многочленом Фурье.

Г Л А В А 4

ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ И ДИФФЕРЕНЦИРОВАНИЕ

§1. Примеры формул численного интегрирования

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

ь	(1)
I = <\j f(x)dx,

а

основанные на замене интеграла конечной суммой

/« = 2 с*/(**),	(2)
6=0

где ск— числовые коэффициенты и хк— точки отрезка [a, b], k =

=0, 1, . . . , п. Приближенное равенство

Ьп

^ / (х) dx ж 2 6V (**)

a k —Q

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

Ьп

Чп = 5 / w dx — 2 (xk)

а6 = о

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

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

161

Рис. 5. Геометрический смысл формулы прямоугольников

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

Введем на [а, Ь] равномерную сетку с шагом h, т. е. множество

точек

аи = {х2=а + Иг, i= 0, 1........ JV, hN = b—а},

и представим интеграл (1) в виде суммы интегралов по частичным отрезкам:

Ь	N	Х 1
=	2	$ f(x)dx.	(3)
a	i=i JC

Для построения формулы численного интегрирования на всем отрезке [а, Ь] достаточно построить квадратурную формулу для

интеграла

Х1
$ / (х) dx				(4)
Х С- 1
на частичном отрезке [х{- ±, х,] и воспользоваться свойством				(3).
2.	Формула прямоугольников.
Заменим интеграл (4)			выраже­
нием f(x ^ 42)h, где х^ч2= х2—0,5 h.
Геометрически		такая		замена
означает, что площадь			криволи­
нейной трапеции ABCD заменяет­
ся	площадью	прямоугольника
ABC'D' (см. рис.		5). Тогда полу­
чим формулу
	5 / (х) dx ~	/ (лщД /г,		(5 )

, .. . „

которая называется формулой прямоугольников на частичном

отрезке [х,.,, х;]. Погрешность метода (5) определяется величиной

Фг= J / (х) dx — f (х;_%) h,

которую легко оценить с помощью формулы Тейлора. Действитель­

но, запишем	в виде	хс
	Ф < =	J	(I (-0 — / (*<-'/,))dx	(6)
	xi-1
и воспользуемся разложением
/ М =	/ (Xi-V,) + ( х	-	г ( * - • /,) +	Г & ) ,
162

где t,i= t,i(x) е [ х ;_ь Xi], Тогда из (6) получим
	♦< =	]	^ ^ - П	М Л .
		*£-1
Обозначая М2,г =	шах		|/" (х )|,	оценим	"ф,- следующим обра-
зом:
Ы	^ м	а>£	j	= f	М2„
		*£-i

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

(7)

т. е. формула имеет погрешность О (Я3) при Я->-0.

Заметим, что оценка (7) является неулучшаемой, т. е. сущест­ вует функция f(x), для которой (7) выполняется со знаком равен­ ства. Действительно, для f(x) = (x—xt- >к)2 имеем М2|<=2, /(х,-./,) =

= 0 и
xi
j / (Х) dx —f (х^А) Я =	M2J.

*£— 1

Суммируя равенства (5) по i от 1 до N, получим составную фор­

мулу прямоугольников
^ ( х ) й х ^ ^ П х ^ А)1г.		(8)
a	i = i
Погрешность этой формулы

¥ == j f (х) dx — ^ f (х<-'/,)h

а*=1

равна сумме погрешностей по всем частичным отрезкам,

N	N	Xir.	tr— г \2	r&)dx.
^ 2 ^ = 2		Г-
£=l	1=1
Отсюда, обозначаяМ2 — шах		|f"(x) |, получим
хе[а.Ь)
|¥ 1 С	м2т а		ft2 (b— а)	М2i	(9)

2424

т.е. погрешность формулы прямоугольников на всем отрезке есть величина О (Я2).

6* 163

В этом случае говорят, что квадратурная формула имеет вто­ рой порядок точности.

З а м е ч а н и е . Возможны формулы прямоугольников и при ином распо­ ложении узлов, например такие формулы:

Ь	N	Ь	N
J /w dx~ 2 hf (xi-1)•		5 /м ^~ 2 А/w•
a	i= l	а	i=i

Однако из-за нарушения симметрии погрешность таких формул является вели чиной 0(h).

3. Формула трапеций. На частичном отрезке эта формула име­ ет вид

j f(x)dx« /K--i)+ /(*,)

(10)

2

*1-1

и получается путем замены подынтегральной функции f(x) интер­ поляционным многочленом первой степени, построенным но узлам х{- 1, хи т. е. функцией

	U.i (х) =			- ( ( х — Х : ^ ) f ( X i )		— (х — Xi) / ( X f - 0 ) .
				h
Для оценки погрешности достаточно вспомнить (см. п. 1 § 2
гл. 3),	что			(х — X; А (х — хЛ
		/ (х) -					/" & (х)).
				Lui (х) = -------	^	--------
Отсюда получим
,	f n w		/ K'-i) + / (*f)		,
ф; =	\ / (х) d x -----------			;--------	h —
	=	j	(f(x) — Lu i (x))dx =			(X--	X; ,) (X--	X,)
						~	---------	- Г Ы * ) ) * х
и, следовательно,						XL-
					12			(П)
Оценка (11)			неулучшаема, так как в ней достигается равенст­
во, например, для f(x) = (х—х{)2.
Составная формула трапеций имеет вид
\ f ( x ) d x z z ^		-------		- ------- h =
а	‘=1			= h ( 0 ,5 f0+		/ 1 + • ■ •	+ h -г +	0 ,5 fN), (1 2 )
где f i =	f { X i ) , i	= 0 ,	1	hN = b—a.

164

Погрешность этой формулы оценивается следующим образом:

1Я ^	\2	м 2,	= шах | Г (х) |.
			х^[а,Ь]

Таким образом, формула трапеций имеет, так же как и форму­ ла прямоугольников, второй порядок точности, Чs = 0{Нг), но ее по­ грешность оценивается величиной в два раза большей (см. (9)).

4.Формула Симпсона. При аппроксимации интеграла (4) заме­

ним функцию f(x)			параболой, проходящей через точки (Xj, f(x})),
j = i—1, i—0,5, i, т. e. представим приближенно f(x) в виде
		f(x) ~ L 2i(x),			х е [х (_ь Xi],
где	L2i(x )— интерполяционный			многочлен		Лагранжа второй	сте­
пени,
1 2Л	( х ) =	- 2 — {(х Х;_у2) (х — X i) /f_ t —
		/г
		— 2 (х — х.ч) (х —х,)			+ (х — хг) (х — x/_Vs) /;}.		(13)
Проводя интегрирование, получим
		xi
		j ‘ L2.i(x)dx = JL{fi-l + 4fi-Vl + fi),				h = xi — Xi-1.
		xi-1
	Таким образом, приходим к приближенному равенству
		J■ч	/ ( X ) r f x ^ A	(/[._1 +4/i.,+/i))v			(14)
		xi-1
которое		называется формулой Симпсона или формулой парабол.
	На всем отрезке [а, Ь] формула Симпсона имеет вид

ьN

Г/ (	х	)	1+ 4/м* +		A) =
а		>'=1
=	~D[/о +	/л/ +	2 (/х +	/2 -+- . . .	+ / JV—i) + 4 (/;/ +	f‘k +	. . . + fJv-vi)].
Чтобы не использовать дробных индексов, можно обозначить
	X i = a + 0,5W,			fi= f(Xi),	i = 0, 1, ... , 2N,	hN = b—a
и записать формулу Симпсона в виде
	ь
	j* / (х) dx ^			[/о + f iN	+ 2 (/2 + / 4 + . . .	+ / 2Л/-2) +
					+ 4 (/j + /3 +	...	-f fzx-i)]. (15)

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

165

степени, т. е. имеет место точное равенство

*i h

j / (A) d x = - ^ (/ t_x + 4 f c - Vt + /у),

Xi~i

если f(x)= a0+ alx+a2x2-\-a3x3. Это утверждение нетрудно прове­ рить непосредственно, что и предоставляется сделать читателю.

Для оценки погрешности формулы Симпсона воспользуемся интерполяционным многочленом Эрмита. Построим многочлен третьей степени Я3(а) такой, что

Н3(A.-I) =	/ (Xi-j),	И, (xc-yj = f (ау-.Д,
я ; (АУ-Д =	/ ' (АУ_%),	Я 3 (а£) = / (А[).

Из § 3 гл. 3 известно, что такой многочлен существует и единствен. Он построен в явном виде в примере из п. 2 § 3 гл. 3. Однако нам даже не потребуется явный вид многочлена Я3(а). Вспоминая, что формула Симпсона точна для любого многочлена третьей степени, получим

Xi

h

(А-;_у2) +

Яз

(ху)) =

J

Нц(х)йх — — (Я3(х^) + 4Яз

=

4

- ( f t_1

+

4 /l-_y1+

/ i).

(16)

Представим теперь /(х) в виде

О

f{x) = Я3(х) + /у(х), хеЕ[х<-и хЛ,

(17)

где

г,(а ) — погрешность

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

многочленом

Эрмита

Я3(а). Интегрируя (17) и учитывая (16), получим

*t

f (A) dx — -J- (/,•_! +

Xi

ri (A) dx.

j*

4/,_I/2-f /,) =

|

(18)

*i-i

x;-i

Согласно (14) из § 3 гл. 3 имеем

ПМ = — 7Гл~~ (х —

— Xl'-%)2 (х —

24

поэтому из (18) для погрешности фу

формулы

(14)

получаем

оценку

XI

Af

•

(*

A('_I/S)2 (а — X i - i ) d x ,

| ф £| < —

(а — X i) (а —

*i-l

где

Miti=

sup

| / IV (A)|.

Вычисляя интеграл, приходим окончательно к оценке

(19)

2880

166

Погрешность составной формулы Симпсона (15) оценивается

так:

~ а ) /И*’ hN = b — a, М4=	sup |/ ‘V )|.
zooU	x^[a,b]

Отсюда видно, что формула Симпсона существенно точнее, чем формулы прямоугольников и трапеций. На частичном отрезке она имеет точность О(ft5), а на всем отрезке — О(ft4).

Приведем вывод формулы Симпсона, основанный на методе экстраполяции. Метод экстраполяции состоит в следующем. Проведем два расчета по формуле трапеций ( 1 2 ), первый расчет с шагом ft, когда вычисляется сумма

N	fi+fi-i ft.			/(■= / (х[) i		*i =	a-f ift.
1=1			h=(b—a)/N,
и второй расчет — с шагом 0,5ft, когда вычисляется сумма
N	\	(	fi-</2+	fi	h
/ 1—! + / i-%
= 2					~2 ~
	=	2	№ -i +	2h	%+ fi)	\ .	fi-Vl = / (*i - 0,5ft).
		i=i

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

rh= I + Clh4 0 (h*),

где 1 — исходный интеграл ( 1 ) и Cj — постоянная, не зависящая от ft. Точно

так же

/ Л/2= / + С1 ^ )

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

т. е. выражение

1

3 hpi — 3 h

совпадает с интегралом 1 с точностью до величин О (А4).

В данном примере не обязательно проводить расчет на двух сетках, так как можно построить явное выражение для суммы /л. Действительно,

N	N	fi+ fi-i
4 / Л/2 >k= 2	(/i-i+ 2A--y, + / t) f t - 2		ft =
i=i	i=i	2
			N
				Д-i + 4/r~%+ fi ht
			- 2
				2
			i = i
	jh= 2 /‘-1 + 4/t'-%+ A'		ft.

Таким образом, снова получим квадратурную формулу Симпсона.

167

5. Апостериорная оценка погрешности методом Рунге. Автома­ тический выбор шага интегрирования. Величина погрешности чис­ ленного интегрирования зависит как от шага сетки h, так и от глад­ кости подынтегральной функции f(x). Например, в оценку (11), наряду с /г, входит величина

M2ti=	шах	| /"(*)!,
	xG[x	хi]

которая может сильно меняться от точки к точке и, вообще говоря, заранее неизвестна. Если величина погрешности велика, то ее мож­ но уменьшить путем измельчения сетки на данном отрезке [х{- и лД. Для этого прежде всего надо уметь апостериорно, т. е. после проведения расчета, оценивать погрешность.

Апостериорную оценку погрешности можно осуществить мето­ дом Рунге, который мы поясним сначала на примере формулы тра­ пеций. Пусть отрезок [а, b] разбит на частичные отрезки [х.-ь лД,

i= 1, 2, . .. ,	N,	х„ = а, xN = b, имеющие,			может быть, разную длину
hi = Xi—Xt-i.	На		каждом	частичном отрезке		применяется формула
трапеций			xi
			Л = J / (х)dx «		in =	Ih,i.
			xi-1
Согласно (11)		имеем
				h - I h,i~Cihl		(20)
где константа		с,	зависит	от гладкости	f(x)	и заранее неизвестна.

Измельчим на отрезке [л:,_ь лД сетку в два раза и повторим расчет с шагом 0,5h, т. е. вычислим сумму

Аг/2.1 — ifС-1 + 2 ft-'/,	+ fi) ~ ~ .
Тогда согласно (20) будем иметь	4
А-//,/2 . 1	(21)

Из соотношений (20), (21) можно исключить константу ct и по­ лучить оценку погрешности, которая содержит лишь известные ве­ личины /м, Ih/2 i:

	Ih.i)i
IГ	lh/i,i — — (li ' Ih.i) ^ {Ih/2,i hi, .
	8
Разумеется,	метод Рунге можно применять и для оценки по­

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

т. е. A—	Тогда
	1	т ~ ( hi\m
	‘ С	'/1/2,1—Ci\	I
168

откуда получим

(2 2 )

г	^h , t	(23)

Возможность апостериорно оценивать погрешность позволяет вычислять интеграл (1) с заданной точностью е>0 путем автома­ тического выбора шага интегрирования h{. Пусть используется со­ ставная квадратурная формула

	N
I ~ I h =	^
	i=i

где Ih,i— квадратурная сумма на частичном отрезке, причем на каждом частичном отрезке используется одна и та же квадратур­ ная формула (например, формула трапеций, Симпсона и др.). Про­ ведем на каждом частичном отрезке [х ^ , xt] все вычисления дваж­ ды, один раз — с шагом h{ и второй раз — с шагом 0,5/i; и оценим погрешность по правилу Рунге (23).

Если для заданного е>0 будут выполняться неравенства

I U - h	/2,1	Ih/2.i	Ih,i I	eh,	N,	(24)
		2m—	1	i = 1, 2,
то получим
				N
				1 = 1
т. e, будет достигнута заданная точность е.					(24)	не бу­
Если же на каком-то из частичных отрезков оценка

дет выполняться, то шаг на этом отрезке надо измельчить еще в два раза и снова оценить погрешность. Измельчение сетки на данном отрезке следует проводить до тех пор, пока не будет достигнута оценка вида (24). Заметим, что для некоторых функций f(x) такое измельчение может продолжаться слишком долго. Поэтому в соот­ ветствующей программе следует предусмотреть ограничение свер­ ху на число измельчений, а также возможность увеличения е.

Таким образом, автоматический выбор шага интегрирования приводит к тому, что интегрирование ведется с крупным шагом на участках плавного изменения функции f(x) и с мелким шагом — на участках быстрого изменения f(x). Это позволяет при заданной точности е уменьшить количество вычислений значений f(x) по сравнению с расчетом на сетке с постоянным шагом. Подчеркнем,

что для нахождения сумм / h/2,( не надо пересчитывать значения
f(x)	во всех узлах, достаточно вычислять f(x) только в новых узлах.
6.	Экстраполяция Ричардсона. Способ повышения точности

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

169

Предположим, что для вычисления интеграла (1) отрезок [а, 6] разбит на Л/' равных отрезков длины /г= (Ь—и)/Л/ и на каждом частичном отрезке применяется одна и та же квадратурная фор­ мула. Тогда исходный интеграл / заменяется некоторой квадратур­ ной суммой Д, причем возникающая погрешность зависит от шага

сетки h. Для некоторых квадратурных формул		удается получить
разложение погрешности Д—I по степеням h.		Предположим, что
для данной квадратурной суммы Д существует разложение
Д = I +	+ а Х г + . . . + a S ' n +	о (ham+\	(25)
где 0 < a i< a 2< . . . < а т< а т+, и коэффициенты		а,- не зависят	от h.

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

Вычислим приближенно значение интеграла I по данной квад­ ратурной формуле на последовательности сеток с шагами h 0 = h, hlt hi, ... t hm. Для определенности будем предполагать, что сетка

измельчается по геометрической		прогрессии,	т. е. h k = qkh 0,	k = 0,
1, . . . , m, где <7^(0,	1). Вычисляя квадратурную сумму Д при раз­
личных значениях h,	получим величины h k, Л = 0, 1, . . . , т , причем
согласно (25) будем иметь
Ihk= I -f- dyhf{' -Д ajift’ -Д		... -Д amhk	+ О (hit +l).	(26)

Обозначим / (0) = /, I = Ihk. Исключая коэффициент а, из соот­

	ношений
	=	/ + aX U + О (htr),
	1%=*1+ * № + (>(№),
	получим					(27)
	/ =	//<1	1 +	0(С -1),
	где обозначено			/(и_	/И)
	И2)	г(1)	I			(28)
				hk	hk- l
	По формулам (28) можно вычислить величины					&= 0, 1, ...

	. . . , m—1. Согласно (27), они дают более точное, чем /j^,					прибли­
	жение к интегралу I. Этот процесс повышения точности можно про­
	должить, вычисляя величины			с помощью рекуррентных соотно­
	шений			/(/> _ /(/)
	/(/+1)	/О')	I
				.	a.- ’	(29)
	Ihk-i	lhk-i	f
				1—q 1
	/ ' = 1 , 2 ,		k = \ , 2 ...........m —j + l ,
	=	V	* =	0.	1. ■■■ , m.
	170