Вычислительная математика лекции

При достаточно малом h бесконечно малыми в(1) и (2) можно

пренебречь.

В этом случае, вычитая из (1) (2), получим

Ih/2 - Ih ≈ (1/2k)Chk(2k-1) ≈ (I-Ih)(2k-1)(1/2k)≈( I - Ih/2 ) (2k-1).

I-I	h/2		I h/2			I h				(3)
				2k		1
I - Ih					I h/2		I h
									2k	(4)
					2k
							1

k =1 для формул левых и правых прямоугольников;

k = 2 для формул центральных прямоугольников и трапеций; k = 4 для формулы Симпсона.

8.6.Адаптивные процедуры численного интегрирования.

Пользователь такой процедуры задает отрезок [a,b], правило вычисления подынтегральной функции и требуемую точность ε > 0.

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

n	b
I | Iihi f (x)dx | . Здесь Iihi приближенное значение
i 1	a
интеграла Ii на отрезке		[xi-1,xi], вычисляемое по выбранной
квадратурной формуле.

Программа подбирает длину очередного отрезка так, чтобы

				n
	выполнялось неравенство Ii . Это условие будет выполнено,
				i 1
	если погрешность на очередном отрезке i				hi	. Действительно, в
					b a
	n	n	hi
	этом случае i			.
			b a
	i 1	i 1
				61

Алгоритм

1.Начальная установка ε, a1=a, b1=b, h=b1-a1, s=0.

2.Используя выбранную составную квадратурную формулу,

рассчитывают значения интеграла Ih и Ih/2, оценивают погрешность εi

по правилу Рунге.

	3. Если i	hi	, полагают h=h/2,b1=a1+h. Новый отрезок есть
		b a

[a1,b1]. Возвращаются к п.2.

4.Если условие п.3 не выполнено, то s=s+Ih/2 , a1=b1, b1=b, h=b- a1 (h длина очередного отрезка).

5.Если h > 0, возвращаются к п.2. Иначе s есть искомый результат.

8.7.Обусловленность квадратурных формул интерполяционного типа.

Значения подынтегральной функции всегда вычисляются с погрешностью. Предположим, что |f(x) – f*(x)| (f*) для всех x [a,b]. Тогда

b	b
| f (x)dx f * (x)dx | (b a) (f *) .
a	a
Для квадратурной формулы имеем
n	n		n
| Ai	f (xi ) Ai	f * (xi ) |	| Ai	|	(f*) .
i 1	i 1	i 1

Таким образом, абсолютное число обусловленности

	n
квадратурной формулы	ν = | Ai | .
	i 1

62

Все ИКФ точны для многочленов нулевой степени и поэтому b-

b	n
a= 1* dx Ai .
a	i 1
	n	n
	Если все веса ИКФ положительны, то ν = | Ai	| = Ai =b-a.
	i 1	i 1
	Если среди весов имеются отрицательные, то
n	n
| Ai	| > Ai =b-a. Например, для формул Ньютона -Котеса при n =
i 1	i 1

20, ν = 8.3(b-a), a при n = 30, ν = 560(b-a). Весовые коэффициенты формулы Гаусса все положительны при любом n.

8.8.Примеры использования квадратурных формул.

3 dx

Вычислить 1 2 x . Точное значение ln(5) = 1.6094379…

1) Используем составную формулу трапеций для n = 4 (число узлов равно 5).

63

x	x0	x1	x2	x3	x4	h = (3- (-1) )/4 = 1
	-1	0	1	2	3
f(x)	1	1/2	1/3	1/4	1/5

I

h

( f

f

)

i

i

i 1

2

f

0

f

n

I h

f1

... fn 1

2

2

n 1	f	0	f	n			1 1 / 5			1		1		1
h fi						1*										1.68333...
i 1			2					2		2		3		4

Погрешность реальная |R|ист = 0.074. Число верных значащих цифр 2.

Оценка погрешности

max | R		( f ) |		(b a)		h2 max | f ' ' ( ) |		(b a)			h2M		.
												2
	сост.тр.			12				12
			2					1							2
	f ''(x)				, M2		| f ''( 1) | 2, |R|			* 4 *12 * 2						.
			(2 x)3						12						3

Гарантирована только одна верная значащая цифра.

Рассчитаем оценку h для заданной погрешности δ.

					Если δ = 10-2, то h 0.122, а n	b a
	h		12	.				4	32.78.
			(b a)M 2			h		0.122

При n = 33, I = 1.61011… Вывод: при уменьшении h оценка погрешности становится более точной.

Оценим погрешность, используя правило Рунге. Для этого

64

сформируем таблицу с шагом h/2=1/2/

x	x0	x1	x2	x3	x4	x5	x6	x7	x8
	-1	-1/2	0	1/2	1	1.5	2	2.5	3
f(x)	1	2/3	½	2/5	1/3	2/7	¼	2/9	1/5

I=0.5((1+1/5)/2+2/3+1/2+1/3+2/7+1/4+2/9)=1.628968.

Согласно правилу Рунге погрешность

Ih/2=1/3(1.6833-1.628968)=0.0195303.

Относительно точного значения погрешность равна

1.628968-1.6094379=0.01446541

2) Используем составную формулу Симпсона.

n = 4 четное. Воспользуемся формулой без половинных узлов.

m =2.

			h								m		m 1				1	1	1	4		1			1						1
	I					( f0		fn		4 f2i 1 2 f2i )																		2						1.62222...
				3															5												3
											i 1		i 1				3				2 4
	max | R										( f ) |	(b a)	h4 max |			f (4) ( ) |					(b a)			h4M						.
																													4
						сост.симп.					180									180
	f (4) (x)			24				, M		24.
									4
				x)5
	(2
	|R|				24			* 4 * (1)4 0.5333...						8	. Оценка превышает
				180										15

соответствующую для формулы трапеций, потому что h = 1.

65

При заданном δ = 10-2

h		180	=		180 *10 2		= 0.37, n	4	10.8 .
		(b a)M 4				4 * 24		0.37

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

Iн.к (b a) 1 n Hk fk . Здесь n степень интерполяционного

N k 0

полинома, заменяющего подынтегральную функцию. Элементы Hk , N рассчитываются в зависимости от n по формулам, полученным при анализе ИКФ. Значения этих элементов сведены в таблицу.

m=n-1 H0 H1			H2 H3 H4 H5 H6 N Rm(f)
1	1	1	2		(b a)	3
							f (2) ( )
					12

2	1	4	1												6			(b a)	5
																						f (4) ( )
																	2880
3	1	3	3	1											8			(b a)	5
																						f (4) ( )
																	6480
4	7	32	12	32		7									90			(b a)		7
																							f (6) ( )
																	1935360
5	19	75	50	50		75			19						288			11(b a)					7
																								f (6) ( )
																	87800000
6	41	216	27	272		27			216			41			840
			3	dx
	Вычислим					, используя 7-ординатную формулу
				2 x
			1
	Ньютона – Котеса. h				b a					4			2	. Вычисления сводим в таблицу.
									6			3
						6
												66

k	xk	fk	Hk	Hkfk
0	-1	1	41	41
1	-1/3	3/5	216	129.6
2	1/3	3/7	27	11.571
3	1	1/3	272	90.666
4	5/3	3/11	27	7.3636
5	7/3	3/13	216	49.846
6	3	1/5	41	8.2
				338.2466

			1	n		1
	Iн.к	(b a)		Hk	fk 4		338.2466 1.61704216.
			N k 0		840
	Погрешность относительно точного значения I=1.27*10-3.
		1
	Вычислим ( 25x4				45x2 8)dx 4 с помощью 5- ординатной
		1

формулы Ньютона – Котеса. Эта формула должна быть точна, если подынтегральная функция есть полином степени не выше 4.

Проверим это. h=(b-a)/n = 2/4=1/2.Вычисления сводим в таблицу как в предыдущем примере.

k	xk	fk	Hk	Hkfk
0	-1	12	7	84
1	-1/2	1.6875	32	54
2	0	-8	12	-96
3	1/2	1.6875	32	54
4	1	12	7	84
			Σ=	180

67

I=2*180/90 = 4.

8.10.Применение формул Гаусса.

n

ИКФ Гаусса имеет вид IG Ak f (xk ) , где n степень

k 0

заменяющего интерполяционного полинома. Значения Ak и xk в

зависимости от n сведены в таблицу.

		0	1								2	3
Узлы	n
и веса
t0		0.0000				1				0.577350	-√0.6	-0.861136
A0		2.0000				3					5/9	0.3478448
			1.0000
t1			1					0.577350			0.0000	-0.3399810
A1					3						8/9	0.6521451543
			1.0000
t2											√0.6	0.3399810
A2											5/9	0.6521451541
t3												0.861136
A3												0.3478448

Осуществляя замену переменных, переходят к стандартному отрезку [-1,1]:

b		b a 1		a b			b a
	f (x)dx			f				t dt
		2			2		2
a			1
							68

	3		1	1		1
	В нашем варианте			dx 2			dt .
		2	x		2	(1 2t)
	1			1

Вычислим предыдущий интеграл с помощью трехординатной формулы Гаусса.

	I 2		5			1				8			1		5			1				1.6026935...
													(1 2*0)
			9 2			(1 2 0.6)				9 2					9 2			(1 2 0.6
	Погрешность относительно точного значения																				I=6.7*10-3.

Для четырехординатной формулы получаем I=1.60843…, для пятиординатной I=1.609288…

8.11. Замечания о точности квадратурных формул.

Погрешность для составных квадратурных формул R=O(hk), h=(b-a)/n. Чем больше k и чем меньше h, тем более точен результат.

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

1

Пример: f(x)=-8+45x2-25x4; I= f (x)dx 4 .

1

При h=1

a) составная формула трапеций Iсост.тр=1/2f(-1)+f(0)+1/2f(1)=6- 8+6=4;

б) составная формула Симпсона Iсост.симп.=1/3(f(-1)+4f(0)+f(1))= =1/3(12-32+12)=8/3. Однако, если использовать формулу трапеций при шаге h/2=1/2, получим

69

Iтр.сост.=1/2f(-1)+f(-1/2)+f(0)+f(1/2)+1/2f(1)= =6+1.6875-8+1.6875+6=7.375. Таким образом, полученное ранее точное значение есть результат случайного совпадения. Дальнейшее уменьшение шага будет уменьшать погрешность в обоих случаях,

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

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

x1	x2		x3			x4		x5
						0
	0		0					0

Использование адаптивных процедур в большинстве случаев

позволяет избежать неадекватных решений.

9.Дополнительные элементы линейной алгебры и теории матриц

70