Metod 1

МОСКОВСКИЙ АВИАЦИОННЫЙ ИНСТИТУТ (национальный исследовательский университет)

Битюков Ю.И.

МЕТОДЫ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ Методические указания к выполнению лабораторных работ в среде MathCad

для студентов, обучающихся по специальности «Cистемы управления летательными аппаратами»

Москва Издательство МАИ

2012

1

Методы численного интегрирования: Методические указания к выполнению лабораторных работ в среде MathCad для студентов, обучающихся по специальности «Cистемы управления летательными аппаратами» / Битюков Ю.И. – М.: Издательство МАИ, 2012.

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

Предназначено для студентов факультета «Системы управления информатики и электроэнергетики».

Рецензенты:

доц. каф. Молекулярной физики Физического факультета МГУ, к.ф.м.н. Иванов И.Э.

Кафедра «Аэродинамика, конструкция и прочность Летательных аппаратов» МГТУГА

© Московский авиационный институт, 2012

2

ПРЕДИСЛОВИЕ

Настоящее пособие содержит общие сведения по методам приближенного нахождения определенных интегралов: формула прямоугольников, трапеций и парабол. Кроме того, приведены примеры выполнения лабораторной работы по указанным методам с использованием среды MathCad.

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

MathCad.

3

ГЛАВА 1. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ

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

1.1. Многочлен Лагранжа и Ньютона

Пусть для дискретных значений аргумента x0,x1,...,xn известны значения

функции y0 f x0 , y1 f x1 ,...,yn f xn .

Построим многочлен Ln 1 x

степени не выше n, удовлетворяющий условию Ln 1 xi yi, i 0,1,...,n. Для начала найдем многочлены i x степени не выше n такие, что

	1, i j;
i xj		при i, j 1,...,n.
Поскольку i xj 0,	0, i	j
	i j, то	i x	делится на x xi	при i j. Таким

образом, нам известны n делителей многочлена степени n, поэтому

i x const x xj .

j i

Из условия i xi 1 получаем

	i x	x xj		.
		xi	xj
	i j
	Как легко видеть, многочлен		x xj
	n
	Ln 1 x yi				,
			xi xj
	i 0	i j

называемый интерполяционным многочленом Лагранжа, удовлетворяет условиям Ln 1 xi yi, i 0,1,...,n.

f x	В предположении непрерывности			f n 1 x				оценим разность между
	и построенным интерполяционным многочленом Ln 1 x .
Теорема 1.1. Если функция f Cn 1 a;b ,					то для всякого x a;b найдется
точка x a;b такая, что			f n 1
	f x L	x			n 1	x		a;b .
							,
			n 1 !
	n 1

Доказательство. Положим

t f t Ln 1 t K n 1 t

4

где n 1 t t x0 t x1 ... t xn , а K выберем из условия x 0, где x точка, в которой оценивается погрешность. Из уравнения x 0 получаем

K	f x Ln 1 x	.
	n 1 x
При таком выборе К функция x обращается в нуль в			n 2	точке
x0,...,xn,x. На основании теоремы Ролля, ее производная x			обращается в
нуль по крайней мере в n 1 точках. Применяя теорему Ролля к				x ,

получаем, что ее производная x обращается в нуль, по крайней мере, в n

точках. Продолжая эти рассуждения дальше, получаем, что n 1 x

обращается в нуль, по крайней мере, в одной точке a;b . Поскольку

n 1 t f n 1

t K n 1 !,

из условия n 1

0 будем иметь

K

f n 1

.

n 1 !

Следовательно, соотношение x 0 можно переписать в виде

f x L

x

f n 1

n 1

x

, a;b ,

n 1 !

n 1

который и дает представление остаточного члена. Теорема доказана.

Пример 1.1. По заданной таблице

x

0

1

2

построить многочлен Лагранжа.

y

0

2

1

Решение. Имеем

x0 0, x1 1, x2 2, f0 0, f1 1, f2 2.

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

x xj

2

2

3

2

7x

L3 x fi

x

.

xi

xj

2

i 0

j 0

2

j i

Понятие разделенной разности является обобщением понятия

f xi .

производной. Пусть в точках xi, i 1,2,... заданы значения функции

Определение. Разделенные разности нулевого порядка f xi совпадает со значением f xi ; разности первого порядка определяются равенством

f

f xj f xi

,

xi

,xj

xj

xi

разности k го порядка

f x1;...;xk

f

;...;xk

f

x2

1

.

x1

;...;xk 1

xk 1

x1

5

Лемма 1.1. Справедливо равенство

k		f xj
f x1,...,xk			.	(1.2)
	n
j 1 xj xi
	i 1
	i j
Доказательство. Доказательство будем проводить по индукции.				При k 1

это равенство превращается в равенство f x1 f x1 . Пусть равенство (1.2)

доказано при k l. Тогда,

f

f

x2;...;xl 1

f x1;...;xl

x1;...;xl 1

xl 1

x1

xj

f xj

1

l 1

f

l

xj

xi

xj xi

xl 1 x1 j 2

j 1

i j

i j

2 i l 1

f xj

1 i l

Если j 1, l 1, то коэффициент при

в правой части есть

1

1

1

xj

xi

xj xi

xl 1 x1

i j

i j

2 i l 1

1 i l

xj

x1 xj

xl 1

1

xl 1 x1

xj

xi

xj

xi

i j

i j

1 i l 1

1 i l 1

значение f xj входит

т.е. имеет требуемый вид;

для

j 1

или

j l 1

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

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

f x Ln x

f x

f xi x xj

n

i 1

j i

xi xj

n

f x

n

f xi

x xi

.

n

xi x xi xj

i 1

x xi

i 1

j i

i 1

Сравнивая с (1.2) можно написать

f x Ln x f x;x1;...;xn n x ,

(1.3)

6

n

где n x x xi .

	i 1
Пусть	Lm x	-	интерполяционный			многочлен Лагранжа			с	узлами
интерполяции x1,...,xm . Интерполяционный							многочлен	Ln x		можно
представить в виде
	Ln x L1 x L2 x L1 x ... Ln x Ln 1 x .
Разность	Lm x Lm 1 x			является		многочленом		степени		m 1,
обращающийся		в	нуль	в	точках		x1,...,xm 1,		поскольку

Lm 1 xj Lm xj f xj при 1 j m 1. Следовательно,

Lm x Lm 1 x Am 1 m 1 x .

Полагая x xm (и учитывая, что Lm xm f xm ) получим

f xm Lm 1 xm Lm xm Lm 1 xm Am 1 m 1 xm ,

С другой стороны, полагая в (1.3) n m 1, x xm , имеем

	f xm Lm 1				xm f		;x1;...;xm 1		m 1	xm .
						xm
Таким образом,	Am 1	f				поэтому
			xm	;x1;...;xm 1

Lm x Lm 1 x f x1;...;xm m 1 x .

Подставляя все эти выражения в Ln x , получим

Ln x f x1 f x1;x2 x x1 ...

... f x1;...

;xn x x1 ...

x xn 1 .

Интерполяционный многочлен, записанный в таком виде, называется интерполяционным многочленом Ньютона с разделенными разностями.

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

Пусть требуется найти определенный интеграл

		b
	I f x dx,			(1.4)
		a	a,b . Выразить интеграл через
где функция	f непрерывна на	отрезке
элементарные	функции удается	редко.	Поэтому обычно	функцию f

заменяют на такую функцию x , чтобы интеграл от нее легко вычислялся

в элементарных функциях. Чаще всего функцию f x	заменяют некоторым
многочленом
n
f x f xi i x r x ,	(1.5)
i 0

7

где r x - остаточный член аппроксимации. Подставляя (1.5) в (1.4), получим формулу численного интегрирования (квадратурную формулу)

n	b	b
I f xi ci R, ci	i x dx,	R r x dx,
i 0	a	a

где величины xi называют узлами, ci - весами, а R - погрешностью или остаточным членом формулы.

Рассмотрим более подробно применение для этой цели

интерполяционного многочлена Лагранжа. Заменяя функцию

f x

полиномом Ln 1 x , получим равенство

b

b

f

x dx Ln 1 x dx Rn f ,

(1.6)

Rn f

a

a

где

-

ошибка

квадратурной формулы (1.6). Отсюда,

используя

выражение для многочлена Лагранжа, получим приближенную формулу:

b

n

f x dx Ak yk

(1.7)

x xj

a

k 0

b

где

Ak

dx, k 0,1,...,n. Выберем

равномерное

разбиение

xk xj

a j k

xi n

,

xi a ih, h

b a

, i 0,1,...,n

отрезка a,b и

положим.

i 0

n

x x0

Выведем явные выражения для коэффициентов Ai . Положим q

тогда

x xi x x0 ih q i h, xj xi

j i h.

h

Следовательно, многочлен Лагранжа можно представить в виде

	n		x x0 ... x xi 1		x xi 1 ... x xn
	Ln 1 x							yi
			xi x0 ... xi xi 1		xi xi 1 ... xi xn
	i 0
	n	q q 1 ... q i 1 x q 1 ... q n
								yi
				i i 1 ...1	1 ... i n
	i 0
	n			q q 1 ... q n
	1 n i					yi
	i 0			i! n i ! q i

Поэтому, для коэффициентов Ak получаем следующее выражение:

Ak	b a		1 n k n				q q 1 ... q n		dq.
			k! n k !n 0
							q k
Введем следующее обозначение
		1 n k		n q q 1 ... q n
	Hk			0				dq.
		k! n k !n					q k

8

Эти постоянные называются коэффициентами Котеса. Формула (1.7) принимает вид:

b	n
f x dx b a Hi yi .		(1.8)
a	i 0

Рассмотрим некоторые частные случаи полученной формулы.

Формула трапеций. Пусть n 1. Тогда

	1	q q 1		1		1	1
	H0		dq		,	H1 qdq		.
				2			2
	0	q				0

Отсюда

b

h

f x dx 2 y0 y1 .

a

При получении этой формулы мы заменили подынтегральную функцию многочленом Лагранжа первой степени на отрезке [a;b]. Это соответствует замене кривой на секущую. Искомый интеграл, равный площади криволинейной фигуры, заменяется на площадь трапеции (рис.

1.1).

Предположим, что f C2 a,b и Рис. 1.1

найдем погрешность полученной формулы трапеций. Будем рассматривать остаточный член R h как функцию от h. Имеем

	R h	x0 h	f x dx	h	f x0 f x0 h
		x0	2

Отсюда, дифференцируя эту формулу по h последовательно два раза, получим:

R h f x0

h

f x0 f

x0

h

h f x0

h

2

2

f x0 h f x0

h f x0 h

.

2

2

R h

hf x0

h

.

2

Учитывая, что R 0 R 0 0, интегрируя по h и используя теорему о среднем, получим

			h						1 h
			0 R		t dt					tf x0 t dt
R	h R								2
			0							0
		1	h			h	2
			f 1 tdt					f 1 , 1 x0,x0 h .
	2		0		4

9

h

1

h

R h R 0 R t dt

t2 f 1 dt

0

4

0

1

h

h

3

f t2dt

f , x0,x0 h .

4

0

12

Рассмотрим теперь общую формулу трапеций. Промежуток

интегрирования a,b разделим на n равных частей

x0

,x1

, x1,x2

,..., xn 1,xn

b a

и к

каждому из них применим формулу

трапеций.

Полагая

h

и

обозначая yi

f xi , i 0,...,n, получим

n

b

h

h

h

y

y

0

n

f x dx

y0 y1

y1

y2 ...

yn 1

yn h

y1 ... yn 1

2

2

2

2

2

a

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

n

n xi

h

3

n

3

n

f i

R

f x dx

h

yi 1

yi

f i

h

i 1

.

n

i 1

2

12 i 1

12

xi 1

n

f i

Среднее арифметическое

i 1

заключено

между

наибольшим и

n

наименьшим значениями функции f x . Так как эта функция непрерывна,

то она принимает все свои промежуточные значения. Следовательно, существует число a,b такое, что f . Таким образом,

R h3n f . 12

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

формулы (1.8). Пусть n 2. Тогда
1 1				2							1		11						2	2
H0				0 q 1 q 2 dq								,	H1						0 q q 2 dq		,
	2		2								6			21						3
					1 1			2					1
				H2				0 q q					1 dq					.
					2		2								6
Так как x2 x0 2h, то
				b					h
				f x dx						y0 4y1 y2 .
				a				3

10