Schisla_1

)2n+1

Слагаемое I(Qn) − Sn(Qn) = 0, так как многочлен Qn(x) степени 2n − 1, и Sn(r) = 0, так как r(xk) = 0 для каждого узла xk. В результате для погрешности Rn(f) получается формула

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

из которой следует оценка

|

R

(f)

6

C

M2n

(

b − a

(2n)!

n

|

2

где

M

max

|

f(2n)(x)

|

,

2n = a6x6b

1

n

−1

a + b

∏

C =

ρ

(

+ t

b − a

) k=1

(t

−

d

)2 dt,

∫

2

2

k

,

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

Замечание 1. Если весовая функция ρ(x) положительна почти всюду на

[a, b], то формула Гаусса является квадратурой наивысшей алгебраической точности, то есть нельзя построить квадратурную формулу точную для любого многочлена степени m 6 2n. Например, формула не точна для

∏

n

многочлена P2n = (x − xk)2 степени 2n. Действительно,

k=1

I(P2n) = ∫b ρ(x)	n	(x − xk)2 dx > 0
	∏

ak=1

∑n

Sn(P2n) = CkP2n(xk) = 0

k=1

Замечание 2. Если весовая функция ρ(x) положительна почти всюду на

[a, b], то коэффициенты Гаусса положительны.

			n	(x−xk)2 степени (2n−
Действительно, возьмем многочлен Q(x) =
2	x		=1;k=j
			[a,kb]∏̸	Q(x) = 0	для любого
). Этот многочлен положителен для			, и

51

x = xk, кроме x = xj. Формула Гаусса точна для Q(x). Учитывая свойства этого многочлена, получим

0 < ∫	b	n	n	n
		ρ(x)Q(x) dx = i=1	Cik=1;k=j(xi − xk)2 = Cj	k=1;k=j(xj − xk)2,
a		∑	∏̸	∏̸

то есть Cj > 0.

Замечание 3. При выводе формул Чебышева с равными квадратурными коэффициентами использовались узлы, вычисленные с помощью корней полинома Чебышева. Полиномы Чебышева образуют систему ортогональ-

ных с весом ρ(x) = √ 1 многочленов, то есть

1 − x2

∫1		√		1		Pn(x) Pm(x) dx = 0, если n ̸= m.
			1	−	x2
−	1

Это означает, что формулы Чебышева являются формулами типа Гаусса.

Пример 1. Построим квадратурную формулу Гаусса по двум узлам для

∫1

вычисления xf(x) dx, где весовая функция ρ(x) = x.

0

Сначала надо построить многочлен второго порядка P2(x), ортогональный с весом ρ(x) = x любому многочлену степени m < 2, и найти его корни. Затем по двум узлам, совпадающими с найденными корнями, вы-

вести формулу Ньютона-Котеса.

Если многочлен второй степени P2 = x2 + bx + c ортогонален с весом

x многочленам P0(x) =		1 и P1(x) = x, то он ортогонален произвольному
многочлену степени m < 2. Константы b и c найдем из условий
∫1	∫1
1	1
	x P0(x) P2(x) dx = x(x2 + bx + c) dx = 0 1/4 + b/3 + c/2 = 0,
0	0
∫0	x P1(x) P2(x) dx = ∫0	x2(x2 + bx + c) dx = 0 1/5 + b/4 + c/3 = 0.

Из этих равенств вычисляем b =			−1.2, c = 0.3. Теперь корни x1 =
√		≈ .8449489743 и x2 = 0.6	√		≈ .3550510257 многочлена
0.6 + 0.1	6		− 0.1	6

52

P2(x) = x2 − 1.2x + 0.3, возьмем в качестве узлов и построим формулу Ньютона-Котеса. Полученная формула

S2(f) = .3180413817 f(x1) + .1819586182 f(x2)

является квадратурой типа Гаусса для весовой функции ρ(x) = x и x

[0, 1].

Пример 2. Построим квадратурную формулу Гаусса с одним узлом для

∫b

вычисления (x − a)2f(x) dx, где ρ(x) = (x − a)2 – весовая функция.

a

Многочлен первой степени P1(x) = x + c должен быть ортогональнен с весом ρ(x) = (x − a)2 многочлену нулевой степени P0(x) = 1. Константу c

найдем из условия

b

b

∫a

(x − a)2 P0(x) P1(x) dx = ∫a

(x − a)2(x + c) dx = 0 c = −

a + 3b

4

Корень x1 =

a + 3b

многочлена P1(x) является искомым узлом квадратур-

4

ной формулы Гаусса, а сама формула имеет вид

b

∫a

(x − a)2f(x) dx ≈ (b − a)f (

a + 3b

)

(4.7)

4

R(f)

6

M2(b − a)5

Погрешность этой квадратурной формулы |

|

160

, где

M

= max

f′′(x)

2

a6x6b |

|.

4.3Составные (обобщенные) квадратурные формулы

Каждая из формул Ньютона-Котеса и других интерполяционных квад-

ратурных формул, построенных по n узлам, может быть записана в виде

	b	n
a		∑
I(f) = ∫		ρ(x) f(x) dx ≈ Sn(f) = k=1 Ckf(xk).

53

При выводе формулы Ньютона-Котеса задаются узлы, как правило равноотстоящие, выбор неопределенных коэффициентов может обеспечить алгебраическую точность n− 1. Если неопределены и коэффициенты и узлы, то их выбор позволяет построить квадратуры типа Гаусса, обеспечивающие

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

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

формул можно записать в виде

b

f(m)(ξ(x))

n

a

∏

Rn(f) = ∫

ρ(x)

m!

k=1(x − xk)j dx,

откуда получается оценка

| 6

(

)

|

n

m!

2

m+1 ,

R

(f)

C

Mm

b − a

M

max f(m)(x)

ξ(x)

(a, b)

где

m = a6x6b

|

|,

, константа

1

n

C =

1

ρ

a + b

+ t

b − a

∏

t

d j dt.

−

2

2

| −

k|

∫

(

) k=1

Значения m и j зависят от типа формулы. Например, для формулы НьютонаКотеса m = n, j = 1, а для формулы Гаусса m = 2n, j = 2.

Из оценки погрешности видно, что она будет малой, если достаточно мала длина (b − a). Этим можно воспользоваться, так как интеграл по отрезку равен сумме интегралов по его частям. Если разбить отрезок [a, b] на

N отрезков длиной (b−a)/N с помощью точек aj = a+jH, j = 0, 1, 2, ..., N,

a0 = a, aN = b, для интеграла I(f) справедливо равенство

	b	N−1	aj+1		N−1
I(f) = ∫			∫
		ρ(x) f(x) dx = j=0		ρ(x) f(x) dx = j=0 Ij(f).
a		∑ aj			∑

На каждом отрезке [aj, aj+1] выберем квадратурную формулу Sjn(f) с

54

1/(1 + 25x2) dx (пример

коэффициентами Cjk, построенную по n узлам xjk [aj, aj+1]. Приближен-

ное значение Ij(f) можно записать в виде

a

n

H

n

Ij(f) ≈ Sjn(f) = (aj+12−

j ) k=1 Cjkf(xjk) =

( 2 ) k=1 Cjkf(xjk),

∑

∑

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

I (f)

S

C

Mjm

b − a

m+1 = C

Mjm

H

m+1 ,

−

jn(f)| 6

( 2N )

j m!

( 2 )

| j

j m!

где

M

max

|

f(m)(x)

|, а константа

C

j вычислена для отрезка

[a

, a

]

.

jm = aj6x6aj+1

j

j+1

Теперь расчетную формулу для вычисления I(f) можно записать в виде

b

N−1

I(f) = ∫ ρ(x)f(x) dx = j=0 Ij(f) ≈

a

∑

≈ Sn =

N−1 H

n

H N−1

n

Cjkf(xjk) =

Cjkf(xjk).

2

=1

2

j=0

j=0 k=1

∑ ∑k

∑ ∑

Эта формула называется составной или обобщенной квадратурной формулой. Оценка её погрешности имеет вид

N−1

N−1

Mjm

H

m+1

|Rn(f)| = |I(f) − Sn(f)| 6 j=0

|Ij(f) − Sjn(f)| 6 j=0

Cj

m!

( 2 )

.

Если

max f(m)(x)

∑

∑

M

C =

max

C

m = a6x6b |

|,

06j6N−1

j, то

R

(f)

N C

Mm

H

m+1

=

b − a

C

Mm

H

m .

m! (

2 )

m!

(

2 )

(4.8)

| n

| 6

2

Константа Mm в оценке Rn(f), как правило, неизвестна, поэтому полученная оценка погрешности не может быть использована в практическом счете. Однако, из этой оценки следует, что желание улучшить точность вычислений за счет увеличения числа m (за счет увеличения количества узлов n) требует, чтобы при возрастании m максимум производной подынтегральной функции f(m)(x) не становился очень большим. На-

пример, при использовании квадратурной формулы c симметрично распо-

∫1

ложенными узлами для вычисления интеграла

−1

55

Рунге) f(m)(x) = 0(25m) при возрастании m, и рассчитанное значение интеграла не стремится к точному. Добиваться требуемой точности расчета можно за счет уменьшения шага H.

Пример 1. Выведем составную формулу Симпсона. Квадратурная формула Симпсона (4.5) имеет вид

b

b − a

1

4

a + b

1

I(f) = f(x) dx

S

(f) =

f(a) +

f(

) +

f(b) .

∫a

≈

3

2 (3

3

2

3

)

R

(f)

M

(b − a)5

M

= max

f(4)(x)

Еe погрешность равна |

2880 , где

|

3

| 6

4

4

a6x6b

|.

Разобьем отрезок [a, b] на четное число (N = 2k) отрезков длиной H = b N− a. Для вычисления интеграла используем формулу Симпсона на каж-

дом таком отрезке [a + jH, a + (j + 1)H], j = 0, 1, ..., N − 1
				Ij(f) =							a+(∫j+1)Hf(x) dx ≈												(4.9)
											a+jH
	H	1			4							1
≈		(		f(a + jH) +			f(a + (j + 1/2)H) +												f(a + (j + 1)H)).
	2		3		3												3
Если обозначить h = H/2 и xj												= a + jh, то aj = x2j, j = 0, 1, ..., k, и
(4.9) примет вид
						1						4						1
				Ij(f) ≈ h (				f(x2j) +						f(x2j+1) +						f(x2j+2)).			(4.10)
						3							3						3
Используя (4.10), получим составную формулу
											b					k−1
										a						∑
				I(f) = ∫							f(x) dx = j=0 Ij(f) ≈
				k−1		1						4						1
				∑		(																f(x2j+2)).
				≈ S3 = h j=0					3	f(x2j) +					3	f(x2j+1) +					3

В силу того, что коэффициенты квадратурной формулы Симпсона не зависят от пределов интегрирования, можно привести подобные, и полу-

56

чить составную формулу Симпсона

h

∫a b f(x) dx ≈

k−1

(4.11)

k−1

(f(a) + f(b) + 4

∑

∑j

≈ S3 =

3

f(x2j+1) + 2

f(x2j)).

j=0

=1

Для оценки погрешности выведенной формулы получим выражение

R (f)

k−1

I (f) R (f)

M

k−1

h5

= M

(b − a)h4

,

∑j

∑

| 3

| j − j

| 6 4

90

4

180

| 6

j=0

=0

где

M

max

f(4)(x)

|.

4

= a6x6b

|

Пример 2. Выведем составную формулу Гаусса с одним узлом. Квадратурная формула Гаусса с одним узлом (4.7) имеет вид

	b
∫a	(x − a)2f(x) dx ≈ (b − a)f (	a + 3b	).
		4

Разобьем отрезок [a, b] на N отрезков длиной H = b N− a. Для вычисления интеграла на каждом таком отрезке [a + jH, a + (j + 1)H], j = 0, 1, ..., N − 1 используем формула Гаусса с одним узлом

a+(∫j+1)H

Ij(f) =

f(x) dx ≈ Hf(x1(j)),

a+jH

()

где узел x1(j) = a +	4j + 3		H.
	4

Составная квадратурная, соответствующая выведенной формуле Гаусса,

имеет вид

	b	N−1
a
		∑
∫		f(x) dx ≈ S1 = H j=1 f(x1(j)).	(4.12)

57

В отличии от формулы Симпсона, здесь не удается привести подобные, так как при каждом разбиении на N отрезков приходится вычислять значение подынтегральной функции в новых узлах x(1j).

4.4Практическая оценка погрешности (правило Рунге)

Пусть шаг H настолько мал, что значение f(m)(x) на каждом отрезке [aj, aj+1] отличается от константы на бесконечно малую порядка O(H). Тогда можно выделить главный член погрешности (4.8), представив её в виде

Rn(f) = DHm + O(Hm+1),

где константа D, как правило, неизвестна.

Информация о существовании главного члена погрешности позволяет: во первых получить способ практической оценки погрешности, во вторых выбрать начальный шаг интегрирования.

Пусть отрезок разбит на N отрезков длиной H = (b − a)/N и выбрана квадратурная формула Sn(f) = SH (f), построенная по n узлам с оценкой погрешности RH (f) = DHm + O(Hm+1). Тогда можно записать равенство

∫b

I(f) = ρ(x)f(x) dx = SH (f) + DHm + O(Hm+1). (4.13)

a

Если же разбить отрезок на 2N частей длиной H/2, то (4.13) преобразуется в равенство

∫b

I(f) = ρ(x)f(x) dx = SH=2(f) + D(H/2)m + O(Hm+1). (4.14)

a

Из этих равенств следует, что

SH (f) − SH=2(f) = D(H/2)m(1 − 2m) + O(Hm+1) и

D(H/2)m = (SH=2(f) − SH (f))/(2m − 1) + O(Hm+1).

58

Теперь (4.14) примет вид

∫b

I(f) = ρ(x)f(x) dx = SH=2(f) + SH=2(f) − SH (f) + O(Hm+1). 2m − 1

a

Отсюда следует, что с точностью до бесконечно малых порядка O(Hm+1)

погрешность

Rn(f) = I(f) − SH=2(f) ≈	SH=2(f) − SH (f)		(4.15)
	2m − 1	.

Такой способ оценки погрешности численного интегрирования называется правилом Рунге. Это правило позволяет построить алгоритм вычисления интеграла с автоматическим выбором шага, обеспечивающего заданную точность.

4.5Алгоритм вычисления интеграла с автоматическим выбором шага

∫b

Пусть требуется вычислить интеграл I(f) = ρ(x)f(x) dx с заданной

a

точностью ε, используя квадратурную формулу Sn(f) = SH (f), погрешность которой Rn(f) = DHm + O(Hm+1). Очевидно, шаг H должен быть таким, чтобы погрешность |Rn(f)| 6 ε.

Константа D неизвестна, поэтому разумно сначала выбрать шаг H =

H =

√

, при котором Rn(f) = O(ε). По шагу H вычисляется начальный

ε

m

шаг H

= (b

−

a)/k

, где

k = [(b

−

a)/H]

– целое. При таком выборе

H

,

e

0

0

e

погрешность имеет порядок Rn(f) = O(ε), и начальный шаг укладывается целое число раз на отрезке [a, b].

Далее, после вычисления приближенных значений интеграла по выбранной составной квадратурной формуле для шагов H0 и H = H0/2, проверя-

ется условие
	|SH0 (f) − SH (f)|	6 ε.	(4.16)
	2m − 1

59

Если условие (4.16) выполнено, то In(f) ≈ SH (f) с точность ε. Если же условие (4.16) не выполнено, то надо в качестве H0 взять H и повторить расчет. Такой расчет (деление шага пополам) повторяется до тех пор, пока условие (4.16) не будет выполнено. При этом, надо помнить, что правило Рунге — это приближенная практическая оценка точности квадратурной формулы, и может случиться, что делением шага пополам невозможно достичь заданную точность. Тогда, если за ограниченное число делений шага (заданное вычислителем) не удается добиться выполнения условия (4.16), требуется сообщить вычислителю, что интеграл не сходится.

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

4.6Вычисление интегралов от сильно осциллирующих функций

Пусть требуется вычислить интеграл

∫b

I(f) = f(x, ω) dx,

a

где параметр ω 1, а производная подынтегральной функции n раз (n > ω) меняет знак на [a, b] (сильно осциллирующая функция). График подобной функции, заданной на отрезке [0, 1], изображен на рис.3. Если применять квадратурную формулу, то шаг должен быть меньше расстояния между двумя переменами знаков.

60