tasks-task_35714_1

В случае многочлена четвертой степени (пять узлов) получим:

f ′(x			) =					1				(−25y								+		48y				− 36y							+16y								−	3y						) +				h4				f (5)	(ξ )	(17)
	0																	0														2							3								4
					12h																				1																							5									0
																												;
								1																																									h4									(18)
f ′(x				) =									(−3y							−10y						+18y								− 6y							+ y						) −									f (5) (ξ ) ;
		1															0								1							2						3						4
						12h																																										20									1
														1																													h4															(19)
				f ′(x					)		=							( y				−		8y				+		8y				− y					) +									f (5) (ξ ) ;
								2													0						1						3				4
													12h																													30														2
								1																																										h4								(20)
f ′(x				) =									(−y						+ 6y					−18y								+10y								+ 3y							) +									f (5)	(ξ )
		3														0							1							2						3									4
						12h																																										20									3
																												;
								1																																											h4							(21)
f ′(x				) =								(3y				−16y								+ 36y							− 48y								+ 25y									) +								f (5)	(ξ )
	4														0														2						3											4
					12h																	1																										5									4

.

Выбор оптимального шага численного дифференцирования

Оценка абсолютной погрешности численного дифференцирования складывается из остаточной погрешности, оцениваемой величиной |Rn′(x)|, и

вычислительной погрешности, определяемой приближенным заданием величин yi, i=0,1,...,n, (погрешностью округления результата пренебрегаем). Рассмотрим для определенности формулу (19).

Приближенное значение производной

f ′(x2 ) ≈

1

(y0

− 8y1 + 8y2 − y4 )

(22)

12h

имеет остаточную погрешность

∆1 =

h4 M

5

≥

h4

f (5)

(ξ)

,

M5

= max| f (5) (x)|,

30

30

[a;b]

и вычислительную погрешность согласно равенству (9) темы I

∆2

=

18∆*

=

3∆*

,

12h

2h

где ∆* - абсолютная погрешность каждого из чисел yi, i=0,1,...,n.

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

(22)-

∆(h) = ∆1 + ∆2 = h4 M5 + 3∆* . 30 2h

Для малости ∆1 необходима малость h, но при уменьшении h растет ∆2 . Из уравнения ∆′(h) = 0 получаем значение h*, при котором погрешность ∆(h) формулы (22) имеет минимальное значение.

∆′(h) =	4M		3∆*	45∆*
	305 h3	−	2h2 = 0,	h* = 5 4M5

Задача.

Функция f(x) задана таблицей своих значений, верных в написанных знаках. Найти первую производную этой функции в точках x1*=0,7 и x2*=1,0. Оценить погрешности результатов. Найти оптимальный шаг h* для каждой из формул численного дифференцирования.

31

Xi	0,5	0,6	0,7	0,8	0,9	1,0
yi	0,4794	0,5646	0,6442	0,7174	0,7833	0,8415

Решение.

Точка x1*=0,7 - центральный узел таблицы. Для вычисления f ′(0,7) в данной

задаче следует воспользоваться одной из формул (8), (14), (15), (18).

1) Воспользуемся формулой (8), обозначив x0=0,6; x1=0,7; x2=0,8. Тогда

1

f ′(0,7) ≈ 2 0,1 (0,7174 − 0,5646) = 0,764.

Остаточная погрешность результата в соответствии с формулой (11) -

∆1 = M63h2 ,

где M 3	= max| f (3) (x)|.
	[a;b]

Чтобы оценить M3, построим для данной функции таблицу конечных разностей.

xi		yi	∆yi		∆2 yi								∆3 yi					∆3 yi
0,5		0,4794	0,0852
0,6		0,5646		-0,0056
0,7		0,6442	0,0796	-0,0064							-0,0008							-0,0001
0,8		0,7174	0,0732	-0,0073							-0,0009							0,0003
0,9		0,7883	0,0659	-0,0077							-0,0006
1,0		0,8415	0,0582
						max|∆3 y						|		0,0009
			M 3 ≈				[a;b]				i		=					= 0,9;
									h3					(0,1)3
				∆1 =				0,9			(0,1)2				= 0,0015.
											6
	Вычислительная погрешность результата -
				∆2	=		∆*			=	0,00005					= 0,0005.
								h
												0,1
где ∆*	= 0,00005 - абсолютная погрешность величин yi.
			∆(h) = ∆1 + ∆2							= 0,0015 + 0,0005 = 0,002.

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

M3h

2

*

2M3h

*

∆(h) =

+ ∆h ; ∆′(h) =

−

∆

= 0,

6

6

h2

откуда

h

*

= 3

3∆*

=

3 0,5 10−4

= 0,119 .

M3

0,9

32

2) Решим теперь данную задачу с помощью формулы (14), обозначив x0=0,6; x1=0,7; x2=0,8; x3=0,9.

′

1

6 0,1 (−2 0,5646

−3 0,6442 +6 0,7174 −0,7833) = 0,7655

f

(0,7) =

∆ =

h3M

4

,

1

12

max |

∆4 y |

M4 = max | f

(4) (x) |≈

[a;b]

i

=

0,0003

= 3;

h4

[a;b]

(0,1)4

∆

=

(0,1)3 3

= 0,000025; ∆

=

2 ∆*

=

2 0,00005

= 0,001;

2

1

12

h

0,1

∆(h) = ∆1 + ∆2 = 0,001025.

Определим для этой формулы оптимальный шаг численного

дифференцирования.

h3M

2∆*

h2M

2∆*

∆(h) =

12 4

+ h ; ∆′(h) =

4 4

− h2 = 0,

h

*

= 4

8∆*

=

8 0,5 10−4

= 0,107.

M4

3

Точка x2*=1,0 является последним узлом таблицы. Для вычисления f ′(1,0)

служат формулы (9), (16), (21). Воспользуемся формулой (16), обозначив x0=0,7; x1=0,8; x2=0,9; x3=1,0.

1

f ′(10,) = 6 01,(−2 0,6442 +9 0,7174 −18 0,7833+11 08415,) =054217,;

M4h

3

max|∆4 y |

0,0003

∆1 =

≈

[a;b]

i

=

=0,00075;

4

4 h

4

01,

∆2 =

40 ∆*

=

20 ∆*

=

20 0,00005 .

6h

3h

3 0,1

Определим соответствующий данной формуле оптимальный шаг таблицы.

∆(h) = ∆1 + ∆2

=

M4h3

+

20∆*

′

=

3M4h2

−

20∆*

= 0,

4

3h

; ∆ (h)

4

3h2

h

*

80∆*

80 0,5 10−4

= 4 9M4

= 4

9 3

= 0,110.

Задача B

Пользуясь таблицей задачи Б2, вычислить первую производную заданной функции в точке x* и оценить погрешность результата. Определить оптимальный шаг таблицы для выбранной формулы численного дифференцирования.

1. x*=1,1	2. X*=1,2		3. X*=1,3		4. X*=2,0
5. x*=2,2	6.	X*=0,50	7. X*=0,52		8. X*=0,56
9. x*=0,60	10.	X*=0,61	11.	X*=1080	12.	X*=1090
13. x*=1100	14.	X*=1110	15.	X*=1120	16.	X*=2,70
17. x*=2,74	18.	X*=2,76	19.	X*=2,80	20.	X*=2,84
21. x*=0,7	22.	X*=0,9	23.	X*=1,1	24.	X*=1,3
25. x*=1,5.

33

4.4. Численное интегрирование

Постановка задачи. Пусть требуется вычислить интеграл

b
J = ∫ f (x)dx.	(1)

a

Если функция f(x) является непрерывной на отрезке [a;b], то интеграл (1) существует и может быть вычислен по формуле Ньютона-Лейбница

b
J = ∫ f (x)dx = F(b) − F(a).	(2)

a

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

Определение 1. Приближенное равенство

b	N
J = ∫ f (x)dx ≈ (b − a)∑ Ai		f (xi ) = J N	(3)
a	i =1			узлами xi [a;b] и
называется	квадратурной	формулой,	определяемой
коэффициентами Ai.
Величина	RN ( f ) = J − J N		(4)
называется остаточным членом квадратурной формулы.
В зависимости от способа		задания подынтегральной		функции f(x) будем

рассматривать два различных в смысле реализации случая численного интегрирования.

Задача 1. На отрезке [a;b] в узлах xi заданы значения fi некоторой f, принадлежащей определенному классу F. Требуется приближенно вычислить интеграл (1) и оценить погрешность полученного значения.

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

Задача 2. На отрезке [a;b] функция f(x) задана в виде аналитического выражения. Требуется вычислить интеграл (1) с заданной предельно допустимой погрешностьюε .

Рассмотрим алгоритмы решения задач 1 и 2.

Алгоритм решения задачи 1.

1. Выбирают конкретную квадратурную формулу (3) и вычисляют JN. Если значения функции f(xi) заданы приближенно, то фактически вычисляют лишь приближенное значение J N для точного JN.

2.Приближенно принимают, что J ≈ J N .

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

∆1 ≥| J − JN |=|RN ( f )|. 4. Определяют погрешность вычисления J N

∆2 ≥| J N − J N |,

по погрешностям приближенных значений f(xi).

34

5. Находят полную абсолютную погрешность приближенного значения J N :

| J − J N |≤ ∆1 + ∆2 = ∆

6. Получают решение задачи в виде

J = J N ± ∆ .

Алгоритм решения задачи 2.

1. Представляют ε в виде суммы трех неотрицательных слагаемых

ε = ε1 + ε2 + ε3 ,

где ε1 - предельно допустимая погрешность метода; ε2 - пре-дельно допустимая погрешность вычисления J N ; ε3 - предельно допустимая погрешность округления результата.

2.Выбирают N в квадратурной формуле так, чтобы выполнялось неравенство

∆1 =| J − J N |=|RN ( f )|≤ ε1 .

3.Вычисляют f(xi) с такой точностью, чтобы при подсчете JN по формуле (3)

обеспечить выполнения неравенства

∆2 =| J N − J N |≤ ε2 .

Для этого, очевидно, достаточно вычислить все f(xi) с абсолютной погрешностью

εT

=

ε2

.

N

(b − a)∑| Ai |

i =1

ε3

≠ 0 )

4.

Найденную в

п.3

величину

N

округляют

(если

J

==

с предельно допустимой погрешностью

ε3 до величины J N .

5. Получают решение задачи в виде

==

J = J N

± ε .

Построение простейших квадратурных формул

Формула прямоугольников. Допустим, что

f (x) C2 [a;b].

Отрезок

[a;b]

_______

разделим

на N равных частичных отрезков [xi-1;xi], где xi=a+ih; i == 0, n

− 1; xN=b;

h =

b − a

.

N

Тогда

xi

b

N

∫ f (x)dx = ∑

∫ f (x)dx .

(5)

a

i=1 xi −1

Обозначим среднюю точку отрезка [xi-1;xi] через

ξ =

xi −1

+ xi

.

(6)

i

2

Запишем для функции f(x)

на каждом из отрезков [xi-1;xi] формулу Тейлора с

остаточным членом в форме Лагранжа

(7)

f

(x) = f (ξi ) +(x −ξi ) f ′(ξi ) +

(x −ξi )2

f ′′(ηi );

ηi (xi−1; xi )

2!

Подставим в правую часть соотношения (5) вместо f(x) ее представление (7)

35

b

N

xi

(x −ξi )

2

∫ f (x)dx = ∑ ∫[ f (ξi ) +(x −ξi ) f ′(ξi ) +

f ′′(ηi )]dx =

2!

a

i=1 xi−1

N

xi

xi

xi

(x −ξi

)

2

∫

dx + f ′(ξi )

∫

(x −ξi )dx +

∫

f ′′(ηi )dx]

(8)

i=1

2

= ∑[ f (ξi )

xi−1

xi−1

xi−1

xi

(x − ξi )2

Используя для вычисления

∫

f ′′(ηi )dx

вторую теорему о среднем

2

xi −1

xi

значении функции и, учитывая, что

∫(x − ξi )dx = 0 , получим, что

xi −1

b

N

h

3

N

∫ f (x)dx = h∑ f (ξi ) +

∑ f ′′(

ηi

);

ηi

(xi−1, xi )

(9)

a

i=1

24 i=1

В силу непрерывности

f ′′(x)

существует такая точка η (a;b) , что

N

(10)

∑ f ′′(ηi ) = Nf ′′(η) .

i=1

Используя (10), получаем

b

N

h

3

∫ f (x)dx = h∑ f (ξi ) +

Nf ′′(η)

b − a

a

i=1

24

или, так как h =

,

N

b

N

1

b − a

2

′′

∫ f (x)dx = (b −a)∑

f (ξi ) +

h

f

(11)

N

24

(η)

a

i=1

Приближенное равенство

b

N

1

∫ f (x)dx ≈ (b − a)∑

f (ξi ) = J Nпр

(12)

N

a

i =1

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

узлами ξi [a;

RN ( f

	b] и коэффициентами A =			1		. Величина
		i		N
	b
		b − a
	) = ∫ f (x)dx − J Nпр =		h2	f	′′(η)		(13)
		24
	a

является остаточным членом формулы прямоугольников (12).

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

|RN ( f )|≤	b − a	h2	M 2 = ∆	1 ,	(14)
24
где			M 2	= max| f ′′(x)| .
				[a;b]

Выражения для остаточного члена (13) и остаточной погрешности (14) показывают, что формула прямоугольников (12) является точной для любой линейной функции, так как вторая производная такой функции равна нулю, и, следовательно, ∆1 = 0 .

Оценим вычислительную погрешность ∆2 формулы прямоугольников,

которая возникает за счет приближенного вычисления значений функции f(x) в узлах

ξi .

36

Пусть, например, значения f( ξi ) в формуле (12) вычислены с одинаковой

абсолютной погрешностью ∆* ,

тогда

(15)

N

1

∆2 =| JN −

J

N |= (b −a)∑

∆* = (b −a)∆*

i=1

N

Пример 1. Вычислить с помощью формулы прямоугольников

1

dx

∫0

с точностью ε

= 10-2.

1 + x

Применяя алгоритм решения задачи 2, представим суммарную погрешность

ε в виде суммы трех слагаемых.

Выберем h из условия

ε =0,01=0,009+0,0005+0,0005.

h2 (b − a)M 2

∆

=

≤ 0,009 .

1

24

Так как M

= max| f ′′(x)|= max

2

= 2 и (b-a)=1, то h ≤ 0,009 12 < 0,32

2

(1 + x)3

[0;1]

[0;1]

и, следовательно,

N =

b − a

> 3,1, т.е. N=4, h=0,25, ∆1 = 0,0052 .

h

Составим таблицу значений функции 1/1+x

с тремя знаками после запятой,

так как ∆2 = ∆* (b − a) ≤ 0,0005,

∆* = 0,0005 .

ξi

0,125

0,375

0,625

0,875

f (ξi )

0,889

0,727

0,615

0,533

Используя формулу (12), получаем

J пр4 = 0,25 (0,899 +0,727 +0,615 +0,533) = 0,25 2,764 = 0,691.

Так как в данном случае погрешность округления равна ε3 = 0,0005,, то

получим

J = 0,691± (0,0052 ± 0,0005) = 0,691± 0,0057 .

Формула трапеций. Предположим, что

f (x) C2 [a;b]. Разделим отрезок [a;b]

на N равных частей, тогда

b

N xi

∫ f (x)dx = ∑ ∫ f (x)dx ,

(16)

b −a

a

i=1 xi−1

где

x

= a +ih;i =

= b; h =

.

0, N

−1; x

N

i

N

Заменим функцию f(x) на каждом из отрезков [xi-1,xi] первой

интерполяционной формулой Ньютона первой степени

f (x) = f (x

) +

x −xi−1

( f (x ) − f (x

)) +

f ′′(ηi )

(x −x

)(x −x ),

i−1

h

i

i−1

2!

i−1

i

(17)

ηi (xi−1; xi ).

Подставляя формулу (17) в правую часть (16), интегрируя и используя вторую теорему о среднем значении функции, получим

b	N	f (x	) + f (x )		h3 N
∫ f (x)dx = ∑h		i−1	i	−		∑ f ′′(ηi );ηi (xi−1; xi ).(18)
			2
a	i=1				12 i=1

37

В силу (10) получаем:

b

f (x0 ) + f

(xN )

N −1

h2

′′

∫ f (x)dx = h

+

∑ f (xi )

−

(b − a) f

2

12

(η).(19)

a

i=1

Приближенное равенство

b

b −a f (x ) +

f (x

N

)

N −1

J = ∫ f (x)dx ≈

0

+

∑ f (xi )

= JNТР

(20)

N

2

a

i=1

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

ТР

h2

′′

RN ( f ) = J − JN

= −

(b −a) f

(21)

12

(η)

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

	RN ( f )		≤	b − a	h2 M 2 = ∆1.	(22)
			12
Формула трапеций, как и формула прямоугольников, является точной для любой

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

∆2 = (b −a)∆ .

(23)

Так как остаточные члены формул прямоугольников и трапеций (13) и (21) имеют противоположные знаки, то формулы (12) и (20) дают двухстороннее приближение для интеграла (1), то есть

J Nпр < J < J NТР , если f ˝(x) > 0,

J NТР < J < J Nпр , если f ˝(x) < 0.

В таком случае можно принять, что

J Nпр + J NТР

~

J ≈

=

J

,

тогда

2

~

J Nпр − J NТР

J − J

<

,

2

т.е. погрешность выражается через приближенные значения интегралов.

Пример 2. Вычислить ∫1

dx

по формуле трапеций, полагая N=4; оценить

0 1+ x

~

полную погрешность результата. Учитывая результаты примера 1, найти

по

J

формуле (24) и оценку (25).

Применяя алгоритм решения задачи 1, находим:

∆1

=

h2

(b − a)M 2 =

0,25

2

1 2 = 0.0104 .

12

12

Составим таблицу значений функции 1/(1+x) с тремя знаками после запятой

∆2 = ∆* (b −a) = ∆* ≤ 0,5 10−5 .

xi

0,00

0,25

0,50

0,75

1,00

f (xi )

1,000

0,800

0,667

0,571

0,500

1

1,000 + 0,500

TP

J4

=

+ 0,800 + 0,667 +

0,571

= 0,697 .

4

2

Суммарная погрешность равна

∆1 + ∆2

= 0,0104 + 0,0005 ≈ 0,011.

(24)

(25)

38

		Если округлить результат до двух знаков, то
ε3	= 0,70 −0,697 = 0,003									и
J = 0,70 ± (0,0109 + 0,003)= 0,70 ± 0,014 .
		Используя формулы (24) и (25) и результаты примера 1, получим
0,691 < J < 0,697 ;
~		0,691+			0,697
J	=						= 0,694			;
				2
		~		0,697 −0,691
J − J			<						= 0,003 ;
					2
J = 0,694 ± 0,003.

Формула Симпсона. Предположим, что f (x) C4 [a,в]. Разделим отрезок [a,в] на

N = 2k	равных частей, тогда			к−1 x2i +2
	а
	∫ f (x)dx =			∑	∫ f (x)dx ,					(26)
	в			i =0	x2i
где x	= а+ih ; i =			x	= в; h =	в − а	=	в − а	.
		0,2к −1 ;
i				2k		N		2к
										] длиной 2h по формуле
Заменим функцию f (x) на каждом из отрезков [x2i , x2i +2

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

	в	h		k	k −1		C
	J = ∫f (x)dx ≈		f (x0 ) + f (x2k ) +4∑f (x2i−1) +2∑f (x2i )				= J2k	(27)
		3
	а			i=1	i=1

	с остаточным членом
R	( f ) = J − J C	= −	h4	(в − а) f (IV )(η), η (а,в).	(28)
N	2k	180

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

R

( f

в − а

h4M

= ∆

,

(29)

N

)

≤

4

1

180

где

M

4

= max

f IV

(х)

.

[a,в]

Вычислительная погрешность формулы Симпсона равна

∆2 = ∆* (в − а).

(30)

Из выражения для остаточного члена формулы Симпсона следует, что она точна для многочленов третьей степени.

1

dx

Пример 3. Вычислить ∫

по формуле Симпсона с точностью ε = 0,0001.

+ x

0 1

Применяя алгоритм решения задачи II, представим суммарную погрешность ε в

виде суммы трех слагаемых ε = 0,0001 = 0,000045 + 0,000005 + 0,00005.

Выберем h из условия

∆ =

h4 (в − а)

M

≤ 0,000045.

4

1

180

Так как M 4 = max

f (4)(x)

= max

24

= 24, то

[0.1]

[0.1] (1+ x)5

39

h ≤ 4				0,45 10−4 180 = 4											3,375 10−4				≈ 0,135 и, следовательно,
								24
N		=		в − а					≥		10			≈ 7,38 ≈ 8.
										1.35
					h																	0,1254
Таким образом,														N = 8 ,		h = 0,125			и ∆	=			24 = 0,0000325.
																				1		180
Составим таблицу значений функций
																				11+x		с пятью знаками после запятой
∆2			= ∆* (b − a) = ∆* ≤ 0,5 10−5
xi								0.000					0.125			0.250			0.375			0.500		0.625	0.750
f (xi )								1.000					0.8888			0.800			0.7272			0.6666		0.6153	0.5714
													9						7			7		8	3
xi								0.875						1.000
f (xi )								0.5333						0.500
										3
Используя формулу (27), получаем:
1			dx					0,125
J = ∫						≈						[1,5+4(0,88889+0,72727+0,61538+0,53333) +2(0,8+0,66667+0,57143)]=
		1+x						3
0
=	0,125						[1,5 + 4 2,7648 + 2,03810 2 ] = 0,6931533 .
		3

Округляя полученный результат, получим

J = 0,69315 ± (0,0000325 + 0,000005 + 0,0000033) = 0,69315 ± 0,000041.

Правило Рунге практической оценки погрешности квадратурных формул.

Пусть

f (x) С4 [а, в]

и интеграл (1) вычисляется по формуле

прямоугольников. Наличие у

f (x) производных

′′′

f

IV

(x) позволяет при

f (x) и

выводе формулы прямоугольников (7)-(13) получить следующее полезное

соотношение

а

J = ∫ f (x)dx = J Nnp + ch2 + O(h4 ),

(31)

в

1

в

где

с =

∫ f ′′(x)dx

(32)

24

а

-

постоянная, не зависящая от h . Величина

ch2 называется главной частью

погрешности формулы прямоугольников.

Если

f (x) C4[а,в]

, то справедливо аналогичное соотношение и для формулы

трапеций

J = J TP + c h2 +O(h4 ),

(33)

N

1

1

в

′′

где

с1 = −

f

(34)

12

(x)dx

∫а

не зависит от h .

40