Учебное пособие_Численные методы -26 -дек_2014 (97-2003)

x0

2h

x0

2h

x0 2h

R(h)

f (x)dx

(x)dx

( f (x) (x))dx

x0

x0

x0

6

x0

2h

1

f '''(ε)(x x

)(x x )(x x

)dx

0

1

2

x

0

6

x0

h

6

x0

2h

1

f '''(ε)(x x

)(x x )(x x

)dx

1

f '''(ε)(x x

)(x

x )(x x

)dx.

0

1

2

0

1

2

x0

x0 h

К каждому из этих двух интегралов уже можно применить теорему о

среднем,

поскольку

f ' ' '

непрерывна на

[x0 , x0

2h]

и функция

(x x0 )(x x1 )(x x2 ) неотрицательна на первом интервале интегрирования и

неположительна на втором (то есть не меняет знака на каждом из этих интервалов). Поэтому

				1					x																					1						x
										1																											2
R(h)					f	'''(η )					(x x					0	)(x x )(x x								2	)dx					f '''(η			2	)			(x x
				6			1														1									6
									x	0																										x
																																					1
													4	x						2			2	x																4
	1						(x x					)		1				h			(x x	)		1				1							(x x				)
			f '''(η )																														)
											1										1									f '''(η								1
																																2
	6				1			4													2							6										4
														x0										x0
	h	4																	h		4
			( f	'''(η ) f						'''(η					))						f ''''(η)(η						η			)
														2												1			2
	24					1													24
(мы воспользовались теоремой о среднем, поскольку
функция;			η1	η η2 ).

0	)(x
	x
	2
	x
	1

f

' ' '

x )(x x				2	)dx
		1
h	2	(x x	)	2	x	2
		1
		2			1
					x

- непрерывная

Дифференцируя

R(h)

дважды и применяя

получим для	R(h) другое выражение:
	R(h)	h5 f ''''(η)	, где η [x
	90

затем теорему о среднем,

0	, x	0	2h].

Из обеих оценок для

R(h)

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

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

b	(b
f (x)dx
a

Если отрезок [a, b]

a)	a b				h5			b a
	( f (a) 4 f			f (b))		f ''''(η) ,	h		.
6		2			90			2

интегрирования слишком велик, то его разбивают

на 2n равных частей (полагая h b a ),

2n

отрезков a, a 2h , a 2h, a 4h ,..., b 2h,

а именно: запишем формулу Симпсона

после чего к каждой паре соседних

b применяют формулу Симпсона, в общем виде

b	h
f (x)dx		[ y		y2n	4( y1 y3 ... y2n 1 ) 2( y2			y4	... y2n 2 )]	(6.4)
	3		0
a
					h	b a	.			(6.5)

2n

Погрешность формулы Симпсона

171

			(b a)	5	f	''''(η)
	Rn							,	η [a,b]
			180 (2n)				4

6.7. Правило Рунге оценки погрешности

(6.6)

Идея метода состоит в том, чтобы, организовав вычисления значений интеграла по нескольким семействам (множествам) узлов, затем сравнить

результаты вычислений	и получить оценку погрешности.	Наиболее
удобное правило связано	с вычислением интеграла дважды:	Ln [ f ], L2n [ f ] ,

где

L	[ f ]
n

b

– приближенное значение f (x)dx при разбиении отрезка

a

интегрирования на n частей.

Правило Рунге оценки погрешности
R		[ f ]	L	2n	[ f
	2n				2

] L		[ f
	n
p	1

]

,

где p – порядок погрешности квадратурной формулы.

Для метода левых прямоугольников p 1,			правых прямоугольников
p 1, средних прямоугольников	p 2	, трапеций	p 2 , Симпсона p 4 .
После подсчета величин Ln	[ f ] и	L2n [ f ] , кроме оценки погрешности по

правилу Рунге можно также дополнительно уточнить приближенное

значение интеграла. Величина			L					[ f ]
*		[ f ]		2n	[ f ] L
							n
L	2n				2	p	1

называется уточненным (или экстраполированным) по Ричардсону

значением искомого интеграла.

Для практического вычисления интеграла	L[ f ] с заданной точностью
выбирается некоторое начальное число n	разбиений отрезка a, b и
вычисляются величины Ln [ f ] и L2n [ f ] . Если	R2n [ f ] ε , то с точностью

полагают	Ln		[ f ] L2n [ f ] . В противном случае вычисляют значение				L4n [ f ]	и
сравнивают		R		4n	[ f ]	и , и т.д.

6.8.Вычисление интеграла с заданной точностью

спомощью остаточного члена

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

Чтобы найти количество отрезков разбиения п, требующееся для вычисления интеграла с заданной точностью , нужно решить уравнение

172

относительно n: Rn ε , где Rn – остаточный член; а – необходимая точность вычисления. После этого вычислить интеграл, используя n точек разбиения.

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

где

		(b a)	2
	Rn		f '(c) ,
		2n
	f '(c) имеет максимальное значение на [a,b] .

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

где

		(b a)	2
Rn				f '(c)	,
		2n
f '(c) имеет максимальное значение на				[a,b] .

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

где

где

где

(b a)

3

Rn

f

,

24n

2

(c)

f

[a,b] .

(c) имеет максимальное значение на

Остаточный член для метода трапеций

(b a)

3

Rn

f

,

12n

2

(c)

f (c) имеет максимальное значение на

[a,b]

.

Остаточный член для метода Симпсона

(b a)

5

Rn

f

( IV )

(c)

,

180n

4

f

[a,b] .

( IV )

(c) имеет максимальное значение на

6.9. Метод Монте-Карло. Реализация простого метода

Рис. 6.6. Метод Монте-Карло

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

173

случайное число из интервала интегрирования. Вследствие случайности узла xi погрешность вычисления также будет носить случайный характер. Проведя n вычислений со случайными узлами xi , усредним результат, который и принимаем за приближенное значение интеграла

b

a

f (x)dx

	b a
		n
	n
		i 1

f

(x	)
i

.

Расчет методом Монте-Карло с заданной точностью

Отличие данного метода заключается в том, что мы увеличиваем количество случайных узлов n (например, в 10 раз), пока не достигнем определенной точности результата (разность между предыдущим и текущим результатом вычисления).

6.10. Квадратуры Гаусса

Рассмотрим функцию y f (t) ,		заданную
Поставим задачу: как нужно подобрать точки					t1 , t2
A1 , A2 , A3 ,..., An , чтобы квадратурная формула
1		n
	f (t)dt	i	i
		A	f (t	),
1		i 1

на

, t	3	,...,

промежутке	1,1 .
tn и коэффициенты
	(6.7)

была точной для всех полиномов				f (t) наивысшей возможной степени			N .
Так как в нашем распоряжении имеется					2n	постоянных	xi	и
Ai (i 1,2,3,..., n) , а полином степени				2n 1 определяется		2n коэффициентами,
то эта наивысшая степень в общем случае равна N 2n 1.
Для обеспечения равенства (6.7) необходимо и достаточно, чтобы оно
было верным при f (x) 1, t, t	2	,..., t	2n 1	.

Действительно, полагая

1				n
					i	i
	t	k	dt		A t	k
1				i 1

k 0,1,2,...,2n 1 и	2n 1
	f (t) Ck t k ,
	k 0

будем иметь

1	2n 1	1	2n 1	n	n	2n 1	n
f (t)dt Ck t k dt Ck Ai tik Ai Ck tik Ai f (ti ).								(6.8)
1	k 0	1	k 0	i 1	i 1	k 0	i 1

Таким образом, учитывая соотношения

174

1				( 1)	k 1		2		k четное,
			1					,
t	k	dt
						k 1
				k 1
1								,	k нечетное,
						0

заключаем, что для решения поставленной задачи достаточно определить ti и Ai из системы 2n уравнений

	n
	i	2,
	A
	i 1
	n
	A t	i	0,
i
i 1
...
	n				2
	i	i				,
	A t	2n 2
i 1					2n 1
	n	i
	i
	A t	2n 1		0.
i 1

(6.9)

Система (6.9) нелинейная и исследование

чтобы упростить решение в качестве узлов											ti
Лежандра n-й степени:
		1		d	n		x		n
P (x)											(n
								2	1
n	2	n	n! dx			n

ее громоздко. Поэтому берут корни полинома

0,1, 2,...) .

(6.10)

	Формула	(6.7), где	ti	–	нули полинома Лежандра	Pn (t)	и
Ai	(i 1, 2, 3..., n) определяются из системы (6.9), называется квадратурной
формулой Гаусса.
	Рассмотрим	использование			квадратурной формулы Гаусса		для

b

вычисления общего интеграла

f (x)dx

. Делаем замену переменной

a

x

a b

b a

t .

2

2

Получим

b

b a

1

b a

b a

f (x)dx

.

f

t dt

a

2

1

2

2

(6.11)

(6.12)

Применяя к интегралу (6.12) квадратурную формулу Гаусса (6.7), будем иметь

b

b a

n

f (x)dx

Ai f (xi ) ,

(6.13)

a

2

i 1

где x

a b

b a

t

(i 1, 2,..., n)

;t

– нули полинома Лежандра

P (t) .

i

i

i

2

2

n

175

Таблица 6.2

Таблица узлов и коэффициентов формулы Гаусса

n	ti	Ai
1	t1=0	A1=2
2	t1=0,577350269	A1=1
	t2=-0,577350269	A2=1
3	t1=0,774596669	A1=0,555555556
	t2=0	A2=0,888888889
	t3=0,774596669	A3=0,555555556
4	t1=-0,861136312	A1=0,347854845
	t2=-0,339981044	A2=0,652145155
	t3=0,339981044	A3=0,652145155
	t4=0,861136312	A4=0,347854845
5	t1=-0,906179846	A1=0,236929885
	t2=-0,538469319	A2=0,478628670
	t3=0	A3=0,56888889
	t4=0,538469319	A4=0,478628670
	t5=0,906179846	A5=0,236929885
6	t1=-0,932469514	A1=0,171324492
	t2=-0,661209386	A2=0,360761573
	t3=-0,238619186	A3=0,467913934
	t4=0,238619186	A4=0,467913934
	t5=0,661209386	A5=0,360761573
	t6=0,932469514	A6=0,171324492
7	t1=-0,949107912	A1=0,129484966
	t2=-0,741531185	A2=0,279705391
	t3=-0,405845151	A3=0,381830051
	t4=0	A4=0,417959184
	t5=0,405845151	A5=0,381830051
	t6=0,741531185	A6=0,279705391
	t7=0,949107912	A7=0,129484966

6.11. Примеры

Пример 6.1
Вычислить интеграл	I

остаточный член формулы,

по формуле Симпсона при n = 8, оценить составив таблицу конечных разностей:

	1,6	sin 2x 2,1
	I				dx.
		x	2	1
	1,2

h

Решение.

b a n 1,6

Согласно					условию,							n		=		8,		поэтому
1,2 8 0,05. Вычислительная формула имеет вид
I	h	y		4 y 2 y		4 y		2 y		4 y		2 y		4 y		y	,
			0		2		3		4		5		6		7
3				1												8

176

где

y

y(x

)

sin 2x

i

2,1

,

x

1,2 ih

(i 0, 1,..., 8).

i

2

i

i

x

1

i

Вычисление значений функции, а также сложение значений функции, имеющих одинаковые коэффициенты в формуле, производим в табл. 6.3.

i

0

1

2

3

4

5

6

7

8

x	i

1,20

1,25

1,30

1,35

1,40

1,45

1,50

1,55

1,60

2x	i	2,1

0,30

0,40

0,50

0,60

0,70

0,80

0,90

1,00

1,10

Вычисление значений функции

sin 2x		2,1	x	2	1	y		, y		y	, y		, y		, y
	i						0		8			3		5		7
				i						1
0,2955			2,44			0,1211
0,3894			2,5625								0,1520
0,4794			2,69
0,5646			2,8225								0,2000
0,6442			2,96
0,7174			3,1024								0,2312
0,7833			3,25
0,8414			3,4025								0,2473
0,8912			3,56			0,2503
						0,3714					0,8305

Таблица 6.3

y	2	, y	4	, y	6

0,1782

0,2176

0,2410

0,6368

Следовательно,

	I		,05	0,3714	4		0,8305	2 0,6368		0,05	4,9670		0,88278 0,8828.
			3							3
	Для оценки остаточного члена формулы составим таблицу конечных
	разностей функций до разностей четвертого порядка (табл. 6.4).
															Таблица 6.4
							Таблица конечных разностей
	i			yi			yi			2			3		4
									yi				yi		yi
	0		0,1211				0,0309		-0,0047				0,0003		-0,0001
	1		0,1520				0,0262		-0,0044				0,0002		0,0000
	2		0,1782				0,0218		-0,0042				0,0002		0,0000
	3		0,2000				0,0176		-0,0040				0,0002		0,0001
	4		0,2176				0,0136		-0,0038				0,0003		-0,0001
	5		0,2312				0,0098		-0,0035				0,0002
	6		0,2410				0,0063		-0,0033
	7		0,2473				0,0030
	8		0,2503

Так как

max yi	0,0001,		то остаточный член формулы
4
	R	b a max 4 yi			0,4 0,0001		0,0000003.
	ост		180		180

177

Поскольку вычисления производились цифрами, то остаточный член на результат полученные десятичные знаки верны. Ответ: I

с четырьмя значащими не влияет, значит, все

0,8828.

Пример 6.2

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

		1,3		dx
	I
			2x	2	0,3
		0,7

.

Решение. Для достижения заданной степени точности необходимо определить значение n так, чтобы

b a
	3
12n	2	M	2	0,0005.

(6.14)

Здесь a 0,7; b 1,3; M 2					max	f	(x) ,		где f (x) 1						2x		0,3.
																2
					[0,7;1,3]
Находим
	2x							8x		2		0,6							8 1,3	2		0,6
f (x)			3	,	f (x)								5	,	max		f (x)					5
		2									2				[0,7;1,3]						2
	2x		0,3					2x				0,3							2 0,7			0,3

6,98

.

Положим M 2 7, тогда неравенство (6.14) примет вид

откуда n	2	252,	то есть n>16, возьмем n = 20.

Вычисление интеграла производим по формуле

	0,6	3
				7	0,0005,
	12n		2

I

где h b a n 0,620 0,003,

Все расчеты приведены вычисленные суммы.

y

0

y

20

y

y

y

... y

h

,

2

3

20

2

1

y

i

y(x

) 1

2x2

0,3 , x

i

0,7 ih

(i 0, 1,

2, . . . , 20) .

i

в

табл. 6.5.

Последняя строчка

содержит

Таблица 6.5

i

1

0

1

2

3

4

5

6

x	i

2

0,7

0,73

0,76

0,79

0,82

0,85

0,88

Результаты расчетов

2	2	2x			0,3
xi	2xi 0,3		2
			i
3	4			5
0,49	1,28	1,1314
0,5329	1,3658	1,1686
0,5576	1,4552	1,2063
0,6241	1,5482	1,2443
0,6724	1,6448	1,2825
0,7225	1,7450	1,3210
0,7744	1,8488	1,3597

y0 , y20	y1 , y2 , y3 ,... y20
6	7
0,88386
	0,85572
	0,82898
	0,80366
	0,77973
	0,75700
	0,73546

178

Окончание табл. 6.5

1	2	3	4	5	6	7
7	0,91	0,8281	1,9562	1,3986		0,71501
8	0,94	0,8836	2,0672	1,4378		0,69551
9	0,97	0,9409	2,1818	1,4771		0,67700
10	1,00	1,0000	2,3000	1,5166		0,65937
11	1,03	1,0609	2,4218	1,5562		0,64259
12	1,06	1,1236	2,5472	1,5960		0,62657
13	1,09	1,1881	2,6762	1,6356		0,61140
14	1,12	1,2544	2,8088	1,6759		0,59669
15	1,15	1,3225	2,9450	1,7161		0,58272
16	1,18	1,3924	3,0848	1,7564		0,56935
17	1,21	1,4641	3,2282	1,7967		0,55658
18	1,24	1,5376	3,3752	1,8372		0,54431
19	1,27	1,6129	3,5258	1,8777		0,53253
20	1,30	1,6900	3,6800	1,9817	0,52129	0,
					1,40515	12,77022

Таким образом, I

Ответ: I 0,404.

1,40515
0,03		12,77022
	2

0,40418

0,404

.

Пример 6.3
Найти точки xi	и коэффициенты Ai
	1
вычислить интеграл		(x 1)dx .
	1

для

n 4

для квадратуры Гаусса и

Решение. Полином Лежандра четвертой степени

P (t)						1				d 4			(x 2		1)4 x 4 3(x 2 1)x 2											3	(x 2 1)2 4x 4 3x 2						3	x 4	6	x 2	3
4					24 4! dx4																				8							8		8			8
	1	(35x 4					30x 2						3).
	8
Приравнивая этот полином к нулю, находим корни:
t				15					120				, t					15 120			,	t			15			120	, t			15 120		,
														2									3							4
1							35													35					35							35
t1 0,861136312,t2															0,339981044,								t3 0,339981044,									t4 0,861136312.
Для определения коэффициентов A1, A2, A3, A4																																имеем систему
A1 A2 A3 A4 2,
			t2 A2						t3 A3 t						4 A4 0,
t1 A1
	2			t		2			t		2	A3 t				2			2
t1 A1						2	A2				3					4	A4			,
																			3
	3 A			t		3 A			t		3 A t					3 A 0.
t
	1	1				2	2		3				3		4		4
																								179

Решая эту систему, получаем

	A 0,347854845,					A 0,652145155,						A 0,652145155,	A 0,347854845,
	1					2						3	4
1
	(x 1)dx A (t	1	1) A	(t	2	1) A	(t	3	1) A	(t	4	1) 0,347854845(-0,861136312 1)
	1		2			3			4
1
0,652145155(-0,339981044 1) 0,652145155(0,339981044 1) 0,347854845(-0,861136312 1)
2(0,347854845 0,652145155) 2 1 2.

6.12. Примеры применения численного интегрирования в экономике

Интегральное исчисление – широко применяемый для экономического анализа математический аппарат [21, 22]. В экономическом анализе часто приходится решать следующую задачу: дана функция f(x), требуется найти функцию F(x) такую, что F’(x)=f(x) (например, найти суммарные издержки, зная предельные). Для решения такой задачи служит операция интегрирования. Вообще говоря, в экономических задачах переменные меняются дискретно. Для использования определенного интеграла нужно составить некоторую идеализированную модель, предполагающую непрерывное изменение зависимых переменных (функций) и независимых переменных (аргумента).

Модель экономического роста, предложенная Е.Д. Домаром

Основные допущения этой модели.				I t влияет
1.	Всякое изменение величины скорости денежного потока
	как на совокупный спрос, так и на совокупное предложение объема
	производства.	Y t
2.	Скорость изменения величины спроса		пропорциональна
	производной скорости денежного потока с коэффициентом
	пропорциональности k 1/ s , где s –	предельная		величина
	накопления. Это предположение можно записать в виде уравнения

dY dt

1 dI s dt

.

(6.15)

3.Экономический потенциал k (т.е. величина стоимости товаров, которые можно произвести) пропорционален объему оборотных средств К с коэффициентом пропорциональности , т.е. k= K. Дифференцируя по t, получим

	dk		dK	I .	(6.16)
	dt		dt
		180