Schisla_1

ности

		b		b − a	n
R (f) =		ρ(x) f(x) dx				C f(x ).
			−		∑
n	a			2		k k
	∫				k=1

Если при выводе квадратурных формул задаются узлы xk (как правило, равноотстоящие), а разыскиваются коэффициенты Ck, то получаются формулы Ньютона - Котеса. К группе этих формул при ρ(x) = 1 относятся хорошо известные формулы прямоугольников, трапеции, Симпсона.

4.1Формулы Ньютона - Котеса

Пусть в заданных n узлах xk [a, b], k = 1, 2, ..., n известны значения подынтегральной функции fk = f(xk). Функцию f(x) с помощью интерполяционного полинома Лагранжа Ln(x) можно записать в виде

	f(x) = Ln(x) + rn(x),							(4.2)
где		n			n
							−
		∑k			∏̸
Ln(x) =		f(xk)				x − xj			,
		=1			j=1;j=k	xk		xj
	f(n)(ξ(x))			n
а rn(x) = f(x) − Ln(x) =				∏k			– погрешность интерполя-
		n!			(x − xk)

=1

ционного многочлена Лагранжа, ξ(x) (a, b). Тогда из (4.2) и (4.1) получается выражение

∫b ∫b ∫b

I(f) = ρ(x) · f(x) dx = ρ(x) · Ln(x) dx + ρ(x) · rn(x) dx,

a a a

первое слагаемое которого определяет квадратурную формулу Sn(f), а вто-

рое ее погрешность Rn(f).
Здесь удобно сделать замену
x = a + b		+ tb − a		,	если	a	6	x	6	b,	то −	1	6	t	6	1.	(4.3)
	2		2

Узлы интерполяции можно записать в виде

41

xk =	a + b	+ dk	b − a	,
			2
2

где dk [−1, 1]. После замены интерполяционный многочлен примет вид

n (

)

n

n

−

a + b

∑

∏̸

L

+ t

b − a

= f

t − dj

,

2

2

k=1

k j=1;j=k dk

dj

а квадратурная формула

b

2

1

(

Sn(f) = ∫a ρ(x) Ln(x) dx =

2

)

∫

2

2

)

n (

2

=

b

− a

ρ

a + b

+ t

b

− a

L

a + b

+ t

b − a

dt =

−1

b − a

n

1

a + b

b − a

n

t − dj

b − a

n

=

f

ρ

+ t

dt =

C f ,

∏̸

∑

2

∑

−1

(

2

2

dk

−

dj

2

k k

k=1

k ∫

)j=1;j=k

k=1

1

a + b + tb − a

n

t − dj dt

C = ρ

(

)j=1;j=k

где k

∫

— коэффициенты Ньютона

2

2

dk

−

dj

1

∏̸

−

-Котеса. Эти коэффициенты зависят от пределов интегрирования a и b,

если ρ(x) ≠ const .

Для погрешности получается выражение

								b
R	f	I f			S	f = ρ(x) r				(x) dx =
	n( ) = b		( ) −			n( )		n ∫a	n
	a			f(n)(ξ)				∏
	= ∫		ρ(x)		n!			k=1(x − xk) dx,

из которого после замены (4.3) следует оценка

∫1

(b − a)n+1 |Rn(f)| 6 n! 2n+1

−1

ρ a + b + tb − a f(n)(ξ)					n	(t d )dt .
					∏
(	2	2		)	k=1	− k		(4.4)

Из формулы (4.4) видно, что, даже если производная f(n)(x) подынте-

гральной функции имеет на [a, b] интегрируемую особенность, то погреш-

42

ность Rn(f) ограничена (отметим, что интерполяционный многочлен сте-

пени n − 1 для такой функции строить нельзя). Если же существует Mn =

max f(n)

(x)

x [a;b]

|

|, то

|

n

| 6

n!

(

2

)

n+1 ,

R

(f)

C

Mn

b − a

1

(

n

1

a + b

b − a

∏

где

C =

ρ

+ t

t

d

dt.

−

2

2

| −

k|

∫

) k=1

Если f(x) является произвольным многочленом степени m 6 n − 1, то

f(n)(x) = 0. Это означает, что квадратурная формула Ньютона– Котеса точна (Rn(f) = 0) для любого многочлена степени m 6 n − 1, если отсут-

ствуют ошибки округления. В этом случае говорят, что формула Ньютона–

Котеса имеет алгебраическую точность n − 1.

Если вес ρ(x) – четная функция относительно середины отрезка (a+b)/2

(ρ (	a + b	+ x) = ρ (	a + b	− x)), тогда узлам, расположеным симмет-
	2		2

рично относительно середины (dk = −dn+1−k), соответствуют равные коэффициенты.

Действительно, если сделать замену t на −t и воспользоваться равен-

ством dk = −dn+1−k, то получится следующая цепочка равенств

1

a + b

b − a

n

t − dj

C = ρ

+ t

dt =

k

1

(

2

2

)j=1;j=k dk

−

dj

1

∫

1

−∫

n

n

∏̸

a + b

∏̸

−

=

− ρ

−

t

b − a

−t + dn+1−j

dt =

−∫

− 1

(

2

2 )j=1;j=k

dn+1−k + dn+1−j

(

∏̸

−

= ρ

a + b

+ t

b − a

t

− dn+1−j

dt = C

.

1

2

2

)j=1;j=k dn+1−k

dn+1−j

n+1−k

Если ρ(x) – четная, а f(x) – нечетная функции относительно середины

отрезка (f (	2	+ x) = −f (	2 − x)), то I(f) = ∫	ρ(x) f(x) dx =
	a + b		a + b	b

a

0, а в силу свойства коэффициентов и Sn(f) = b − a ∑n Ckf(xk) = 0.

2

k=1

43

Такие квадратурные формулы, точны (Rn(f) = 0) для любой нечетной относительно середины отрезка функции.

Произвольный многочлен степени n можно записать в виде

Pn(x) = anxn + an−1xn−1 + ... + a1x + a0 = an (x −	a + b	)	n
			+ Pn−1(x),
	2
где многочлен Pn−1(x) степени n − 1.

Такие "симметричные" квадратурные формулы, построенные по нечетному числу n симметричных узлов, имеют алгебраическую точность n (на единицу большую, чем в других случаях). Уточненная оценка погрешности получается, если использовать интерполяционный полином с двукратным узлом x = (a + b)/2. Эта оценка выражается через (n + 1)-ую производную

и имеет вид

(

n+2

R

(f)

|

6

C

Mn+1

b − a

,

(n + 1)!

где

| n

2

)

M

n+1

= max f(n+1)(x)

x [a;b] |

|,

1

) · |

a + b

b − a

n

C =

ρ

+ t

t

d

t

d

dt.

( 2

−

(n+1)=2| j=1

−

j|

∫

2

|

−

1

∏

Пример. Вывод формулы Симпсона. Пусть ρ(x) = 1. Интерполяционный многочлен Лагранжа, построенный по трем узлам x1 = a, x2 = (a + b)/2 и x3 = b, имеет вид

	L3(x) = f(a)				(x − (a + b)/2)(x − b)	+
					(a − (a + b)/2)(a − b)
	+f (		)
		a + b			(x − a)(x − b)

2((a + b)/2 − a)((a + b)/2 − b)

+f(b)

(x − a)(x − (a + b)/2)

=

2

(f(a) (x

−

(a + b)/2)(x

−

b)

(b − a)2

(b − a)(b − (a + b)/2)

−

−2f((a + b)/2) (x − a)(x − b) + f(b) (x − a)(x − (a + b)/2))

b

Если вычислить ∫a

L3(x) dx, то получится формула Симпсона

I(f)

≈

S

(f) =

b − a

1

f(a) +

4

f

a + b

+

1

f(b) ,

(4.5)

2

3

3

2 (

3

3 (

)

)

44

погрешность которой

R3(f) = I(f) − S3(f) =

b

′′′(3!( ))(x − a)(x − b)(x −

2

) dx

∫a

f

ξ x

a + b

f′′′(x)

max

f′′′(x)

| 6

M

Если

существует и ограничена ( a6x6b

|

3), то

|

R

(f)

6

M

(b − a)4

.

196

3

|

3

При выводе формулы Симпсона использовались расположенные симметрично относительно середины отрезка три узла, значит можно получить уточненную оценку погрешности. Пусть у подынтегральной функции существует f(4)(x) на отрезке [a, b]. Запишем интерполяционную формулу Ньютона с разделенными разностями для узлов x1 = a, x2 = b и двукрат-

ного узла x3 = x4 =	a + b
2
													a + b
Q4(x) = f(a) + f(a, b) (x − a) + f(a, b,														) (x − a)(x − b)+
													2
+f(a, b,		a + b		,		a + b				) (x − a)(x − b)(x −					a + b	) =
		2				2									2
					a + b						a + b						a + b
= L3(x) + f(a, b,									,			) (x − a)(x − b)(x −								).
					2						2					2
Погрешность такого полинома имеет вид
				1				f(4)(ξ(x)) (x − a)(x − b) (x										a + b			)	2
r(x) = f(x) − Q4(x) =																−						.
							4!											2

Квадратурная формула, построенная с помощью многочлена Q4(x), име-

ет вид

		b				b
b	Ib(f) = ∫a		f(x) dx ≈ ∫a				Q4(x) dx =
= ∫a	L3(x) dx + ∫a	f(a, b,	a + b	,	a + b	) (x − a)(x − b)(x −		a + b	) dx
			2		2			2

Первое слагаемое совпадает с формулой Симпсона, а второе равно нулю, так как подынтегральная функция нечетная относительно середины

45

отрезка. Если подынтегральная функция имеет ограниченную четвертую

производную max |f(4)(x)| 6 M4, то для этой формулы получается оценка

a6x6b

погрешности

|R(f)| 6

a

b

f(a, b, a + b, a + b)(x a)(x b)(x

a + b)

dx M (b − a)5 ,

6

∫

2

2

−

−

−

2

6

4

2880

Итак,

используя n произвольных узлов можно вывести

квадратурную

формулу Ньютона - Котеса, имеющую алгебраическую точность (n − 1). Если же значение n нечетное, весовая функция четная относительно се-

редины отрезка интегрирования (x =	a + b	) и узлы выбраны симметрич-
	2

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

Оказывается для некоторых весовых функций (например, для положительных почти всюду на [a, b]) существуют квадратурные формулы с мак-

симальной алгебраической точностью 2n − 1.

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

x		=	a + b	+	a − b	cos	π(2k + 1)		, k = 0, 1, ..., (n	−	1).
	k						2n
		2		2

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

|

f(x)

−

L

(x)

Mn

1−2n

b

−

a n,

M

max f(n)(x),

n! 2

n

| 6

(

)

n = a6x6b

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

Если использовать эти узлы для вычисления

I(f) = ∫

b

f x

√

( )

dx,

1

x2

a

−

46

1
с весовой функцией ρ(x) = √1 − x2	, то получаются квадратурные фор-

мулы с равными коэффициентами — квадратурные формулы Чебышева.

Для отрезка [−1, 1] формула Чебышева (Эрмита) имеет вид

	1					π n
	I(f) = ∫	f x
		√		( )	dx ≈ Sn(f) =		k=1 f(xk).
						n
			1	x2
	−1			−			∑

Формула Чебышева относится к группе формул типа Гаусса, имеющих наивысшую алгебраическую точность 2n − 1.

4.2Формулы типа Гаусса

Будем разыскивать такие узлы xk и такие коэффициенты Ck квадра-

турной формулы

∑n

Sn(f) = Ckf(xk),

k=1

чтобы она имела как можно более высокую алгебраическую точность m. Если выполняются равенства I(xi) = Sn(xi), i = 0, 1, ..., m, то для произвольного многочлена Pm(x) степени m также справедливо равенство I(Pm) =

Sn(Pm). Искомые коэффициенты Ck и узлы xk должны удовлетворять системе уравнений

µ0	≡	∫a b		ρ(x) dx = C1 + C2 + ... + Cn,
µ1	≡ ∫a b ρ(x)x dx = C1x1 + C2x2 + ... + Cnxn,
	≡ ∫		b		(4.6)
				ρ(x)x2 dx = C1x12 + C2x22 + ... + Cnxn2
µ2					,

a

− − − − − − − − − − − − − − − − − − −−

∫b

µm ≡ ρ(x)xm dx = C1xm1 + C2xm2 + ... + Cnxmn .

a

47

Находить 2n неизвестных Ck и xk, k = 1, 2, ..., n как решение системы

(m + 1) уравнений (4.6) затруднительно (уже при n = 4). Число неизвестных и число уравнений у неоднородной системы (4.6) должны быть равными, значит m не может быть больше 2n − 1. Если бы были известны узлы, обеспечивающие алгебраическую точность 2n − 1, то соответствующую им формулу можно было бы построить как формулу Ньютона - Котеса. Докажем две леммы, которые помогут находить такие узлы.

Лемма 1. Пусть существуют узлы xk, k = 1, 2, . . . , n, и построенная по ним квадратурная формула Ньютона-Котеса обеспечивает алгебраическую точность 2n − 1. Тогда для произвольного многочлена Pi(x) степени i 6 n − 1 справедливо равенство

∫b

ρ(x) Pi(x) ωn(x) dx = 0,

a

∏n

где ωn(x) = (x − xk) .

k=1

Доказательство. Узлы выбраны так, что квадратурная формула точна для любого многочлена степени m 6 2n − 1. Под знаком интеграла стоит многочлен Pi(x) ωn(x) степени не более 2n − 1. Формула точна для него, и

	b	n
a		∑
∫		ρ(x) Pi(x) ωn(x) dx = k=1 Ck Pi(xk) ωn(xk) = 0,

так как ωn(xk) = 0 для любого узла.

Лемма 2. Пусть для некоторого многочлена ωn(x), имеющего n различных вещественных корней на [a, b], справедливо равенство

∫b

ρ(x) Pi(x) ωn(x) dx = 0,

a

48

где Pi(x) – произвольный многочлен степени i 6 n−1. Тогда квадратурная формула Ньютона - Котеса, построенная по узлам, совпадающим с корнями многочлена ωn(x), точна для любого многочлена степени не более 2n − 1.

Доказательство. Выберем произвольный многочлен Qm(x) степени m 6

2n − 1 и представим его в виде

Qm(x) = pi(x) ωn(x) + qj(x),

где pi(x) и qj(x) многочлены степеней i 6 n − 1 и j 6 n − 1. Из условия

леммы следует, что

b	b	b
∫a ρ(x) Qm(x) dx =	∫a ρ(x) pi(x) ωn(x) dx +	∫a ρ(x) qj(x) dx =
	∫b
	= ρ(x) qj(x) dx.
	a

Если построить квадратурную формулу Ньютона– Котеса по n узлам xk

(корням многочлена ωn(x)), то она будет точна для многочлена qj степени j 6 n − 1

	b	n
a		∑
∫		ρ(x) qj(x) dx = k=1 Ck qj(xk).

В узлах xk выполняется равенство

Qm(xk) = pi(xk) ωn(xk) + qj(xk) = qj(xk), k = 1, . . . , n.

Квадратурная формула точная для произвольного многочлена Qm(x)

степени 2n − 1 построена, действительно

	b		b	n	n
a		a		∑	∑
∫		ρ(x) Qm(x) dx = ∫		ρ(x) qj(x) dx = k=1 Ck qj(xk) = k=1 Ck Qm(xk).

Если весовая ρ(x) ≡ 1, то такие формулы называются просто формулами Гаусса.

49

Итак, чтобы построить квадратурную формулу с алгебраической точностью 2n − 1, надо найти многочлен ωn(x), удовлетворяющий условию

∫b

ρ(x) Pi(x) ωn(x) dx = 0,

a

для произвольного многочлена Pi(x) степени i 6 n − 1 и имеющий n вещественных корней на [a, b]. Если весовая функция ρ(x) положительна почти всюду на [a, b], то такой многочлен существует, известен алгоритм его построения и вычисления его корней (см. стр. ??). Используя вычисленные корни ωn(x) в качестве узлов, надо построить формулу Ньютона - Котеса, которая и является формулой типа Гаусса.

Для квадратур Гаусса (ρ(x) ≡ 1) существуют таблицы узлов и коэффициентов, рассчитанные для больших значений n с большим числом знаков.

Оценка погрешности формул типа Гаусса.

Построим интерполяционный многочлен Qn(x) с n двукратными узлами, в качестве которых возьмем узлы квадратуры Гаусса xk, k = 1, 2, ..., n. Степень такого многочлена 2n − 1, и верна формула

f(x) = Qn(x) + r(x)

где r(x) = f(2n)(ξ(x)) ∏n (x − xk)2, ξ(x) (a, b). (2n)!

k=1

Применим для вычисления интеграла I(f) квадратурную формулу Гаусса и выпишем два равенства

I(f) = I(Qn) + I(r),

Sn(f) = Sn(Qn) + Sn(r).

Из этих равенств следует, что погрешность квадратуры Гаусса

Rn(f) = I(f) − Sn(f) = (I(Qn) − Sn(Qn)) + I(r) − Sn(r).

50