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

4.2. Интерполирование сплайнами. Пусть задано разбие-

ние a = x0 < x1 < · · · < xn−1 < xn = b отрезка [a; b]. Функция Sm(f, x) называется сплайном порядка m, если Sm(f, x) является

многочленом степени m на каждом из отрезков [xi−1; xi], т.е.

Sm(x) = Pim(x) = ai0 + · · · aimxm, x [xi−1; xi],

и удовлетворяющую условиям непрерывности производных до порядка m − 1 в точках x1, . . . xn−1:

Pi,m(k)(xi) = Pi(+1k) ,m(xi), k = 0, 1, . . . m − 1, i = 1, . . . , n − 1.

Всего в слайне (m + 1)n переменных. Из условий равенства m −1 производной в n − 1 точке, получаем m(n − 1) уравнений. Кроме того, из условия равенства сплайна функции f (x) в узлах интерполяции получаем еще n + 1 уравнение. Следовательно, мы имеем nm+ n переменные и nm−m+ n+ 1 уравнение. Таким образом, нам необходимо еще m − 1 условие. Их мы можем задать в граничных точках.

5.Лекция 5

5.1.Системы линейных уравнений. Метод Гаусса. Начнем с самой простой задачи. Дано линейное уравнение

AX = B,

где A матрица n × n с определителем, отличным от нуля, X = (x1, . . . , xn)T искомый вектор, B = (b1, . . . , bn)T заданный вектор. Будем решать систему уравнений методом Гаусса. Т.е.

a11 a12

a21 a22

. . . . . .

an1 an2

1 a′120 a′22

. . . . . .

0 a′n2

. . . . . . a1n

. . . . . . a2n

. . . . . . . . .

. . . . . . ann

. . . . . . a′1n

. . . . . . a′2n

. . . . . . . . .

. . . . . . a′nn

1

0

. . .

0

b2	→		a21
b1					1
.b.n.			an1
					. . .
b1′					1
b2′		→			0
.b.′ .				.	0
n					. .

0. . . . . . 0

1 . . . . . . 0

. . . . . . . . . . . .

0 . . . . . . 1

a′	. . . . . . a′
12	1n
a22	. . . . . . a2n
. . . . . . . . . . . .
an2	. . . . . . ann
a′	. . . . . . a′
12	1n

1. . . . . . a′2n

. . . . . . . . . . . .

0. . . . . . 1

b(1n)

b(2n)

. . .

b(nn)

b′1

b2 →

. . . bn

b′1

b′′

2 →

. . .

b(nn)

10

Заметим, что здесь мы предполагаем, что все a(iii−1) 6= 0.

5.2. Подсчет числа операций. Поскольку операции умножения и деления на ЭВМ требуют гораздо большего времени чем операции сложения и вычитания, ограничимся подсчетом числа умножений и делений.

Число операций при делении на a(iii−1) равно

(n − 1) + (n − 2) + · · · + 1 = n(n − 1) . 2

Число операций при вычислении a(ijk) j > i равно

(n − 1)2 + (n − 2)2 + · · · + 12 = n(n − 1)(2n − 1) . 6

Вычисления правых частей при делении на a(iii−1) требуют n операций. Вычисления правых частей при вычитании уравнений требу-

ют

(n − 1) + (n − 2) + · · · + 1 = n(n − 1) . 2

Для обратного хода метода Гаусса требуется

(n − 1) + (n − 2) + · · · + 1 = n(n − 1) . 2

Таким образом, общее число операций равно

n(n − 1)

+

n(n − 1)(2n − 1)

+n+

n(n − 1)

+

n(n − 1)

=

n(n2 + 3n − 1)

.

2

6

2

2

3

5.3. Теорема об LU -разложении. Найдем условия примени-

мости метода Гаусса. Заметим, что

b1′ = a11b1, b2′′ = (b2

− b1

a21

)/a22′ , . . . .

a11

b(ii) = c1b1 + c2b2 + · · · + cibi.

Теорема 5.1. Пусть все угловые миноры квадратной матрицы A отличны от нуля. Тогда имеет место единственное разложение

A = LU,

где L нижнетреугольная матрица с ненулевыми элементами на главной диагонали, U верхнетреугольная матрица с единицами на главной диагонали.

11

Доказательство. Докажем по индукции. Для n = 2

a b	a 0	1	e2	,
c d =	c e1	0	1

где e2 = ab , e1 = d − c · ab . Предположим утверждение справедливо для n − 1. Докажем для n. Пусть

		a11	a12	. . .	a1,n−1	.
An−1	=	a21	a22	. . .	a2n−1
		...	...	... ...
		an−1,1	an−1,2	. . . an−1,n−1

Тогда, по предположению индукции, An−1 = Ln−1Un−1. Мы полу-

чаем				¯
		Ln−1				Un−1	u¯
				0
¯	A =		¯l	ln 0¯			1 ,
		T	= (u1				¯
где l = (l1, l2, . . . , ln−1), u¯				, u2, . . . , un−1) 0 = (0, 0, . . . , 0). Век-
¯	находятся из соотношений
тора l, u¯
Ln−1u¯ = (a1n, a2n, . . . , an−1,n)					T	¯
						, lUn−1 = (an1, an2, . . . , an,n−1).
Число ln находится по формуле
						¯
			ln = ann − l · u¯.

Покажем теперь единственность. Пусть A = LU = L′U ′. Заметим, что обратная матрица к верхнетреугольной (нижнетреугольной) снова верхнетреугольная (нижнетреугольная) матрица. Более того, обратные матрицы к U и U ′ также содержат единицы на главной диагонали. Тогда U U ′−1 = L′L−1. Следовательно, U U ′−1 L′L−1 диагональные матрицы. Поскольку на диагонали матрицы U U ′−1 (а следовательно и матрицы L′L−1) стоят единицы, то U U ′−1 = L′L−1 = E. Отсюда, в силу единственности обратной матрицы, U = U ′ и L = L′.

5.4. Обращение матрицы. Будем обращать методом Гаусса. Т.е. решим систему уравнений

n
X	aik xkj = δij ,
k=1
где	(0, i =6 j .
δij =
	1, i = j

12

j=i+1

Сначала, напишем LU -разложение. На это требуется n(n2−1)

3

ствия. Таким образом, нам необходимо решить две системы нений

Ly¯i = ei,

U x¯i = y¯i.

дейурав-

Для решения второй системы, при каждом i, требуется n(n−1) дей-

2

ствия. Запишем подробно первую систему для произвольного i. Первые i − 1 уравнений выглядят так

Отсюда, y1i имеют вид

Получаем

	l21y1i + l22y2i = 0		.
	l11y1i =	0
	. . . . . .

li−1,1y1i + li−1,2y2i + · · · li−1,i−1yi−1,i = 0
		2i =	· · ·	=	y	i−1,i	= 0. Тогда оставшиеся уравнения
= y
			li+1,iyii + li+1,i+1yi+1,i = 0					.
						liiyii = 1
						. . . . . .

ln,iyii + ln,i+1yi+1,i + · · · ln,nyn,i = 0

1

yii = lii

Pj−1

yji = − k=i ljkyki . lii

Подсчитывая число действий получаем

n

1 + X (j − i + 1) = (n − i + 1)(n − i + 2) . 2

Суммируя по всем i получаем

1

n

n

(n − i + 1)(n − i + 2) = 1

(k(k + 1) = n(n + 1)(n + 2) .

X

X

2

i=1

2

k=1

6

Окончательно,

n(n2 − 1)

+

n2(n − 1)

+

n(n + 1)(n + 2)

= n3.

3

2

6

13

5.5. Метод прогонки. Рассмотрим один частный случай системы линейных уравнений. Пусть система линейных уравнений имеет трехдиагональный вид, т.е.

an−2

	x1 = µx2 + ω1
a1x1 − b1x2 + c1x3 = d1
a2x2	− b2x3 + c2x4	= d2	.
	. . . . . .

xn−2 − bn−2xn−1 + cn−2xn = dn−2 xn = νxn−1 + ω2

Будем искать решения в виде xi = αixi+1 + βi. Находим,

xi−1 = αi−1xi+βi−1 = αi−1(αixi+1+βi)+βi−1 = αi−1αixi+1+αi−1βi+βi−1.

Подставляем в уравнение

ai−1xi−1 − bi−1xi + ci−1xi+1 = di−1.

Получаем

ai−1(αi−1αixi+1 + αi−1βi + βi−1) − bi−1(αixi+1 + βi) + ci−1xi+1 = di−1.

Отсюда,

[ai−1αi−1αi−bi−1αi+ci−1]xi+1+[ai−1αi−1βi+ai−1βi−1−bi−1βi−di−1] = 0.

Для того, чтобы это равенство превратилось в тождество, выражения в квадратных скобках должны обращаться в ноль, т.е.

ai−1αi−1αi −bi−1 αi + ci−1 = 0, ai−1αi−1βi + ai−1βi−1 −bi−1 βi −di−1 = 0.

Отсюда,

αi =		ci−1	, βi =	ai−1βi−1 − di−1	.
	bi−1	− ai−1αi−1
				bi−1 − ai−1αi−1

Заметим, что α1 = µ, β1 = ω1. Таким образом, мы можем найти все αi.βi. Из последнего уравнения,

xn = ν(αn−1xn + βn−1) + ω2.

Тогда

xn = νβn−1 + ω2 . 1 − ναn−1

Теперь мы можем найти все xi.

14

6. Лекция 6

Пусть имеется система линейных уравнений

Ax = b.

Приведем эту систему к виду

xi = − i−1	aij	xj −	n	aij	xj +	bi	.
			jX
X	aii			aii		aii
j=1			=i+1

Перейдем к изучению итерационных методов.

6.1. Методы Якоби и Зейделя. В методе Якоби, итерации определяются следующим образом:

xim+1 = − i−1	aij	xjm −	n	aij	xjm +	bi	.
			jX
X						aii
j=1 aii			=i+1 aii

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

max |xim − xim−1| < ε.

1≤i≤n

Итерационный метод Зейделя имеет вид

xim+1 = − i−1

aij

xjm+1

−

n

aij

xjm +

bi

.

X

jX

aii

j=1 aii

=i+1 aii

Пусть

a11

a12

. . . a1,n

.

A = a21

a22

. . . a2n

... ...

... ...

an,1

an,2

. . . an,n

Рассмотрим разложение A = A1 + D + A2, где

0

0 . . . 0

a11

0

. . .

0

A1

=

a21

0 . . . 0 , D =

0

a22

. . .

0

,

...

... ... ...

...

...

... ...

an,1

an,2 . . . 0

0

0

. . .

an,n

		0	a12	. . . a1,n
A2	=	0	0	. . .	a2,n
		...	...	... ...
		0	0	. . .	0

15

¯

в виде

Представим систему Ax¯ = b

x¯ = −D

−1

A1x¯ − D

−1

A2x¯ + D

−1¯

b.

Тогда метод Якоби имеет вид

x¯

m+1

= −D

−1

A1x¯

m

− D

−1

A2x¯

m

+ D

−1¯

b.

Соответственно, метод Зейделя имеет вид

x¯

m+1

= −D

−1

A1x¯

m+1

− D

−1

A2x¯

m

+ D

−1¯

b.

Запишем эти выражения в другой форме:

m+1

− x¯

m

) + Ax¯

m

¯

D(¯x

= b;

(D + A1)(¯x

m+1

− x¯

m

) + Ax¯

m

¯

= b.

6.2. Норма матрицы.

Определение 6.1. Нормой матрицы A называется число

||A|| = sup |Ax¯| = sup |Ax¯|.

x¯ |x¯| |x¯|=1

Свойства нормы:

(1) ||A|| ≥ 0;

(2) ||αA|| = |α| · ||A||;

(3) ||A + B|| ≤ ||A|| + ||B||;

6.3. Методы Якоби и Зейделя. Сходимость итерационных методов. Теперь мы исследуем вопрос о сходимости итерационных методов.

Пусть нам дана система уравнений

¯

Ax¯ = b.

Перепишем эту систему в виде

x¯ = Bx¯ + c¯.

Тогда метод простой итерации записывается в виде x¯m+1 = Bx¯m + c¯.

Пусть r¯n = x¯n − x¯, где x¯ точное решение. Тогда

r¯m+1 = x¯m+1−x¯ = Bx¯m+¯c−x¯ = Bx¯m+¯c−(Bx¯+¯c) = B(¯xm−x¯) = Br¯m.

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

r¯m = Bmr¯0.

Получаем,

|r¯m| = |Bmr¯0| ≤ ||B|| · |r¯0|.

16

Тогда имеет место следующая теорема.

Теорема 6.2. Пусть ||B|| < 1. Тогда итерационный метод простой итерации сходится при любом начальном x¯0.

Рассмотрим более общий случай итерационного процесса. Пусть

¯

дана система линейных уравнений Ax¯ = b. Тогда рассмотрим следующую итерацию

	x¯m+1	− x¯m		m	¯
B	τ		+ Ax¯		= b.	(1)

Замечание 6.3. Заметим, что если это процесс сходится, то

m ¯

x¯ → x¯, где x¯ решение уравнения Ax¯ = b. Действительно, для любого ε > 0 существует такое N N, что для любого m > N |x¯m − x¯| < 2ε . Отсюда

|x¯m+1 − x¯m| = |(¯xm+1 − x¯) − (¯xm − x¯)| ≤ |x¯m+1 − x¯| + |x¯m − x¯| < ε.

Тогда

B	x¯m+1	− x¯m	¯
	τ		→ 0.

Погрешностью метода на m-й итерации называется r¯m = x¯m − x¯,

¯

где x¯ решение уравнения Ax¯ = b. Запишем уравнение (1) в виде

	(¯xm+1	− x¯) − (¯xm − x¯)				m	¯
B					+ A(¯x − x¯ + x¯) = b.
		τ
Отсюда,
		B	r¯m+1 − r¯m	+ Ar¯m = 0.
			τ
Для вещественной матрицы				C,		условие C	> 0 означает
(Cx,¯ x¯) > 0, x¯.

Утверждение 6.4. Пусть C > 0. Тогда существует δ > 0 такое, что (Cx,¯ x¯) ≥ δ|x¯|2.

Доказательство. Если C симметричная матрица, то в качестве δ можно взять минимальное собственное значение. Заметим, что (Cx,¯ x¯) = (¯x, C x¯), где C транспонированная матрица. Поскольку (¯x, y¯) = (y,¯ x¯) имеем

(Cx,¯ x¯) =	1	((Cx,¯ x¯)+(¯x, C x¯)) =	1	((Cx,¯ x¯)+(C x,¯ x¯)) = ((		1	(C+C ))¯x, x¯).
			2			2
2
Далее, утверждение следует из того, что матрица 1					(C + C ) сим-
метрична.				2

17

Теорема 6.5. Пусть A > 0 симметричная матрица. Пусть выполнено неравенство

B − 12 τ A > 0.

Тогда итерационный метод (1) сходится.

Доказательство. Покажем, что числовая последовательность Km = (Ar¯m, r¯m) является не возрастающей. Заметим, что

r¯m+1 = (E − τ B−1A)¯rm.

Тогда

Ar¯m+1 = (A − τ AB−1A)¯rm.

Отсюда, поскольку (Ax,¯ y¯) = (¯x, Ay¯),

(Ar¯m+1, r¯m+1) = ((A − τ AB−1A)¯rm, (E − τ B−1A)¯rm) = (Ar¯m, r¯m)− τ (Ar¯m, B−1Ar¯m) − τ (AB−1Ar¯m, r¯m) + τ 2(AB−1Ar¯m, B−1Ar¯m) = (Ar¯m, r¯m) − 2τ (Ar¯m, B−1Ar¯m) + τ 2(AB−1Ar¯m, B−1Ar¯m) =

(Ar¯m, r¯m) − 2τ ((A − τ2 AB−1A)¯rm, B−1Ar¯m) =

(Ar¯m, r¯m) − 2τ ((BB−1A − τ2 AB−1A)¯rm, B−1Ar¯m) =

(Ar¯m, r¯m) − 2τ ((B − τ2 A)B−1Ar¯m, B−1Ar¯m).

Из условия теоремы ((B − τ2 A)B−1Ar¯m, B−1Ar¯m) > 0. Следовательно, (Ar¯m+1, r¯m+1) < (Ar¯m, r¯m). С другой стороны, (Ar¯m, r¯m) > 0 для любого m. Тогда существует

K = lim Km.

m→∞

Из утверждения 6.4 следует, что существует δ > 0 такое, что

((B − τ2 A)B−1Ar¯m, B−1Ar¯m) ≥ δ|B−1Ar¯m|2.

Отсюда,

Km+1 = Km − 2τ ((B − τ2 A)B−1Ar¯m, B−1Ar¯m) ≤ Km − 2τ δ|B−1Ar¯m|2

Получаем,

Km+1 − Km + 2τ δ|B−1Ar¯m|2 ≤ 0.

Переходя к пределу, получаем

q¯m = B−1Ar¯m → 0.

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

r¯m = A−1Bq¯m → 0.

18

Следствие 6.6. Пусть A > 0 симметричная матрица.

Пусть выполнено условие
aii >	|aij |.
j=6 i
X

Тогда итерационный метод Якоби сходится.

Доказательство. Условие теоремы 6.5, в этом случае, запи-

шется A < 2D (т.е. 2D − A > 0). Запишем
	X
(Ax,¯ x¯) =			aij xixj .
		i,j
Воспользуемся оценкой |xixj | ≤ 21 (xi2 + xj2). Получаем
X		1	X		X
(Ax,¯ x¯) ≤ |aij ||xixj | ≤			( |aij |xi2	+	|aij |xj2).
i,j	2		i,j		i,j

Поскольку aij = aji, имеем

(Ax,¯ x¯) ≤ i=1	j=1 |aij |! xi2
n	n

X X

n

X

< 2aiix2i

i=1

n

!

XX

=	aii + |aij | xi2 <
i=1	j=6 i

= (2Dx,¯ x¯).

Следствие 6.7. Пусть A > 0 симметричная матрица. Тогда

метод верхней релаксации
	x¯m+1	− x¯m		m	¯
(D + τ A1)	τ		+ Ax¯		= b

сходится при 0 < τ < 2. В частности, сходится метод Зейделя (τ = 1).

Доказательство. Так как матрица A симметрична, то

(Ax,¯ x¯) = (Dx,¯ x¯) + (A1x,¯ x¯) + (A2x,¯ x¯) = (Dx,¯ x¯) + (A1x,¯ x¯) + (¯x, A1x¯)

		= (Dx,¯ x¯) + 2(A1x,¯ x¯).
Тогда условие теоремы 6.5 примет вид
(((D + τ A1) −	1	τ A)¯x, x¯) = ((D + τ A1)¯x, x¯) −			1	τ (Ax,¯ x¯) =
	2				2
		1					τ
(Dx,¯ x¯) + τ (A1x,¯ x¯) −				τ ((Dx,¯ x¯) + 2(A1x,¯ x¯)) = (1 −				)(Dx,¯ x¯).
			2				2

19