Численные методы (лекции)2013

Это условие выполнено при 0 < τ < 2.

Теорема 6.8. Итерационный метод (1) сходится тогда и только тогда, когда все собственные значения матрицы S = E−τ B−1 A по модулю меньше единицы.

7. Лекция 7

Рассмотрим способы численного вычисления определенных интегралов. Пусть

Z b

I = f (x)dx.

a

Будем искать приближенное значение в виде

n

X

In = In(f ) = cif (xi),

i=0

где a = x0 < x1 < · · · < xn = b. Такая сумма называется квадратурной суммой. Точки xi называются узлами квадратуры, а числа ci коэффициентами квадратуры. Разность

Rn = I − In

называется погрешностью квадратуры.

Утверждение 7.1. Квадратура In(f ) является линейным функционалом, т.е. In(f + g) = In(f ) + In(g), In(Cf ) = CIn(f ).

Доказательство. Действительно,

n	n	n
X	X	X
In(f +g) = ci(f (xi)+g(xi)) =	cif (xi)+ cig(xi) = In(f )+In(g).
i=0	i=0	i=0
Аналогично,
n		n
X		X
In(Cf ) = ciCf (xi) = C		cif (xi) = CIn(f ).
i=0		i=0

Предположим, что узлы квадратуры расположены равномерно, т.е. xi = a + ih, где nh = b − a. Тогда

b	n	xi	f (x)dx.
I = Za	f (x)dx = i=1	Zxi−1
	X

20

Для построения квадратуры на отрезке [a; b] достаточно построить квадратуру для интеграла

	xi
	Zxi−1	f (x)dx.
7.1. Формула прямоугольников. Пусть xi−1/2			= xi−1	+ h .
xi				2
Заменим интеграл Rxi−1 f (x)dx выражением f (xi−1/2)h. Формула
Z xi	f (x)dx ≈ f (xi−1/2)h

xi−1

называется формулой прямоугольников. Найдем погрешность формулы прямоугольников. Применим формулу Тейлора

f (x) = f (xi−1/2) + f ′(xi−1/2)(x − xi−1/2) + f ′′(ζ) (x − xi−1/2)2. 2

Тогда

xi

xi

ri = Zxi−1 f (x)dx − f (xi−1/2)h = Zxi−1 (f (x) − f (xi−1/2))dx =

xi

f ′′(ζ)

Zxi−1 (f (xi−1/2) + f ′(xi−1/2)(x − xi−1/2) +

(x − xi−1/2)2 − f (xi−1/2))dx =

2

xi

f ′′(ζ)

xi

f ′(xi−1/2) Zxi−1 (x − xi−1/2)dx +

Zxi−1 (x − xi−1/2)2dx =

2

f ′′(ζ)

(x − x

)3

(x

i−1

− x

)3

=

f ′′(ζ)h3

i

i−1/2

−

i−1/2

.

2

3

3

24

i − xi−1/2

= xi−1/2

i

=

2

xi

( )

что

Ri,2 = maxx [xi−1,xi]

и

Здесь мы

воспользовались

тем,

xi−1 (x − xi−1/2 )dx

=

0

x

− x

h . Пусть M

|f ′′

x | и

M2 = maxx [a,b] |f ′′(x)|. Получаем

|ri| ≤

Mi,2h3

≤

M2h3

.

24

24

Рассмотрим квадратуру

n

X

In =

f (xi−1/2)h.

i=1

21

Оценим ее погрешность

	n			n	ri| ≤
|Rn| = |I − In| = | Za b f (x)dx − i=1 f (xi−1/2)h| = |				i=1
n	X			X
	= M2h2(b − a) .
|ri| ≤ n M2h3
X
24	24
i=1

Здесь мы воспользовались тем, что nh = b − a.

7.2. Формула трапеций. Формула трапеций получается, если заменить f (x) интерполяционным многочленом Лагранжа пер-

вой степени, т.е. мы заменяем f (x) на

1

((x − xi−1)f (xi) − (x − xi)f (xi−1)) .

h

Тогда

xi

f (x

+ f (x )

Zxi−1 f (x)dx ≈

i−1)

i

h.

2

Оценим погрешность метода трапеций. Согласно ??,

1

f ′′(ζ)

f (x) −

((x − xi−1)f (xi) − (x − xi)f (xi−1)) =

(x −xi−1 )(x −xi).

h

2

Получаем

xi

f (x

+ f (x )

ri = Zxi−1 f (x)dx −

i−1)

i

h =

2

xi

1

Zxi−1

(f (x) −

((x − xi−1)f (xi) − (x − xi)f (xi−1)))dx =

h

xi f ′′(ζ)

f ′′(ζ)

h/2

h

h

Zxi−1

(x − xi−1)(x − xi)dx =

Z−h/2(t −

)(t +

)dt =

2

2

2

2

2

−h/2 t2dt − −h/2

4 dt! =

2

( 12 −

4 ) = −

12 .

f ′′(ζ)

Z

h/2

Z

h/2

h2

f ′′(ζ) h3

h3

f ′′(ζ)h3

Здесь t = x − xi−1/2 = x − xi−1

+ h . Таким образом,

2

Mi,2h3

M2h3

|ri| ≤

≤

.

12

12

Рассмотрим квадратуру

In = n

f (xi−1) + f (xi) h = h( f (x0)

+ f (x1) + · · ·+ f (xn−1) + f (xn) ).

X

i=1

2

2

2

22

Оценим ее погрешность

	n				n
|Rn| = |I − In| = | Za b f (x)dx − i=1		f (xi−1)2+ f (xi) h| = | i=1 ri| ≤
X					X
n	=		M2h2(b − a) .
|ri| ≤ n M2h3
X
12		12
i=1

Здесь мы воспользовались тем, что nh = b − a.

7.3. Формула Симпсона. В формуле Симпсона мы заменяем f (x) на параболу, проходящую через точки

(xi−1, f (xi−1)), (xi−1/2, f (xi−1/2)), (xi, f (xi)) т.е. приближаем f (x)

интерполяционным многочленом Лагранжа второй степени. Имеем

		xi			xi
где		Zxi−1		f (x)dx ≈ Zxi−1		L2,i(x)dx,
L2,i(x) =	2	((x − xi)(x − xi−1/2)f (xi−1) − 2(x − xi)(x − xi−1)f (xi−1/2)
	2
	h
		+(x − xi−1)(x − xi−1/2 )f (xi)).
Интегрируя L2,i(x), получим
xi	2				xi
					f (xi−1) Zxi−1 (x − xi)(x − xi−1/2)dx−
Zxi−1 L2,i(x)dx =
			h2
		xi				xi

Z Z

2f (xi−1/2) (x − xi)(x − xi−1)dx + f (xi) (x − xi−1)(x − xi−1/2 )dx

xi−1 xi−1

2

= h2 f (xi−1)

Z h/2

f (xi) t(t +

−h/2

h/2

h

h/2

h

h

Z−h/2

(t −

)tdt − 2f (xi−1/2) Z−h/2

(t −

)(t +

)dt+

2

2

2

2 )dt!

=

h2

(f (xi−1) 12

− 2f (xi−1/2)( 12 −

4 ) + f (xi) 12 )

h

2

h3

h3

h3

h3

h

2h

h

=

f (xi−1) +

f (xi−1/2) +

f (xi).

6

3

6

Формула

xi

h

2h

h

Zxi−1 f (x)dx ≈

f (xi−1) +

f (xi−1/2) +

f (xi)

6

3

6

называется формулой Симпсона или формулой парабол.

Лемма 7.2. Формула Симпсона точна для любого многочлена третей степени.

23

Доказательство. Согласно утверждению 7.1 достаточно доказать лемму только для функций 1, x, x2, x3. Для функций 1, x, x2 утверждение очевидно, поскольку формула Симпсона приближает функцию многочленом второй степени. Докажем для функции x3. Имеем

b

b − a

2(

b − a

a + b

3

b − a

Za

x3dx ≈

a3 +

)

+

b3 =

6

3

2

6

b

− a

a3

a3

3

a2b

3

ab2

b3

b3

+

+

+

+

+

=

3

2

4

4

4

4

2

b − a

a3

3

a2b

ab2

b3

=

3

a2b

+

ab2

+

b3

)

3

+

+

3

+

3

(b − a)(a +

3

4

4

4

4

4

=

b4 − a4

.

4

С другой стороны,

Za

b

b4

a4

x3dx =

−

.

4

4

Отсюда следует утверждение леммы.

Пусть H3(x) многочлен Эрмита третей степени такой, что

H3(xi−1) = f (xi−1), H3(xi) = f (xi),

H3(xi−1/2) = f (xi−1/2), H3′ (xi−1/2) = f ′(xi−1/2 ).

Представим

f (x) = H3(x) + ri(x),

где r(x) погрешность интерполирования. Тогда

xi

h

2h

h

xi

Zxi−1

f (x)dx −

f (xi−1) +

f (xi−1/2) +

f (xi) = Zxi−1

ri(x)dx.

6

3

6

Мы знаем, что

ri(x) =

f (IV )(ζ)

(x − xi−1)(x − xi−1/2)2(x − xi).

24

Отсюда,

xi

xi

V )

Zxi−1 ri(x)dx = Zxi−1

f (I (ζ)

(x − xi−1)(x − xi−1/2)2(x − xi) =

24

f (IV )

ζ

)

h/2

h

h

f (IV ) ζ

)

h/2

h2

(

Z−h/2 t2(t −

)(t +

)dt =

(

Z−h/2(t4 −

t2)dt =

24

2

2

24

4

f (IV )(ζ)

h5

h5

= −

f (IV )(ζ)h5

−

.

24

80

48

2880

24

Аналогично, как и для методов прямоугольников и трапеций,

n	M4h4(b − a) ,
|Rn| ≤ |ri| ≤
X
i=1	2880

где M4 = maxx [a;b] |f (IV )(x)|.

В конце запишем составную формулу Симпсона

		Za	b	n	h		2h			h
			f (x) =	i=1 (					f (xi−1/2) +
					6 f (xi−1) +		3			6 f (xi)) =
	h			X
		(f (x0) + 4f (x1/2) + 2f (x1) + 4f (x3/2) + · · · + 2f (xn−1)
	6

+4f (xn−1/2) + f (xn)).

7.4. Метод экстраполяции Ричардсона. Приведем метод вывода формулы Симпсона, основанный на методе экстраполяции. Рассмотрим формулу трапеций

Ih =	n	f (xi−1) + f (xi)	h.
	2
	i=1
	X
Тогда, используя разложение в ряд Тейлора, получаем
Ih = I + ch2 + O(h4).				(2)

Здесь мы воспользовались тем, что коэффициенты при членах с нечетными степенями равны нулю. Рассмотрим шаг h/2. Тогда

n

i−1

2

i−1/2

+

i

2

i−1/2

)

2 .

Ih/2 = i=1

X

f (x

) + f (x

)

f (x ) + f (x

h

С другой стороны,

						h2
	Ih/2 = I + c									+ O(h4).	(3)
						4
Умножим (2) на 1	и вычтем из (3). Получим
4
	1					3			I + O(h4).
	Ih/2 −				Ih =
				4			4
Тогда
	4					1		Ih + O(h4).
	I =		Ih/2 −
		3				3

25

Отсюда

4

1

+ O(h4) =

Jh =

Ih/2

−

Ih

3

3

n

f (x

) + f (x

) f (x ) + f (x

)

f x

+ f x

h

i=1

i−1

i−1/2

−

3

+

i

3

i−1/2

(

i−1)6

( i)

X

n

i−1

) + 4f (x

i6

h.

= i=1

X

f (x

−1/2) + f (xi−1/2)

Это и есть формула Симпсона. Пусть Ih квадратура с шагом h,

Ih = I + a1hα1 + a2hα2 + · · · + amhαm + O(hαm+1 ),

где 0 < α1 < α2 < . . . < αm < αm+1 и ai не зависят от h. Рассмотрим последовательность шагов hk = qkh0, где q < 1. Обозначим Ih1 = Ih.

Тогда

1

α1

α2

αm

αm+1

),

Ihk = I + a1hk

+ a2hk

+ · · · + amhk

+ O(hk

1

α1

α2

αm

αm+1

Ihk−1

= I + a1hk−1

+ a2hk−1

+ · · · + amhk−1

+ O(hk−1

).

Исключим a1. Получим

1

− q

α1

1

α1

α2

αm

αm+1

).

Ihk

Ihk−1 = (1 − q

)I + b2hk−1 + · · · + bmhk−1 + O(hk−1

Отсюда,

I1

− qα1 I

1

′

α2

′

αm

αm+1

hk

hk−1

I =

+ b h

+ · · · + b h

+ O(h

k−1

).

1 − qα1

2

k−1

m

k−1

Здесь мы пользуемся тем, что O(hk) = O(hk−1). Обозначим

2

I

1

− qα1 I1

hk

hk−1

Ihk−1 =

.

1 − qα1

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

квадратура, приближающая с точностью

до O(hα2

hk−1

). Этот процесс повышения точности можно продолжить.

k−1

Получим

j+1

Ij

− qαj Ij

I

=

hk

hk−1

.

hk−1

1 − qαj

Такой метод повышения точности называется методом экстраполяции Ричардсона.

26

b−a

n

7.5. Метод Рунге. Оценку погрешности можно осуществлять методом Рунге. Пусть какая-то квадратура имеет на отрезке [xi−1; xi] порядок точности m т.е.

	Ii − Ii,h ≈ cihm,
xi	f (x)dx, Ii,h квадратураhнаmотрезке [xi−1; xi]. Тогда
где Ii = Rxi−1
	Ii − Ii,h/2 ≈ ci		.
		2

Отсюда,

Ii − Ii,h ≈ 2m(Ii − Ii,h/2).

Получаем

	Ii,h/2	− Ii,h
Ii − Ii,h/2 ≈			.
	2m	− 1

Таким образом, мы можем задавать точность интерполирования. Пусть ε > 0 такое, что

	|Ii − Ii,h/2| ≈	|Ii,h/2 − Ii,h|		< εh.
		2m − 1
	Рассмотрим квадратуру
		n	Ii,h = Ih.
	I = Za b f (x)dx ≈ i=1
	Тогда	X
	n		n
	|I − Ih/2| ≤ |Ii − Ii,h/2| ≈		εh = ε(b − a).
	i=1		i=1
	X		X

8.Лекция 8

8.1.Решение нелинейных уравнений. Пусть задана функ-

ция f (x). Требуется найти (приближенно) решение уравнения f (x) = 0. Обычно, для этого используют итерационные методы.

Предположим, что f (x) непрерывна и имеет корень на отрезке [a; b]. Рассмотрим разбиение отрезка на n равных частей

a = x0 < x1 < x2 · · · < xn = b, xi = a + ih, где h = шаг. Вычислим значение функции в точках x0, x1, . . . , xn. Предположим,

что существуют два значения f (xi) и f (xi+1) различные по знаку (если таких значений нет, то мы можем взять разбиение с большим n). Тогда на отрезке [xi; xi+1] функция f (x) имеет нечетное число корней (т.е. по крайней мере, один).

Далее, будем искать корень на [xi; xi+1] методом деления отрез-

ка пополам (методом бисекции). Пусть y0 = xi+xi+1 . Рассмотрим

2

27

f (y0). Предположим f (y0) =6 0 (иначе мы уже нашли корень). Поскольку f (xi) и f (xi+1) имеют различные знаки, то либо f (y0) и f (xi) имеют различные знаки, либо f (y0) и f (xi+1). Предположим, что f (y0) и f (xi+1). Тогда пусть y1 = xi+1. Рассмотрим отрезок [y0; y1]. На концах этого отрезка функция f (x) имеет различные значения. Тогда f (x) имеет корень на [y0; y1]. Возьмем y2 = y0+2 y1 . Аналогично, либо f (y2) и f (y0) имеют различные знаки, либо f (y1) и f (y0). Снова, возьмем середину отрезка, на котором функция f (x) имеет различные знаки и т.д. Таким образом, мы построили последовательность отрезков, длина которых стремится к нулю. Пусть z1, z2, . . . , zn, . . . последовательность середин этих отрезков. Тогда {zm} сходится к корню уравнения f (x) = 0. Действительно, так как длина m-го отрезка 2mh−1 , то |zm − zm+p| < 2hm . Следовательно, {zm} сходится по критерию Коши. В силу непрерывности f (x),

	lim f (zm) =	lim f (zm −			h	) =		lim f (zm +		h	).
					2m					2m
	m→∞	m→∞						m→∞
	С другой стороны, f (zm −		h	) и f (zm +				h	) имеют разные знаки.
			m					m
			2				2

Рассмотрим метод простой итерации. Уравнение f (x) = 0 заменяем на уравнение x = s(x). Выберем начальное значение x0. Далее, xn+1 = s(xn). Для сходимости итерационного метода важно правильно выбрать s(x). Мы можем взять

s(x) = x + τ (x)f (x),

где τ (x) непрерывна и не обращается в ноль на том отрезке, где ищется корень. Говорят, что итерационный метод сходится, если существует предел последовательности {xn}. Заметим, что если метод простой итерации сходится, то он сходится к корню уравнения f (x) = 0. Действительно, пусть xn → x . Тогда

xn+1 = s(xn) = xn + τ (xn)f (xn).

Переходя к пределу в обеих частях, получаем

x = x + τ (x )f (x ).

Отсюда,

τ (x )f (x ) = 0.

Поскольку τ (x ) =6 0, то f (x ) = 0.

Определение 8.1. Функция f (x) называется липшицнепрерывной с постоянной q на множестве Υ, если для любых точек x1, x2 Υ выполнено неравенство

|f (x1) − f (x2)| < q|x1 − x2|.

28

Теорема 8.2. Пусть s(x) липшиц-непрерывная функция с па-

раметром 0 < q < 1 на отрезке Ur (a) = {x | |x − a| ≤ r}, причем |s(a) −a| ≤ (1 −q)r. Тогда уравнение x = s(x) имеет единственное

решение x на отрезке Ur (a) и метод простой итерации сходится к x при любом начальном x0 Ur (a).

Доказательство. Докажем по индукции, что xi не выходят за пределы отрезка Ur (a). Пусть xi Ur (a). Тогда

|xi+1 − a| = |s(xi) − a| = |s(xi) − s(a) + s(a) − a| ≤ |s(xi) − s(a)|+ +|s(a) − a| ≤ q|xi − a| + (1 − q)r ≤ qr + (1 − q)r = r.

Оценим теперь

|xi+1 − xi| = |s(xi) − s(xi−1)| ≤ q|xi − xi−1| ≤ · · · ≤ qi|x1 − x0|.

Тогда

|xi+p − xi| = |xi+p − xi+p−1 + xi+p−1 − xi+p−1 + · · · + xi+1 − xi| ≤

|xi+p − xi+p−1| + |xi+p−1 − xi+p−1| + · · · + |xi+1 − xi| ≤ qi+p−1|x1 − x0| + qi+p−2|x1 − x0| + · · · + qi|x1 − x0| =

qi	1 − qp	|x1 − x0| < qi	1	|x1	− x0|.
	1 − q
			1 − q

Последнее выражение стремится к 0 при i → ∞. Тогда последовательность {xn} сходится. Пусть xn → x . Переходя к пределу в равенстве xn+1 = s(xn), получаем x = s(x ).

Следствие 8.3. Пусть s(x) функция, дифференцируемая на отрезке Ur (a), и |s′(x)| ≤ q < 1 на Ur (a). Более того, |s(a) − a| ≤ (1 − q)r. Тогда метод простой итерации сходится при любом начальном x0 Ur (a).

Доказательство. Пусть x, y Ur (a) и x > y. Тогда, по теореме Лагранжа, s(x) − s(y) = s′(ξ)(x − y), где ξ [y; x]. Отсюда,

|s(x) − s(y)| = |s′(ξ)(x − y)| ≤ q|x − y|.

Следовательно, s(x) липшиц-непрерывная функция с параметром 0 < q < 1 на отрезке Ur (a).

8.2. Метод Ньютона. Вернемся к задаче отыскания корня уравнения f (x) = 0. Пусть x0 начальное приближение. Разложим f (x) в ряд Тейлора до первого члена, т.е.

f (x) ≈ f (x0) + f ′(x0)(x − x0).

29