книги / Численные методы решения задач строительства. Ч. 1

y
y = f(x)
		Si
			x
x0 = a	xi	xi+1	xn = b
Рис. 4.1. Геометрическая интерпретация
численного интегрирования

Задачу вычисления определенного интеграла (4.1) заменяем задачей вычисления площади этой криволинейной трапеции. Однако задача нахождения площади криволинейной не является простой.

Отсюда идея численного интегрирования [3, 6] будет заключаться в замене криволинейной трапеции фигурой, площадь которой вычисляется достаточно просто.

Для этого отрезок интегрирования [a, b] разобьем на n равных элементарных отрезков [xi, xi+1] (i = 0, 1, 2, ….., n – 1), с шагом h = (b – a)/n. При этом криволинейная тра-

пеция разобьется на n элементарных криволинейных тра-

пеций с основаниями, равными h (см. рис. 4.1).

Каждая элементарная криволинейная трапеция заменяется фигурой, площадь которой вычисляется довольно просто. Обозначим эту площадь Si. Сумма всех этих площадей называется интегральной суммой и вычисляется по формуле

n
n Si .	(4.3)

i 1

Тогда приближенная формула вычисления определенного интеграла (4.1) имеет вид

81

J b f x dx n .	(4.4)
a

Точность вычисления по формуле (4.4) зависит от шага h, т.е. от числа разбиений n. С увеличением n интеграль-

ная сумма n приближается к точному значению интеграла

J lim n .

(4.5)

n

Это хорошо проиллюстрировано на рис. 4.2.

σn

J Точное значение интеграла

n

Рис. 4.2. Зависимость точности вычисления интеграла от числа разбиений

В математике доказывается теорема: если функция y = f(x) непрерывна на [a, b], то предел интегральной суммы σn существует и не зависит от способа разбиения отрезка [a, b] на элементарные отрезки.

Формулу (4.4) можно использовать, если известна степень точности такого приближения. Существуют различные формулы для оценки погрешности выражения (4.4), но, как правило, они достаточно сложны. Будем проводить оценку точности приближения(4.4) методом половинного шага.

Алгоритм метода половинного шага

Циклически повторяется следующая последовательность действий:

82

1.Разбиваем отрезок интегрирования [a, b] на n равных отрезков с шагом h = (b – a)/n.

2.Строим интегральную сумму σп по формуле (4.3).

3.Повторяем эти вычисления (п. 1, 2) для шага h/2, т.е. для 2n. Строим интегральную сумму σ2п.

4.Если два соседних приближения (две итерации) близ-

ки, т.е.

| n 2n | ,

(4.6)

то за приближенное значение интеграла (4.1) с точностью принимаем интегральную сумму σ2п:

J 2n

(4.7)

5. Если условие (4.6) не выполняется, то надо вернуться на пункт 3, т.е. еще раз уменьшить шаг вдвое. Итерационный процесс продолжим до тех пор, пока условие (4.6) не будет выполнено.

Рассмотрим несколько методов численного интегрирования.

4.1. Квадратурные формулы прямоугольников

Отрезок интегрирования [a, b] разбиваем на n равных отрезков и получаем множество (n+1) равноудаленных уз-

лов (равномерную сетку Ωn):

Ωn = {x0 = a, xi+1 = xi +h, i = (0, 1, …, n – 1), xn = b, h = (b – a)/n},

где h – шаг разбивки.

Введем обозначение уi = f(xi), i = (0, 1, …, n). Площадь каждой элементарной криволинейной трапе-

ции заменим площадью прямоугольника с основанием h

ивысотой f( i), где i [xi, xi+1], i = 0, 1, 2, …, n–1 (рис. 4.3).

Взависимости от выбора i существует несколько

формул прямоугольников.

83

Рис. 4.3. Схема метода прямоугольников

Формула «левых» (входящих) прямоугольников, когда

i = xi,

n 1
J h f (xi ).	(4.8)

i 0

Формула «правых» (выходящих) прямоугольников, когда i = xi+1,

n
J h f (xi ).	(4.9)

i 1

Формула «средних» прямоугольников, когда i = xi + h/2,

J h f xi h			.	(4.10)
n 1
i 0		2

Пример 4.1. С помощью формул левых и правых пря- мо-угольников вычислить интеграл

	J 9	dx		, полагаяn 4.
		x	2
	1

Решение. Зная пределы интегрирования а =1, b = 9, находим шаг h = (b – a)/n = 2.

Тогда точками разбиения (узлами) служат х0 = 1, х1 = 3, х2 = 5, х3 = 7, х4 = 9, а значения подынтегральной функции

	f x	1	в этих узлах равны соответственно:
		x 2
	84

x	1	3	5	7	9
y = f(x)	1	1	1	1		1
	3	5	7	9	11

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

		4	h y	0	y		y	2	y		2 1		1	1		1	1,6024.
						1				3			5	7		9
												3
Вычислим интеграл по формуле правых прямоуголь-
ников:
			h y		y		y		y		2	1	1	1		1				1,1053,
	4					2				4
			1				3						7	9		11
												5

что достаточно близко совпадает со значением интеграла, вычисленного по формуле Ньютона–Лейбница:

	9	dx	9
J	1		ln | x 2 |	1,299.
		x 2	1

4.2. Квадратурная формула трапеций

Площадь каждой элементарной криволинейной трапеции заменим площадью линейной трапеции с основа-

ниями yi = f(xi) и yi + 1 = f(xi + 1) и высотой h (i = 0, 1, 2, ..., n–1) (рис. 4.4).

Рис. 4.4. Схема метода трапеций

85

Формула трапеций имеет вид

			n 1		y		y
J h					i		i 1		(4.11)
							2
			i 0
или
	y		0	y		n	n 1
J h							yi .		(4.12)
				2
							i 1
Пример 4.2. С помощью формулы трапеций вычис-
лить интеграл примера 4.1.
J 9	dx		, полагаяn 4.
	x 2
1

Решение. Здесь f x x 1 2 . Шаг h = 2. Точки раз-

биения (узлы): х0 = 1, х1 = 3, х2 = 5, х3 = 7, х4 = 9. Тогда по формуле (4.12) получим

			1		1
								1		1		1
	4	2		3	11									1,3322,
								5		7		9
					2

что практически совпадает со значением интеграла, вычисленного по формуле Ньютона-Лейбница:

	9	dx	9
J	1		ln | x 2 |	1,299.
		x 2	1

4.3. Квадратурная формула Симпсона

Точность приближенного интегрирования заметно возрастает, если подынтегральную функцию y = f(x) на отрезке [xi–1, xi+1] длиной 2h заменить квадратичной

86

функцией,	проходящей		через три точки A(xi–1, yi–1),
B(xi, yi), C(xi+1, yi+1) (рис. 4.5).
y
		Погрешность		2
				y = ax + bx + c
		y = f(x)	B
				C
		A

xi–1

xi

xi+1

x

Рис. 4.5. Схема формулы Симпсона

Площадь элементарной криволинейной трапеции с основанием [xi–1, xi+1] по формуле Симпсона имеет вид

Si h3 yi 1 4yi yi 1 .

Для определенного интеграла (4.1) отрезок [a, b] разбивается на четное число отрезков n = 2m c шагом h = = (b – a)/2m. Формула Симпсона в общем случае может быть записана так:

. m h y0	4 y2k 1	2 y2k y2m		, .	(4.13)
	m	m
3	k 1	k 1

где yi = f(xi), i = 0, 1, 2, …, n.

Формула Симпсона обладает повышенной точностью.

Алгоритм для приближенного вычисления интеграла (4.1) по формуле Симпсона:

1.n = 2m;

2.h = (b – a)/2m;

3.x0 = a, xi+1 = xi +h (i = 0, 1, …, n – 2), xn = b;

87

4.уi = f(xi), (i = 0, 1, …, n);

5.M0 = y0 + yn = f(a) + f(b);

m

6. M1 y2k 1 ;

k1 m

7.M2 y2k ;

k1

8.J h3 M0 4M1 2M2 .

Для вычисления интеграла с заданной степенью точности ε надо использовать метод половинного шага, изложенный выше.

Пример 4.3. Вычислить по формуле Симпсона интеграл:

				J 9				dx			, полагая n 4, шагh 2.
							x 2
				1
	Решение. По формуле (4.13) имеем
				9		dx		h y 4y 2y							4y y					.
														2					4
				1 x 2					3		0	1			3
	Подставляя							значения				подынтегральной функции
	f x	1		вузлах х0 = 1, х1 = 3, х2 = 5, х3 = 7, х4 = 9, получим
		x 2
		9	dx				2 1				1	1	1			1
		1					3			4 5 2			4 9					1,3029,
			x 2						3			7				11

что практически совпадает со значением интеграла, вычисленного ранее по формуле Ньютона–Лейбница.

Формула Симпсона, в частности, используется в строительной механике стержневых систем при вычислении интегралов Мора, где подынтегральной функцией является

88

произведение эпюр изгибающих моментов (М1Мр) или других внутренних усилий. Результат получается точным, если обе перемножаемые эпюры прямолинейны (их произве-

дение – квадратная парабола) или одна из эпюр – параболическая, а другая линейная (произведение – кубическая парабола). Формула Симпсона применима и в тех случаях, когда стержень имеет переменное сечение или криволинейное очертание.

Пример 4.4. Перемножить две эпюры Мр и М1 (рис. 4.6), используя формулу Симпсона и полагая EJ = Const.

Рис. 4.6. Перемножаемые эпюры

	l	M p M1		l	ap a1		bpb1 .
			dx			4СpС1
		EJ		EJ
	0

Перемножать эпюры по формуле Симпсона следует на участках, где эпюры меняются плавно, без скачков и переломов. При наличии же таковых (в местах приложения сосредоточенных моментов или сил) перемножение надо производить на каждом отдельном участке, где функции меняются плавно.

89

4.4. Реализация методов численного интегрирования средствами приложения MS Excel

Вычислим определенный интеграл
J 2 x2dx	(4.14)
1

методом прямоугольников (входящих) и методом трапеций.

Последовательность действий:

На отрезке x [a, b] построим разностную сетку

n{x0 = a, xi = xi–1 + h, i = 1, 2, ..., n – 1, xn = b, h = (b – a)/n}

и создадим таблицу по образцу рис. 4.7.

Рис. 4.7.Схема численного интегрирования

1.В ячейки В1, В2 и D1 введем значения нижнего

иверхнего пределов интегрирования а, b и количество разбивок n.

2.В ячейке В3 вычислим шаг разбивки h: В3 = (В2 –

– В1)/D1.

90