[ Буслов, Яковлев ] Часть 2 - Численные методы

áûë k-й шаг) переставляют строки с номерами k ; k + 1 ; : : : ; N таким образом, чтобы на месте kk оказался элемент a(mkk) , наибольший из всех в k-ом столбце при m > k (при этом, естественно, переставляются и компоненты вектора b).

Можно для максимальной точности переставлять также и столбцы преобразуемой матрицы, чтобы на месте kk оказался максимальный элемент из всех с индексами больше, либо равными k. Эта процедура

называется методом Гаусса с выбором главного элемента. Она несколько повышает точность по сравнению с частичным выбором главного элемента, но весьма неудобна, в том числе и для программирования, поскольку при перестановке строк компоненты искомого вектора x переставлять не надо, тогда как при

перестановке столбцов надо переставлять и соответствующие компоненты вектора x.

Опишем обратный ход метода Гаусса в несколько иной форме (треугольное разложение). Введем матрицы Mk по правилу

	0	1	0	¢ ¢ ¢		0
	B	0	1	¢ ¢ ¢		0
	B			...
	B			¢ ¢ ¢
M = B		¢ ¢ ¢	¢ ¢ ¢		¢ ¢ ¢
	B
k	B	0	0	¢ ¢ ¢		1
	B
	B
	B
	B	0	0		ck+1;k
	B
	B	0	0
	B				ck+2;k
	B			¢ ¢ ¢
	B
	B	¢ ¢ ¢	¢ ¢ ¢	¢ ¢ ¢	¢ ¢ ¢
	B	0	0	¢ ¢ ¢	c
	@					N;k
	B

¢ ¢ ¢	0	1
¢ ¢ ¢	0	C
		C
		C
¢ ¢ ¢	0	C	:
¢ ¢ ¢	¢ ¢ ¢	C
	0	C
¢ ¢ ¢		C
		C
	0	C
¢ ¢ ¢		C
		C
...		C
		C
		C
	¢ ¢ ¢	C
	1	C
		C
¢ ¢ ¢		A

На каждом шаге метода Гаусса получается некоторая промежуточная матрица Ak+1 =MkMk¡1 : : : M1A , и вектор fk+1 = MkMk¡1 : : : M1b . Нетрудно видеть, что

	N¡1	N¡1	N
U =	á MiA ; f =	á Mib ; Ux = f ; det U =	Uii = det A :
	Y	Y	iY
	i=1	i=1	=1

Вопрос. Почему det U = det A?

Если производить также выбор главных элементов, то необходимо использовать оператор P переста-

новки индексов l è m, матричные элементы которого равны: pij = 0 , i; j =6 l; m ; pim = pmi = 0 , i =6 l ; pli = pil = 0 , i =6 m ; pml = plm = 1 . При применении оператора перестановки индексов к матрице слева, меняются местами строки матрицы и компоненты свободного вектора (P Ax = P b) , если же его применить

справа к матрице, то меняются местами ее столбцы и компоненты решения (A P P x = b).

|{z}

=I

1.2.3 L-R разложение

Для решения задачи Ax = b несколько модифицируем ее. Именно, введем N £ (N + 1) матрицу

0	A		b1	1
C = B			b2	C
B			.	C
B				C
B			bN	C
B				C
B				C
@				A

11

и вектор

X = (x1; x2; : : : ; xN ; ¡1)T

размерности (N + 1), тогда исходная задача эквивалентна следующей

CX = 0 :

Представим C â âèäå C = LR, ãäå L нижнетреугольная

N £ N матрица

l11

0

0

1

0 l21

l22

¢ ¢ ¢

0

L = B .

.

¢.¢..¢ .

C ;

B

C

B

C

BlN1 lN2

lNN C

N £ (N + 1)-матрица:

@

¢ ¢ ¢

A

à R

0

1

1

r12

r13

¢ ¢ ¢

r1N

r1;N+1

0 1 r23

r2N

r2;N+1

R = B

0 0

1

¢ ¢ ¢

r3N

r3;N+1

C :

B

¢ ¢ ¢

C

B

.

.. .

.

C

B . . .

C

B

C

B

0

0

0

1

C

B

rN;N+1 C

@

¢ ¢ ¢

A

Как находить матрицы L è

R?

1-й шаг) а) Умножим каждую строку матрицы

L

на первый столбец матрицы R, откуда li1 = ci1. Таким

образом, мы определили первый столбец матрицы

L.

б) Умножим первую строку

L

на каждый столбец

R, откуда r1i = c1i=l11, то есть определена

первая строка R.

2-й шаг) a) Умножим каждую строку

L (начиная со второй) на второй столбец R и определим второй

столбец

L: li2 = ci2 ¡ li1r12.

L на каждый столбец R, определяем вторую строку R: r2i =

б) Умножая вторую строку

(c2i ¡ l21r1i)=l22.

m ¡ 1

m ¡ 1

m-й шаг) Пусть известны первые

столбец

L è

строка

R, тогда при i ¸ m

lim = cim

likrkm ; rmi =

kP

:

m¡1

m¡1

cmi ¡

lmkrki

X

=1

¡

lmm

k=1

Теперь заметим, что вовсе нет необходимости решать задачу CX = 0, а достаточно решить систему RX = 0. Действительно, ранг матрицы R равен N, таким образом исходная матрица A è L вырождены или невырождены одновременно. Компоненты xi находим последовательно, начиная с N-îé:

	N
xN = rN;N+1 ; xi = ri;N+1 ¡	kX
	rikxk :
	=i+1

Вычисления по изложенному методу требуют в два раза меньший объем памяти, чем по методу Гаусса.

1.2.4 Метод прогонки

Пусть A трехдиагональная матрица, которую мы представим в виде:

12

0

c1

¡b1

0

0

0

¢ ¢ ¢

0

0

1

A =

B

¡a2

c2

¡b2

0

0

¢ ¢ ¢

0

0

C

B

0

¡

a3

c3

¡

b3

0

¢ ¢ ¢

0

0

C

B

C

B

C

B

¢ ¢ ¢

¢ ¢ ¢

¢ ¢ ¢

¢ ¢ ¢

¢ ¢ ¢

¢ ¢ ¢

¢ ¢ ¢

¢ ¢ ¢

C

B

C

B

0

0

0

0

0

aN

cN

C

B

¢ ¢ ¢

C

B

¡

C

@

A

Çíàê ¡ перед bi ; ci поставлен для удобства. Для решения задачи

At = s

в этом случае применяется

метод прогонки.

Положим a1 = bN = 0 , тогда трехдиагональная система может быть записана в виде

¡tk¡1ak + tkck ¡ tk+1bk = sk ; k = 1 ; 2 ; : : : ; N :

Рассмотрим эту систему подробнее. Выразим из первого уравнения

t1

через

t2 :

b1

s1

t1c1 ¡ t2b1 = s1 ) t1 =

t2 +

:

c1

c1

Теперь из второго уравнения выразим t2

через t3 : ¡t1a2 + t2c2 ¡ t3b2 = s2, èëè

s

1

b

b

t

3

s2

+ a2 s1

¡ µ

1

t2¶a2 + t2c2

2

c1

+

¡ t3b2 = s2 ) t2 =

+

:

c1

c1

c2 ¡ cb11 a2

c2 ¡ cb11 a2

Аналогично,

tk = ®ktk+1 + ¯k ;

ãäå

bk

sk + ¯k¡1ak

®k =

; ¯k =

:

ck ¡ ®k¡1ak

ck ¡ ®k¡1ak

Убедимся в справедливости этого представления по индукции. Действительно, ®1 = b1 ; ¯1

= s1

c1

c1 , таким

образом база индукции верна. Теперь осуществим собственно индукционный переход. Пусть tk = ®ktk+1 + ¯k , тогда

¡ak+1tk + ck+1tk+1 ¡ bk+1tk+2 = sk+1 ;

¡ak+1(®ktk+1 + ¯k) + ck+1tk+1 ¡ bk+1tk+2 = sk+1 ;

	откуда	bk+1tk+2		sk+1 + ¯kak+1
	tk+1 =		+		= ®k+1tk+2 + ¯k+1 ;
		ck+1 ¡ ®kak+1
				ck+1 ¡ ®kak+1

то есть индукционный переход также имеет место.

Рассмотрим теперь каким образом применяется метод прогонки. На первом этапе (прямой ход прогонки) мы определяем коэффициенты ®k ; ¯k через известные нам элементы матрицы A (bk ; ck ; ak) , заданные

значения sk и предыдущие			®k¡1 ; ¯k¡1:
®1 = b1 ;		¯1 = s1			; начало прямого хода,
	c1			c1
®k =		bk		;	¯k = sk+¯k¡1ak ;	прямой ход.
	ck¡®k¡1ak				ck¡®k¡1ak
После того как определены коэффициенты ®k						è ¯k начинается обратный ход прогонки собственно

определение компонент tk . Имеем

tN = ®N tN+1 + ¯N ;

ïðè ýòîì ®N = 0 , ò.ê. bN = 0 , à ®N = bN

cN ¡®N¡1aN . Таким образом,

13

tN = ¯N (начало	обратного	õîäà) ,
tk = ®ktk+1 + ¯k	(обратный	õîä) .

Утверждение (Достаточное условие разрешимости прогонки): Пусть jckj > jbkj + jakj , k = 1 ; : : : ; N , тогда det A 6= 0.

Доказательство. Необходимо убедиться, что знаменатель в формулах прямого хода не обращается в нуль. Для этого достаточно убедиться в том, что j®kj < 1. Âåäü åñëè ýòî òàê, òî

jck ¡ ®k¡1akj ¸ jckj ¡ j®k¡1jjakj > jckj ¡ jakj > jbkj ¸ 0

и не происходит деления на нуль. Имеем :

®

=

b1

< 1 ; ®

=

jbkj

<

jbkj

= 1 :

j

jc1 j

jck ¡ ®k¡1akj

jbkj

1j

j kj

1.2.5 Метод итераций для решения линейных систем

Система линейных уравнений Ax = b :

N
Xj
aijxj = bi ; i = 1; 2; : : : ; N ;	(1)
=1

может быть решена не только прямыми методами, но также и итерационными. Разумеется, мы предполагаем, что система имеет единственное решение, т.е., что det A =6 0.

Представим матрицу A â âèäå A = B+D , ãäå D = diagfa11 ; : : : ; aNN g . Предположим, что det D =6 0 , что равносильно тому, что aii =6 0 ; i = 1 ; : : : ; N (если исходно это не так, то перестановкой строк и столбцов этого всегда можно добиться при det A 6= 0 ). Тогда (1) переписывается в виде Dx = b ¡ Bx, èëè

x = D¡1b ¡ D¡1Bx :

Предложим следующую итерационную процедуру

xs+1 = D¡1b ¡ D¡1Bxs ;

x0 произвольный начальный вектор. В развернутой форме

xis+1 = aii¡1bi ¡ aii¡1	n
	aijxjs ; i = 1; 2; : : : N :
	=1;j
j	X6=i
Обозначим D¡1b = u, D¡1B = T , тогда итерационный процесс принимает вид
xs+1 = u ¡ T xs :		(2)

Теорема1. Процесс (2) сходится, для любого начального вектора, если jjD¡1(A ¡ D)jj = jjT jj < 1 :

Доказательство. Для доказательства достаточно заметить, что отображение x ! u ¡ T x является сжа-

òèåì.		T x¤, èëè,
Таким образом, последовательность xs имеет предел. Пусть x¤ = lim xs, тогда x¤ = u	¡
s!1

возвращаясь к исходной формулировке, Ax¤ = b. Итак, для сходимости изложенного метода, называемого

методом простых итераций, необходимо чтобы

jjD¡1(A ¡ D)jj < 1 :

14

1.2.6 Метод Зейделя

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

s+1	1					s	s	] ; i = 1; 2; : : : ; N :
xi	= aii¡	[bi ¡ aijxj					¡ aijxj
		j<i					j>i
		X					X
		¤
		|		{z		}

Заметим, что если последовательно вычислять компоненты (s + 1)-го приближения xs+1, начиная с первой

xs+1	i-ой компоненты xs+1	xs+1	; : : : ; xs+1
1 , то к моменту вычисления конкретной	i	, координаты 1	i¡1 уже опреде-

лены и их можно было бы использовать для определения более точного последующего приближения xs+1. Модифицируем соответствующим образом метод простых итераций, заменив в сумме ¤ компоненты xsj íà

xjs+1	. Таким образом, мы получаем новую итерационную процедуру
	xis+1	= aii¡1bi ¡ aii¡1	X	aijxjs+1 ¡ aii¡1	X
					aijxjs ; i = 1; 2; : : : ; N :
			j<i		j>i

Такой итерационный процесс называется методом Зейделя. Представим его в матричной форме. Пусть Lнижнетреугольная матрица с элементами

L :		8 aij
à U верхнетреугольная матрица с	lij =			; j < i ;
		:	0	; j ¸ i
		<
	элементами
	uij =	8 aij		; j > i :
		:	0		·
		<		; j		i

Как и раньше введем матрицу D = diagfa11 : : : aNN g , тогда A = D + L + U . В матричном виде метод

Зейделя имеет вид:

xs+1 = D¡1b ¡ D¡1Lxs+1 ¡ D¡1Uxs :

Сходимость метода Зейделя

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

Dxs+1 = b ¡ Lxs+1 ¡ Uxs ;

èëè

xs+1 = (D + L)¡1b ¡ (D + L)¡1Uxs :

Отображение x 7!(D + L)¡1b ¡ (D + L)¡1Ux является сжатием, если jj(D + L)¡1Ujj < 1 , òî åñòü

справедлива

Теорема. Метод Зейделя сходится, если jj(D + L)¡1Ujj < 1.

Условия этой теоремы довольно трудно проверяемы, так как матрица (D + L)¡1U должна еще и вычисляться. Существует достаточно простой признак сходимости метода Зейделя, который связан с понятием

15

называется положительно определенным, если

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

hAx; xi ¸ °hx; xi ; ° > 0 :

Если оператор положительно определен, то у него существует и обратный и он также положительно определен. Также важно отметить, что если оператор A положительно определен и симметричен в RN , òî

форма

hx; yiA = hAx; yi

удовлетворяет всем свойствам скалярного произведения. В дальнейшем, факт положительной определенности оператора A будем обозначать: A > 0 . Заметим, что в комплексном евклидовом пространстве

положительная определенность оператора A автоматически влечет за собой эрмитовость: A = A¤.

Теорема (достаточный признак сходимости метода Зейделя). Метод Зейделя сходится в вещественном евклидовом пространстве, если A симметричная положительно определенная матрица.

Для доказательства этой теоремы нам потребуется следующая

Лемма. Пусть последовательность векторов zk â RN определена рекуррентным соотношением

B(zk+1 ¡ zk) + Azk = 0 ;

(3)

ãäå B ¡ 21 A > 0 , A > 0 и симметрична, тогда

zk ! 0 .

Доказательство. Представим zk â âèäå

1

1

zk =

(zk+1

+ zk) ¡

(zk+1 ¡ zk) ;

2

2

и подставим это представление в (3), тогда

1

1

B(zk+1 ¡ zk) +

A(zk+1 + zk) ¡

A(zk+1 ¡ zk) = 0 ;

2

2

èëè

1

1

(B ¡

A)(zk+1 ¡ zk) +

A(zk+1 + zk) = 0 :

2

2

Умножим это равенство скалярно на

zk+1 ¡ zk , тогда

1

0 = jzk+1 ¡ zkjB¡A=2 +

hA(zk+1 + zk); zk+1 ¡ zki =

2

1

= jzk+1 ¡ zkjB¡A=2 +

fjzk+1jA ¡ jzkjAg = 0 ;

2

ãäå j ¢ jA = hA¢; ¢i ; j ¢ jB¡A=2 = hfB ¡ A=2g¢; ¢i

нормы, определяемые операторами A è B ¡ A=2, ñîîò-

ветственно. Из последнего равенства в силу положительной определенности оператора (B ¡ A=2) следует,

÷òî jzk+1jA ¡jzkjA · 0 , т.е. последовательность

jzkjA невозрастающая: jzk+1jA · jzkjA . При этом, после-

довательность чисел jzkjA		ограничена снизу поскольку jzkjA ¸ 0 . Таким образом, существует конечный
предел	klim jzkjA = a . Но тогда из того же равенства следует, что норма jzk+1 ¡ zkj(B¡21 A) стремится к
	!1	! 0 ; k ! 1 . Вернемся теперь к определению последовательности zk :
нулю, а значит и zk+1 ¡ zk
		Azk = ¡B(zk+1 ¡ zk) ;
откуда	zk = ¡A¡1B(zk+1 ¡ zk) и, следовательно,
		jjzkjj · jjA¡1Bjj £ jjzk+1 ¡ zkjj ! 0 ;

16

и, таким образом, zk ! 0, ïðè k ! 1.

Приступим теперь собственно к доказательству достаточного признака сходимости метода Зейделя. Как нетрудно видеть, метод Зейделя (D + L)xs+1 + Uxs = b может быть представлен в виде

(D + L)(xs+1 ¡ xs) + Axs = b :

Пусть u точное решение уравнения Au = b , оно существует, так как A положительно определенный оператор и, следовательно, обратим. Положим также zs = xs ¡ u , тогда

(D + L)(zs+1 ¡ zs) + Azs = 0 :

Убедимся в том, что (D + L ¡ 12 A) положительно определенная матрица, если A симметрична и положительно определена. Действительно

D + L ¡ 12A = D + L ¡ 12(D + L + U) = 12(D + L ¡ U) :

Рассмотрим соответствующую квадратичную форму

h(D + L ¡ U)x; xi = hDx; xi + hLx; xi ¡ hUx; xi :

Заметим, что поскольку A симметричная матрица, следовательно LT = U è

hLx; xi = hx; LT xi = hx; Uxi = hUx; xi ;

поэтому
h(D + L ¡ U)x; xi = hDx; xi = Xi	aiixi2 > 0 ;

так как у положительно определенной матрицы все диагональные элементы больше нуля (почему?): aii > 0 . Таким образом, мы находимся в условиях Леммы, следовательно, последовательность zs стремится к нулю, откуда следует, что последовательность xs = u + zs стремится к истинному решению u .

17

hei; e˜ji = ±ij

.

e1; e2; : : : ; eN

(базиса) суще-

Глава 2

Алгебраические спектральные задачи

2.1 Некоторые сведения из матричной теории

Пусть A линейный оператор, действующий в вещественном RN или в комплексном CN евклидовом пространстве: A : RN (CN ) ! RN (CN ) .

Число ¸ и вектор x называются, соответственно, собственным числом (значением) и собственным

вектором оператора A, отвечающим собственному числу ¸ , åñëè Ax = ¸x .

В частности, справедливы следующие теоремы.

Теорема 1.Всякий линейный оператор в CN имеет по крайней мере одно собственное значение.

Теорема 2. Собственные векторы, отвечающие различным собственным значениям, линейно независимы.

Теорема 3. Для любого набора из N линейно независимых векторов ствует единственный дуальный базис e˜1; e˜2; : : : ; e˜N , такой что

Заметим, что всякий ортонормированный базис самодуален. Пусть A имеет N различных собственных векторов xi , тогда они образуют базис, и, следовательно, существует дуальный базис x˜i . В этом случае, как нетрудно убедиться, сопряженный оператор A¤ (в случае вещественного евклидова пространства

просто транспонированная матрица A	T	) имеет в качестве собственных значений числа	¯
			¸i , а в качестве

собственных векторов векторы дуального базиса:

i i ¤˜i ¯ ˜i

Ax = ¸ix ; A x = ¸ix :

i ¯ ˜i i ˜i i ˜i i ¤˜i

Действительно: hx ; ¸ix i = h¸x ; x i = hAx ; x i = hx ; A x i и, аналогично,

j ¤˜i j ¯˜j

åñòü hx ; A x i = ±ijhx ; ¸x i : Кроме того, нетрудно показать справедливость

разложения оператора A :

A¢ = XN ¸ih¢; x˜iixi = XN ¸iPi¢ ;

hxj; A¤x˜ii = 0 ; i 6= j , òî

следующего спектрального

i=1 i=1

где операторы Pi¢ = h¢; x˜iixi суть собственные проекторы оператора A . В самом деле, произвольный вектор f можно разложить по собственным векторам оператора A : f = Phf; x˜iixi . Тогда

Af = Xhf; x˜iiAxi = Xhf; x˜ii¸ixi = X¸iPif :

18

Пусть A = A¤ эрмитова матрица (в RN симметричная). В этой ситуации собственные значения

вещественны, алгебраическая и геометрическая кратности любого собственного значения совпадают, соб- ственные векторы xi , отвечающие различным собственным числам ортогональны, и существует ортонор-

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

отвечающий ему собственный проектор P одномерен и имеет вид P = h¢; xix (всюду считаем, что собствен-

ный вектор x нормирован на единицу). Если подпространство решений Ax = ¸x более чем одномерно, то в нем выбирается произвольный ортобазис xi и собственный проектор, отвечающий собственному числу

¸, представляет собой сумму соответствующих одномерных проекторов P = Ph¢; xiixi . Отметим (легко проверяемое) важное свойство ортогональных проекторов:

PiPk = ±ikPi :

Степень оператора имеет следующую запись через ортогональные проекторы

X
Am = ¸m(A)P	k	:
k
k

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

X

f(A) = f(¸k)Pk :

k

Собственные функции у оператора и у функции от оператора совпадают, тогда как собственные значения функции от оператора есть числа f(¸k).

2.2 Собственные числа эрмитовых матриц

2.2.1 Интерполяционный метод

Поскольку собственные числа ¸i матрицы A являются корнями характеристического полинома FA(¸) = det(A ¡ ¸I) ; то можно вычислить FA(¸) â (n + 1)-ом наугад выбранном значении ¸ (их естественно выбирать в промежутке (¡jjAjj; jjAjj), если границы спектра известны; оценить их можно по максимальному

по модулю элементу матрицы) и построить по ним интерполяционный полином степени n, который совпа-

дает собственно с характеристическим, после чего определяются его корни. Этот метод применим и для неэрмитовых матриц (при соответствующем выборе метода определения корней).

2.2.2 Нахождение максимального по модулю собственного значения

Для удобства будем считать, что собственные числа пронумерованы в порядке убывания их модуля.

а) Метод итераций

Пусть g = g0 произвольный начальный вектор. Определим последовательность

	gn = A	gn¡1		=		Ang		;
		jjg(n¡1)jj			jjAn¡1gjj

тогда

lim jjg(n)jj = j¸maxj :

19

одномерно, то

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

g = Pn

Pkg ; kjng

n

jjAngjj

n

jj = jjAn¡1gjj , è

n

X

X

A g =

¸k Pkg = ¸max(Pmaxg +

f¸k=¸maxg Pkg) :

k6=1

Пусть ¸0 собственное число, следующее за максимальным по модулю. Тогда

jj

Ang

jj

2

=

h

Ang; Ang

i

= ¸2n (

P

max

g; g

i

+ O([¸0

=¸

max

]2n)) ;

max

h

è

=

jj

Ang

jj

=

h

Pmaxg; g + O([¸0=¸max]2n)

=

gn

¸

jj

jjAn¡1gjj

j

maxjshPmaxg; gii+ O([¸0=¸max]2n¡2)

jj

= j¸maxjf1 + O([¸0=¸max]2n¡2)g :

Таким образом, если стартовый вектор g имел ненулевую проекцию на собственное подпространство,

отвечающее максимальному по модулю собственному значению (то есть Pmaxg =6 0 ), то приведенная итерационная процедура приводит к нахождению ¸max . Однако, хотя формально, предыдущее рассмотрение

верно лишь в случае ненулевой проекции, в действительности, из-за ошибок округления при вычислениях эта проекция наверняка появится на некотором шаге и дальнейшее применение метода итераций приведет к желаемому результату. Попутно заметим, что если подпространство отвечающее ¸max

метод итераций одновременно приводит к нахождению собственного вектора xmax , отвечающего ¸max . Этим вектором с точностью до нормировки является

xmax = lim gn : n!1

Замечание. Для нахождения ¸max можно применять метод итераций и в более простой постановке. Пусть l-ая компонента в максимального собственного вектора в стандартном евклидовом базисе не равна нулю (хотя бы одна такая существует), тогда

¸max = lim (Ang)l :

n!1 (An¡1g)l

б) Метод следов

Известно, что след матрицы (сумма диагональных элементов) равен сумме е¼ собственных значений с

P P

учетом кратности: ¸i = T rA , таким образом, ¸mi = T rAm и, следовательно,

T rAm = ¸mmax[1 + (¸0=¸max)m + : : :] ;

ãäå ¸0 следующее по модулю за максимальным собственное значение. Таким образом, ¸max можно искать как следующий предел

j¸maxj = mlim	m
	pT rAm ;
!1

	или, например, в виде	T rAm+1
	¸max = lim		:
		T rAm
	m!1

Процедуру возведения матрицы в степень можно оптимизировать:

A £ A £ A £ A ;

| {z } | {z }

|	A2	{z		A2	}
		A4

и так далее, в частности: A16 = (A8)2 = A2222 .

20