ГЛАВА 1. ОШИБКИ ВЫЧИСЛЕНИЙ

1.1 Источники и классификация погрешностей

В настоящее время разработана технология исследования сложных проблем, основанная на построении и анализе с помощью ЭВМ математических моделей изучаемого объекта. В соответствии с ней решение большинства практических задач включает следующие этапы:

1)создание математической модели задачи;

2)интерпретация математической модели в виде, доступном для реализации на ЭВМ (построение численного метода решения);

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

Этот процесс неизбежно носит приближенный характер, так как на каж-

дом этапе вносятся те или иные погрешности. Они обусловлены следующими причинами:

а) математическое описание задачи является упрощенным; б) недостаточно точно заданы исходные данные, являющиеся, как правило,

результатом проведенных экспериментов; в) всякий численный метод, используемый для решения математической за-

дачи, является приближенным в том смысле, что он не дает точного решения задачи, (для точного решения, например, может потребоваться неограниченно большое число арифметических операций);

г) в самом процессе решения задачи в ЭВМ при вводе исходных данных, выполнении арифметических операций и при выводе результатов производятся округления.

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

ции и на возникающую в процессе решения (пункт «г») вычислительную по-

грешность или погрешность округления.

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

6

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

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

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

1.2. Приближенные числа, их абсолютные и относительные погрешности

Точные и приближенные числа. При численном решении задач приходится оперировать двумя видами чисел - точными и приближенными. К точным относятся числа, которые дают истинное значение исследуемой величины. К приближенным относятся числа, близкие к истинному значению, причем степень близости определяется погрешностью вычислений. Например, в утверждениях "в аудитории 102 слушателя", "в январе 31 день", "стоимость костюма 75 долларов" числа 102, 31 и 75 - точные. В высказываниях типа "сто метров спортсмен пробежал за 10,1 сек.", "масса мороженого 100 грамм", "расстояние от Бреста до Москвы равно 1054 км." числа 10,1; 100; 1054 - приближенные. Это связано с тем, что при измерениях появление приближенных чисел связано с несовершенством измерительных приборов. Кроме того, часто нет необходимости знать точное значение интересующей величины. Ясно, что пассажиру, едущему из Бреста в Москву неважно, что данное в справочниках расстояние в 1054 км. от Бреста до Москвы может быть короче на 200 м или длиннее на

150 м.

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

Если x - точное, вообще говоря, неизвестное значение некоторой величины, а x$ - его приближение, то разность x − x$ называется ошибкой, или погрешностью приближения. Часто знак ошибки x − x$ неизвестен, поэтому используется так называемая абсолютная погрешность ∆(x$) приближенного числа x$ ,

определяемая равенством
∆(x$) =		x − x$		,	(1.1)
откуда имеем
x = x$ ± ∆(x$).					(1.2)

Изучаемая числовая величина есть величина именованная, т.е. определяемая во вполне определенных единицах измерения, например, в сантиметрах, килограммах и т.п. Погрешность (1.1) имеет ту же размерность. Однако часто

7

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

	∆(x )						x − x
δ(x$) =	$					=	$	.	(1.3)
			$				$
			x				x

Относительную погрешность часто выражают в процентах: δ(x$) = ∆(x$x$) 100% .

Интервалы приближения. Поскольку точное значение x определяемой величины неизвестно, то и погрешности ∆(x$) и δ(x$) в соотношениях (1.1) и (1.3)

тоже неизвестны. Однако,						бывают известны оценки ∆ и δ сверху и снизу этих
величин, т.е. для чисел ∆ и δ выполнены неравенства.
						( $)	≤ ∆ , δ		( $)	≤ δ .	(1.4)
						∆ x			x
В этом случае говорят, что абсолютная и относительная погрешности
приближенного					числа x$	не превосходят ∆ и δ					соответственно. Если
∆(x$) =		x − x$		≤ ∆ ,	то говорят, что число x$			является приближением числа x с точ-
ностью ∆ . В этом случае, согласно (1.1) и (1.4),											(1.5)
						x − ∆ ≤ x ≤ x + ∆.
						$			$

Числа ∆ и δ называются предельными абсолютной и относительной погрешностями приближенного числа.

Практическое применение изложенного состоит в нахождении для числа x приближения ∆ с недостатком и приближения ∆ с избытком, для которых ∆ ≤ x ≤ ∆ . В этом случае говорят, что имеется интервал приближения величины x. Согласно (1.5), если известна предельная абсолютная погрешность ∆ , то границы интервала приближения имеют вид: ∆ = x$ − ∆ , ∆ = x$ + ∆.

Предположим, что имеется два интервала приближения величины x: ∆1 ≤ x ≤ ∆1 и ∆2 ≤ x ≤ ∆2 . Тогда можно, очевидно, без потери информации заме-

нить эти два интервала приближения на один: max(∆1 , ∆2 )≤ x ≤ min(∆1 , ∆2 ).

Запись приближенных чисел. Верные значащие цифры. Запись при-

ближенных чисел должна подчиняться правилам, связанным с понятиями о верных значащих цифрах.

Первая слева отличная от нуля цифра числа x и все расположенные справа от нее, называются значащими. Например, числа 0.00028745 и 200.37500 имеют соответственно 5 и 8 значащих цифр. Цифра ai числа, записанного в десятич-

ной системе, называется верной, если ∆(x$) <10i , т.е если абсолютная погрешность числа x$ не превосходит одной единицы соответствующего разряда десятичного числа. Цифры, не являющиеся верными, называются сомнительными. Если верная цифра ai - значащая, то она называется верной значащей цифрой. Вид записи приближенного числа должен показывать его абсолютную погрешность, которая не должна превосходить единицы последнего разряда, сохраняемого при записи. Другими словами, приближенные числа принято записы-

8

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

Пример 1. Числа 27,0563 и 27,0563000, как приближенные, различны.

Абсолютная погрешность первого числа не превосходит 10−4 (i = −4), а второго

10−7 (i = −7).

Пример 2. Числа x1=34,174562 и x2=375,16342 содержат пять верных

цифр. Определить их абсолютную погрешность.
	∆ В числе x1 последняя верная значащая цифра 4, считая слева, находится
на третьем месте					после	запятой, т.е.		i=-3. Следовательно, по определению
∆(x$1 )=10−3 = 0,001. Во втором числе x2								последняя значащая цифра 6 имеет ин-
декс i=-2. Следовательно,						$	−2	= 0,01 .
						∆(x 2 )=10
	Пример 3. Абсолютная погрешность чисел x1=28,9356 и x2=183,45123 со-
ставляет			∆ = 0,03. Определить, какие цифры являются верными, и округлить
числа, оставив только верные цифры.
	∆ При заданной максимальной абсолютной погрешности ∆ = 0,03 верной
значащей			цифрой, по			определению, должна быть та			ai , для	которой
ˆ	≤	10	−i ≤	0,03, т.е. i=-2. Значит, в заданных числах все цифры являются вер-
∆(x)
ными. Округлив					эти числа до указанной точности,				получим	x1 = 28,94 ,
										$

x$2 =183,45.

Замечание. При записи чисел имеются определенные тонкости. Например, запись числа 5,1813 означает, что абсолютная погрешность этого приближенного числа не превосходит 0,0001, а запись числа 680 - что абсолютная погрешность его не превосходит единицы. Если же второе число имеет большую точность, следует писать 680,0. Тем самым числа 68 10; 680; 680,0; 680,00 имеют различную степень точности. Если надо показать, что у числа 680,00 первые три цифры верные значащие, то его следует писать в виде 680, а для числа

600028 с тремя верными слева значащими цифрами следует писать 600 103 . В соответствии с этим округлением, например, число 28065 при 1, 2 или 3 верны-

ми значащими цифрами будет записана в виде 3 104 , 28 103 , 280 102 соответственно.

Связь между числом верных знаков и погрешностью числа. Пусть число x$ является приближенным значением точного числа x и имеет вид

x$ = an 10n + an−1 10n−1 +K+ai 10i +K+a−m 10−m , an ≥1,

где цифры an , an−1 ,K, ai - верные. По определению, число верных знаков числа

x$ определяется неравенством ∆(x$) = x − x$ <10i . Разделив обе части этого неравенства на x$ , получим

9

δ(x$) =

x

−

x$

≤

10i

x$

an 10n + an−1 10n−1 +K+ai 10i +K+a−m 10−m

.

Это неравенство еще больше усилится, если число x$

заменить заведомо мень-

шим числом an 10n . Тогда будем иметь

δ (xˆ )≤

10i

1

i−n

.

=

10

an 10n

an

Таким образом, относительная погрешность числа характеризуется коли-

чеством верных цифр числа l, отсчитываемых от an

- первой значащей цифры

числа до ai - последняя значащая цифра этого числа ( l = i – n ). За предельную

относительную погрешность δ приближенного числа x$					можно принять вели-
чину	1
δ(x ) ≤		10	i−n	.	(1.6)
$	an

Пример 4. Какова предельная относительная погрешность в процентах чисел x1 = 358,563 и x 2 = 0,047 , если все их цифры верные?

∆ Для числа x1 имеем n = 3, i = −3. По формуле (1.6) получаем

δ1 = 13 10−3−3 100% = 0,00003% .

Для числа x 2 при n = −2 , i = −3 будем иметь

δ1 = 14 10−3−(−2) 100% = 2,5% .

Пример 5. Со сколькими верными десятичными знаками надо взять 3 28 , чтобы погрешность не превышала 0,25%?

∆ Так как n 1 + α ≈1 + αn , при достаточно малых α, то

3	28 = 3		+	1		+	1	≈ 3,037037
		27 +1 = 3 27 1			= 3 1
				27			81	.

По условию δ ≤ 0,25% , т.е. δ ≤ 0,0025. Для числа 3 28 имеем an = 3 , n = 0 . Согласно формуле (1.6), имеем

0,0025 ≥ 13 10i−0 10i ≤ 0,0075 i = −4 .

Таким образом, число 3 28 ≈ 3,037037 имеет пять верных цифр, считая их слева, т.е. 3 28 ≈ 3,0370 .

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

10

Пусть x и y - точные числа, x$ и y$- их приближения, ∆(x$) и ∆(y$) - их абсолютные, а δ (xˆ ) и δ (ˆy) - их относительные погрешности соответственно.

1) Погрешность суммы и разности. Согласно определению (1.1), имеем:
∆(xˆ + yˆ )= (x + y) − (xˆ + yˆ) = (x − xˆ) + ( y − yˆ) ≤ x − xˆ + y − yˆ = ∆(xˆ) + ∆( yˆ) .
Таким образом,	(x	y)	x	(y),	(1.7)
∆	$	+ $	≤ ∆( $) + ∆	$

т.е. абсолютная погрешность суммы не превосходит суммы их абсолютных погрешностей.

Аналогично получаем оценку погрешности для разности двух чисел:

(x y)

(x

y)

(x y)

x

x

(y

y)

x

x

y

y

x

(y),

∆

$ − $

=

−

−

$ −

$

=

(

− $) −

− $

≤

−

$

+

−

$

= ∆( $) + ∆

$

т.е.

(x

y)

x

(y).

(1.8)

Пусть теперь

x1, x 2 ,..., x n

∆

$

− $

≤ ∆( $) + ∆

$

1

, x

2 ,..., x n

- соответствующие

- точные числа, а

x

$

$

$

приближения. Тогда обобщая (1.7), будем иметь

(1.9)

∆(x

1 + x

2 +...+x n )≤ ∆(x1 )+ ∆(x 2 )+...+∆(x n ).

$

$

$

$

$

$

Из полученных результатов вытекает следующее правило сложения приближенных чисел различной абсолютной точности:

1)выделяются числа, имеющие наибольшую абсолютную погрешность (числа наименьшей точности);

2)более точные числа округляются таким образом, чтобы сохранить в них на один знак больше, чем в выделенном числе (запасной знак);

3)производится сложение или вычитание всех чисел с учетом сохраненных знаков;

4)полученный результат округляется на один знак.

Пример 6. Вычислить a = 01732. +17.45 − 0.00354, где все цифры чисел вер-

ные.

∆ Поскольку цифры всех чисел верные, то погрешность вычисления первого слагаемого не превышает 0.0001, второго – 0.01, третьего - 0.00001. Наибольшую погрешность имеет второе слагаемое 17.45, поэтому числа 0.1732 и 0.00354 округляются до трех знаков после запятой. В результате получим приближенное число a$ = 0173. +17.45 − 0.004 =17.619 . Допускаемая при этом погреш-

ность ∆(a) не превышает величины ∆ = 0.01 + 0.001 + 0.00001 = 0.01101. Округлив a

$

$

до двух знаков после запятой, окончательно получим a =17.62 .

$

Относительная погрешность суммы и разности приближенных чисел, со-

гласно (1.3), (1.7) и (1.8), определяются соотношениями:

δ(xˆ + yˆ )=

∆(xˆ + yˆ )

≤

∆(xˆ)+ ∆(yˆ )

,

(1.10)

xˆ + yˆ

xˆ + yˆ

δ (xˆ − yˆ )=

∆

(xˆ − yˆ )

≤

∆(xˆ)+ ∆(yˆ )

,

(1.11)

xˆ − yˆ

xˆ − yˆ

11

соответственно.

2) Погрешность произведения. Согласно (1.2), x = x$ ± ∆(x$), y = y$ ± ∆(y$).Тогда

$$

$$

=

$

( $)

)

$

$

)

$$

=

$$ $

$ $ ( $)

( $)

$ $$

≤

∆(xy)=

xy − xy

(x

+ ∆ x

y +

∆(y)

− xy

xy + x

∆(y)+ y∆ x

+ ∆ x

∆(y)− xy

(

≤

x$

∆(y$)+

y$

∆(x$) + ∆(x$)∆(y$)

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

∆( xˆ ˆy ) ≤

xˆ

∆( ˆy ) +

ˆy

∆( xˆ ) + ∆( xˆ )∆( ˆy ) ,

(1.12)

δ (xˆ ˆy)=

∆(xˆ ˆy)

=

xˆ

∆(ˆy)

+

ˆy

∆(xˆ)+ ∆(xˆ)∆(ˆy)

= δ (xˆ)+δ (ˆy)+δ (xˆ) δ (ˆy).(1.13)

xˆ ˆy

xˆ

ˆy

Из соотношения (1.12) получаем правило умножения приближенных чисел различной абсолютной точности.

1)выделяется число с наименьшим количеством верных значащих цифр;

2)оставшиеся сомножители округляются таким образом, чтобы они содержали на одну значащую цифру больше, чем количество верных значащих цифр в выделенном числе;

3)в произведении сохраняется столько значащих цифр, сколько верных значащих цифр имеет выделенное число.

Пример 7. Найти произведение приближенных чисел x1 = 3.75 и x 2 = 0.73788 , все цифры которых верные.

∆ Число x1 содержит три верные значащие цифры, а число x 2 - пять. Поэтому второе число x 2 округляем до четырех значащих цифр: x 2 = 0.7379 . Тогда

a$ = x1 x 2 = 3.75 0.7379 = 2.767125.

В этом произведении сохраняем три значащие

цифры, так что a = 2.76 . Заметим,

что если считать числа x1 и x 2

точными, то

$

x1 x 2 = 375. 0.73788 = 2.76705.

3) Погрешность частного. При y

≠

0 и y

≠

0 имеем

ˆ

xˆ

x

xˆ

xˆ + ∆( xˆ )

xˆ

xˆ

ˆy

+ ˆy ∆( xˆ ) − xˆ ˆy − xˆ ∆( ˆy )

∆

=

−

=

−

=

=

y

ˆy

ˆy + ∆( ˆy )

ˆy

ˆy( ˆy + ∆( ˆy ))

ˆy

=

ˆy ∆( xˆ ) − xˆ ∆( xˆ )

≤

xˆ

∆( ˆy ) +

xˆ

∆( xˆ )

ˆy2

1+

∆

ˆ

ˆy2

∆

ˆ

( y )

1+

( y )

ˆy

ˆy

Так как

a + b

≥

a

−

b

,

то

1+ ∆( ˆy )

≥

1−

∆( ˆy )

=

1−δ( ˆy )

. Отсюда получаем,

что

ˆy

ˆy

xˆ

∆( ˆy ) +

ˆy

∆( xˆ )

∆

xˆ

≤

,

(1.14)

ˆy

2

1 −δ

( ˆy )

ˆy

12

δxˆ =ˆy

xˆ

∆

ˆ

δ

( xˆ

)

+δ

( ˆy )

y

≤

.

(1.15)

xˆ

1

−δ

ˆ

( y )

ˆy

Таким образом, получаем следующее правило деления приближённых чисел различной абсолютной погрешности:

1)выделяется число с наименьшим количеством верных значащих цифр;

2)оставшиеся числа округляются таким образом, чтобы они содержали на одну значащую цифру больше, чем выделенное число;

3)в полученном результате деления сохраняется столько значащих цифр, сколько верных значащих цифр имеет выделенное число.

На практике при работе с числами достаточно хорошей точности обычно считают, что ∆( xˆ )∆( ˆy ) ≈ 0 , δ( xˆ ) δ( ˆy ) ≈ 0 , 1 −δ( ˆy ) ≈1. Поэтому вместо

формул (1.12)-(1.15) пользуются более упрощёнными формулами:

∆( xˆ ˆy ) ≤ xˆ ∆( ˆy ) + ˆy ∆( xˆ ) , δ( xˆ ˆy ) ≤δ( xˆ ) +δ( ˆy ) ,

xˆ

xˆ

∆( ˆy ) +

ˆy

∆( xˆ )

xˆ

≤δ( xˆ

) +δ( ˆy ).

∆

≤

, δ

ˆy

2

ˆy

ˆy

1.3. Вычислительная погрешность

Представление чисел в ЭВМ. Одним из источников вычислительной погрешности являются особенности представления чисел в ЭВМ и арифметические операции с ними. Рассмотрим эти особенности.

Целые числа, представимые в ЭВМ (целые машинные числа) - это все целые числа, принадлежащие замкнутому интервалу [N¯,N+]. В языке Паскаль, например - это числа диапазона [–32768, 32767]. Арифметические операции с целыми машинными числами выполняются на ЭВМ точно, если результат принадлежит интервалу [N¯,N+]. В противном случае фиксируется ошибочная ситуация, а результат операции не определен или не имеет смысла. Поэтому при вычислениях следует заботиться о том, чтобы результаты арифметических операций не выходили за пределы интервала [–32768, 32767].

Вещественные числа могут быть записаны в двух формах. Обычную форму записи числа в виде

x = ±an an −1 La1a0 , a−1a−2 La−m

называют записью с фиксированной запятой. В позиционной системе счисле-

ния с основанием β эта запись означает, что

x = ±an β n + an−1 β n−1 +K+ a1 β + a0 + +a−1 β −1 + a−2 β −2 +K+ a−m β −m .

Каждое из чисел an может принимать одно из значений {0, 1, …, β -1}. Числа an называются разрядами. Например, в десятичной системе запись 1358,7604 определяет число

1 103 + 3 102 + 5 10 + 8 + 7 10−1 + 6 10−2 + 0 10−3 + 4 10−4 .

13

β−1 ≤ m <1.

Третий разряд перед запятой этого числа равен 3, второй разряд после запятой – 6.

О числах, записанных в виде 0,63750 106 ; 637,50 103 ; 6,3750 105 ; 6375,0 10 ,

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

x =m β p

(1.16)

Здесь p - целое число (положительное, отрицательное или нуль), называемое порядком числа x, m - число с фиксированной запятой, называемое мантиссой числа x. Первая после запятой цифра числа m всегда отлична от нуля:

(1.17)

Так, в рассмотренном примере, нормализованной формой записи числа 63750 будет 0,63750 106 . Заметим, что в этой записи все цифры мантиссы верные. Для числа - 0,00384 нормализованной будет форма − 0,384 10−2 .

Для записи вещественных чисел в ЭВМ отводится фиксированное число разрядов (разрядная сетка), в которой выделены разряды для записи мантиссы, порядка, знаков мантиссы и порядка. Таким образом, множество машинных вещественных чисел характеризуется следующими четырьмя целыми константами: основанием счисления β, β ≥2, точностью t, нижней границей экспоненты e¯ и верхней границей экспоненты e+. Ненулевые машинные вещественные числа имеют вид

d2

p

= ±.d1d2 ...dt β

p

d1

dt

x = ±

+

+ ... +

β

β

β

2

β

t

где 0 ≤ di < β, di ≠ 0, e− ≤ p ≤ e+ .

Вещественный машинный нуль нуль представляется в ЭВМ специальным образом. Обычно это число, у которого d1 = d2 =... = dt = 0 или p < e− .

В языке Паскаль реализованы следующие типы вещественных чисел: real, single, double, extended, некоторые характеристики которых приведены в табли-

це 1.1.

Минимальное положительное число ε0 , которое может быть представлено в ЭВМ, называется машинным нулем (ε0 = βe−−1 ). Максимальное положитель-

ное число ε∞ - машинной бесконечностью (ε∞ = βe+ (1− β−t )). Таким образом,

машинные вещественные числа	принадлежат замкнутому интервалу
−ε∞ ≤ x ≤ε∞. На интервале −ε0 ≤ x ≤ε0	лежит только одно машинное веществен-
ное число – машинный нуль.
	Таблица 1.1

14

Тип	Наименова-	Занимаемая	Число зна-	Диапазон допустимых
	ние	память	чащих цифр	значений
real	веществен-	6 байт	11-12	2.9 10 −39	÷1.7 10 38
	ный
single	с одинарной	4 байта	7-8	1.5 10 −38	÷ 3.4 10 38
	точностью
double	с двойной	8 байт	15-16	5.0 10 −324	÷1.7 10 308
	точностью
extended	с повышен-	10 байт	19-20	3.4 10 −4932	÷1.1 10 4932
	ной точно-
	стью

Округление чисел в ЭВМ. Из-за конечности разрядной сетки в ЭВМ можно представить точно лишь конечное подмножество вещественных чисел. Рассмотрим два машинных вещественных числа f1 и f2 (f1 < f2) таких, что между f1 и f2 нет других машинных чисел. Справедливо неравенство

f2 − f1	≤ε	1	= β1−t .
f1

Число ε1 = β−t+1 называют машинным эпсилоном. Оно зависит от разряд-

ности мантиссы и характеризует точность представления чисел в ЭВМ. Для любого вещественного числа x, f1 ≤ x ≤ f2 выполняется неравенство

max( f1 −x, f2 −x )≤ε1 x .

Таким образом, любое вещественное число x, не представимое точно, заменяется ближайшим к нему числом x$, представимым точно. Такая замена называется округлением числа x. Величина погрешности замены не превышает ε1 x и зависит от способа округления.

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

δ( xˆ )≤ε1 .

(1.18)

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

δ( xˆ ) ≤	1	β−t+1	=	1	ε .	(1.19)
	2			2	1

Таким образом, можно считать, что машинные константы ε0 , ε∞ и ε1 характеризуют погрешность представления вещественных чисел в ЭВМ.

15

Можно считать, что точное число x из диапазона ε0 ≤ x ≤ε∞и отвечающее ему округленное число x$ связаны равенством

xˆ = x(1+ε ),		ε		≤ε1.	(1.20)

Если же число превышает по абсолютной величине ε∞, то оно не приближается

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

ются машинным нулем с погрешностью ε0 .

Накопление погрешностей округления. Ошибочно полагать, что в силу малости величин ε0 и ε1 , ошибки округления практически не влияют на резуль-

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

Пусть a и b - два машинных числа, а знак обозначает одну из четырех арифметических операций. Поскольку ЭВМ выполняет промежуточные вычисления с двойной точностью, т.е. с мантиссой, содержащей 2t разрядов, то округлению до t разрядов подвергается лишь окончательный результат. Таким образом, выполнение каждой арифметической операции вносит относительную погрешность, не превышающую ε1 :

								δ (a b)≤ε1,		(1.21)
или									+ε) +ϕ,
								(a b)=(a b) (1
где		ε		≤ε1,		ϕ		≤ε0 , εϕ =0.

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

Для оценки влияния погрешностей округления на результат вычислений по некоторому алгоритму (прямая задача теории погрешностей), предполо-

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

y = f (x1, x2 ,...,xn ) = f (x).

Поскольку входные данные xˆ1,xˆ2 ,...,ˆxn неточны,	ˆr
	ˆy = f ( xˆ1,xˆ2 ,...,xˆn ) = f ( x ) - прибли-

женное значение результата. Пусть известна область G в пространстве пере-

менных x1,x2 ,...,xn :

G ={xr xi −xˆi ≤∆xi ,i =1,...,n},

которой принадлежат значения входных параметров. Если f (x1, x2 ,...,xn ) - непре-

рывно дифференцируемая функция своих аргументов, то по формуле конечных приращений Лагранжа имеем следующую оценку погрешности:

16

r

r

n

∂f (θ )

n

∆

=

y

−

ˆy

=

f (x)

−

ˆ

=

∑

∆

≤

∑mi

∆

.

(1.22)

y

f (x)

∂xi

xi

xi

i−1

i−1

) и ( xˆ ,xˆ

,...,xˆ

),

Здесь (θ

,θ

2

,...,θ

n

) - точка, расположенная между ( x ,x

2

,...,x

n

2

n

1

r

1

1

mi = sup

∂f (x)

∂xi

G

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

Если f (x , x ,...,x ) ≠ 0		и			x ≠ 0,i =1,...,n, справедлива следующая оценка от-
1 2	n				i
носительной погрешности:					xir			∂f (θ )
		n
	δ( ˆy ) = ∑									δ( x ).	(1.23)
		i −1 f ( x )								i
Коэффициенты								∂xi
					xir		∂f (θ )
	Ki	=
								∂xi
				f (x)

являются коэффициентами усиления относительных погрешностей входных данных. Если значения всех коэффициентов усиления ограничены по абсолютной величине единицей, то накопления погрешностей не происходит. Коэффициенты Ki часто называют числами обусловленности задачи.

Часто возникает обратная задача теории погрешностей: с какой точно-

стью должны быть заданы входные данные, чтобы получить результат требуемой точности?

Применительно к поставленной выше прямой задаче оценивания погрешности результата вычисления значения функции u = f (x1, x2,..., xn)при заданных оценках погрешностей аргументов, обратная задача заключается в оценивании ∆(xi) (или δ(xi)) по известной величине приращения∆y . Поскольку главной, т.е. линейной частью этого приращения является полный дифференциал dy , за границу абсолютной погрешности результата приближенно может быть принята величина

	n	∂u		(1.24)
	∆(u) = ∑ |		| ∆(xi)
		∂xi
	i=1

Из формулы (1.23) легко получается формула приближенной оценки относительной погрешности значения u :

	∆	n	∂u		∆	i	n	∂u	n	∂ln u
δ(u) =	u = ∑ |			|		(x ) = ∑ |			| ∆(xi) = ∑ |		| ∆(xi) (1.25)
	|u|	i=1	∂xi		|u|		i=1 u∂xi		i=1	∂xi

17

Для случая дифференцируемой функции одной переменной грубое решение обратной задачи очевидно:

если y = f (x) , то ∆(u) ≈| dy |=| f ′(x)| ∆x , откуда ∆x ≈ ∆(u) | f ′(x)|

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

Например, можно применить принцип равных влияний, состоящий в

предположении, что частные дифференциалы | ∂u | ∆(xi) в (1.24) одинаково

∂xi

влияют на погрешность значения функции, тогда ∆(u) =n | ∂u | ∆(xi) , откуда

∂xi

∆(xi) = ∆(u) . n |∂u |

∂xi

В качестве другого довольно естественного допущения можно принять равен-

ство

относительных погрешностей

всех

аргументов, т.е. считать

∆(xi)

при всех i=1,n.

δ(xi) =

|xi| = p

n

Тогда ∆ i

=

i

, и значит

∂u

.

(x )

p

| x |

∆(u) = p ∑ |

xi|

∂xi

i=1

, характеризующую

Из последнего равенства получаем величину

∆(u)

p

=

n

∂u

∑ |xi

|

∂xi

i=1

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

границы абсолютных погрешностей принимаем	∆(xi) =	|xi| ∆(u) .
		n	∂u
		∑ |xi		|
			∂xi
		i=1

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

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

В частности, если n>10 и все слагаемые округлены до m-го десятичного разряда, то для подсчета абсолютной погрешности суммы S применяют правило Чеботарева

∆(S) = 3n 0,5 10−m

18

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

Пусть x=x1+x2+...+xn - среднее арифметическое n (n>10) приближенных n

чисел (например, результатов измерений), имеющих одинаковый уровень абсолютных погрешностей ∆(xi) =0,5 10−m . Тогда классическая оценка абсолютной

погрешности величины x есть

∆(x) =1n (∆(x1)+∆(x2)+...+∆(xn)) = 1n n 0,5 10−m =0,5 10−m =∆xi

т.е. такая же, как у исходных данных. В то же время по формуле Чеботарева

∆(x) = 1	3n 0,5 10−m =	3	0,5 10−m =	3	∆(xi) →0.
n		n		n	n→∞

Как видим, применение правила Чеботарева приводит к естественному выводу о том, что арифметическое усреднение результатов измерений увеличивает точность, чего нельзя сказать на основе классической теории погрешностей. Прямое применение статических оценок погрешностей также является достаточно сложным делом. Однако именно такие оценки подкрепляют практические правила работы с приближенными числами, составляют основу технического подхода. Этот подход связывают с именем известного кораблестрои-

теля, математика и механика А.Н. Крылова (1863-1945г).

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

Устойчивость и корректность постановки задачи. Пусть в результате решения задачи по исходному значению величины х находится значение искомой величины у. Если исходная величина имеет абсолютную погрешность ∆(x), то решение у имеет погрешность ∆(y). Задача называется устойчивой по исходному параметру х, если решение у непрерывно зависит от х, т.е. если малое приращение исходной величины х приводит к малому приращению искомой величины у. Другими словами, малые погрешности в исходной величине приводят к малым погрешностям в результате расчётов.

Задача называется поставленной корректно, если для любых значений исходных данных её решение существует, единственно и устойчиво по исходным данным.

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

Пример 8. Пусть задана система

19

3,0000 x −7,0001 y = 0,9998,3,0000 x −7,0000 y =1.

Легко проверить, что х=5, у=2 являются решениями этой системы. Рассмотрим теперь систему, близкую к ней:

3,0000 x − 7,0001 y =1,3,0000 x − 7,0000 y =1.

Эта система, как нетрудно убедиться, имеет решение х= 13 , у=0.

Как видно, изменение исходных данных на 2 10−4 одной из правых частей значительно изменило решение системы.

Рассмотрим снова систему, близкую к исходной:

3,0000 x −7,0000 y = 0,9998,3,0000 x −7,0000 y =1.

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

часто встречаются в приложениях при решении практических задач.

20

ГЛАВА 2. ЧИСЛЕННЫЕ МЕТОДЫ ЛИНЕЙНОЙ АЛГЕБРЫ 2.1. Основные сведения из линейной алгебры

К численным методам линейной алгебры относятся численные методы решения следующих задач:

-решение систем линейных алгебраических уравнений,

-обращение матриц,

-вычисление определителей квадратных матриц,

-нахождение собственных значений и векторов матриц (решение спектральной задачи линейной алгебры).

Напомним основную терминологию - определения, обозначения, понятия

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

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

a	a	...	a
11	12		1n
A = a21	a22	...	a2n
	...	...	...
...
am1	am2	...	amn

называют прямоугольной матрицей размера m×n или m×n - матрицей. Числа aij

называют элементами этой матрицы. Первый индекс (i) - номер строки, а второй индекс (j) - номер столбца, на пересечении которых расположен данный элемент aij . Матрица, состоящая из одной строки, называется вектор-строкой и

обозначается xr = (x1, x2 ,..., xn ). Матрица, состоящая из одного столбца - век-

тор-столбцом. Матрица А размера п×n называется квадратной, а число n - порядком этой матрицы.

Квадратная матрица вида

α1		0	...	0
	0	α2	...	0
D =		...	...
...				...
	0	0	...
				αn

называется диагональной. Если все элементы диагонали равны между собой, матрица называется скалярной, и если они равны единице – единичной (обозначается E или I).

Квадратная матрица U называется верхней треугольной, если ее элементы, расположенные выше диагонали, не равны нулю:

21

a11		a12	...	a1n
	0	a22	...	a2n	.
U =		...	...	...
...
	0	0	...
				ann

Квадратная матрица L называется нижней треугольной, если ее элементы, расположенные ниже диагонали, не равны нулю:

			a11		0			...	0
			a		a			...	0
			L =	21		22		...	...
			...		...
					an2			...
			an1						ann
Переставив	в матрице строки				и столбцы, получают так называемую
транспонированную к матрице А матрицу размера n×m,
			a11		a21			...	am1
			a		a		22	...	a	m2
			AT = 12					...
			...		...				...
					a2n			...
			a1n						amn
с элементами a/ = a		ji	.
ij
Квадратная	матрица называется				симметричной, если она равна своей

транспонированной матрице: aij = a ji .

Квадратная матрица называется ленточной, если ее элементы удовлетворяют соотношениям: aij = 0, m1 < i − j; j −i > m2 для некоторых чисел m1, m2.

Шириной ленты называют число s= m1+ m2+1.

Матрица называется разреженной, если отношение количества ненулевых элементов к общему числу ее элементов (коэффициент плотности) меньше или равно 1/n.

Сравнивать, складывать и вычитать можно только матрицы одинаковых размеров. Указанные операции выполняются над одноименными элементами матриц.

Произведение матрицы Am×n справа на матрицу Bl×k определено только в том случае, когда число столбцов матрицы А равно числу строк матрицы B( l = n ) . Произведением С=АВ называется матрица Cm× k с элементами

n

cij = ∑aipbpj , i =1,...,m; j −1,...,k.

p=1

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

22

матрице ВА (даже в случае квадратных матриц) т.е. операция умножения матриц не коммутативна.

Матрица A называется нормальной, если она перестановочна со своей транспонированной матрицей: ATA = AAT.

Определитель квадратной матрицы. Миноры. Ранг матрицы. Опре-

делителем квадратной матрицы порядка п называется сумма n! членов вида (−1)k a1k1 a2k2 ...ankn , каждый из которых равен произведению п элементов мат-

рицы, взятых по одному из каждой строки и каждого столбца. При этом знак члена определяется количеством инверсий k в перестановке k1, k2, …, kn из чисел

1,2, …,n .

Перестановкой из п чисел 1, 2, 3, ..., п называется любое упорядоченное расположение этих чисел. Количество всех перестановок равно n!. Два числа в перестановке образуют инверсию, если большее число стоит перед меньшим. Например, в перестановке 2413 число инверсий равно 3.

Определители второго и третьего порядков, например, равны:

det A =

a11

a12

= (−1)0 a

a

22

+(−1)1 a a

21

= a

a

22

−a a

,

a21

a22

11

12

11

12

21

a11

a12

a13

det A =

a21

a22

a23

= a11a22a33 +a12a23a31 +a13a21a31 −a13a22a31 −a11a23a32 −a12a21a33.

a31 a32 a33

Вычисление определителей более высокого порядка можно свести к вы-

числению определителей второго и третьего порядков.

Легко заметить, что определитель треугольной и диагональной матриц ра-

вен произведению диагональных элементов.

Минором элемента aij квадратной матрицы A называется определитель

матрицы Mij порядка n-1, полученной вычеркиванием

i-й строки и

j-го столбца

матрицы. Алгебраическим дополнением элемента

называется

выражение

Aij

=

−

i+ j

det Mij

.

( 1)

Минором k–го порядка матрицы A называется определитель матрицы M, образованной элементами, стоящими на пересечении любых k строк и столбцов матрицы A. Минор матрицы называется главным, если он образован элементами, стоящими на пересечении строк и столбцов с одинаковыми номерами. Рангом матрицы A называется максимальный порядок отличных от нуля миноров матрицы. Ранг прямоугольной матрицы Am×n не может превышать min{m,n} и

равен r, если найдется, по крайней мере, один ее минор порядка r, отличный от нуля, а все миноры порядка r+1 и выше равны нулю.

Обратная матрица. Квадратная матрица называется невырожденной, если ее определитель отличен от нуля, в противном случае матрица называется

23

вырожденной. Для невырожденной матрицы существует обратная матрица A-1 такая, что A-1 A = A A-1 = E.

Для матриц небольшого порядка обратную матрицу можно найти, использую матрицу, составленную из алгебраических дополнений (присоединенную матрицу):

		A	A	...	A	T
		11	12		1n
A−1 =	1	A21	A22	...	A2n .		(2.1)
				...
	det A ... ...				...
		A	A	...	A
		n1	n2		nn

Матрица, обратная произведению С=АВ квадратных невырожденных матриц, определяется следующим образом: С-1 = В-1А-1.

В численных методах линейной алгебры большую роль играют ортогональные матрицы. Квадратная матрица является ортогональной, если ATA = E, т.е. для ортогональной матрицы AT = A-1. Ортогональные матрицы обладают следующими свойствами:

-строки (столбцы ) ортогональных матриц попарно ортогональны;

-сумма квадратов элементов каждой строки (столбца) ортогональной матрицы равна 1;

-определитель ортогональной матрицы равен 1.

Норма вектора и норма матрицы. Множество Rn всех n -мерных веще-

ственных векторов образует линейное пространство E . Нормой вектора x Rn называется действительное число xr , которое ставится в соответствие вектору

и обладает свойствами (аксиомами нормы):

1. xr ≥ 0 , причём xr = 0 тогда и только тогда, когда xr- нулевой вектор

( x = (0,0,...,0);

2.αxr = α xr для любого действительного числа α ;

3. xr + yr ≤ xr + yr для любых x и y из E .

Свойство 3 нормы вектора называется неравенством треугольника. Пространство E с определённой на нём нормой называется линейным

нормированным пространством.

Норма в пространстве E вводится неоднозначно несколькими способами: 1. Норма x1 определяется соотношением

n

xr

1 = ∑

xi

=

x1

+

x2

+... +

xn

(2.2)

xr

j=1

2. Норма

2 имеет вид

xr

=

x2

+ x2

+... + x2

(2.3)

2

1

2

n

24

и определяется с помощью скалярного произведения векторов xr = (x1 , x2 ,..., xn ) и

yr = (y1 , y2 ,..., yn )

из E = R n :

(x, y)=x1 y1 + x2 y2

+ ... + xn yn .

(2.4)

3. Норма

xr

3

вводится с помощью соотношения

xr

∞

= max{

x

,

x

2

,...

x

n

,}= max

x

j

..

(2.5)

1

1≤ j ≤n

Легко проверить, что величины

xr

1 ,

xr

2 ,

xr

3 , определённые равенствами

(2.2), (2.3), (2.5) соответственно, удовлетворяют аксиомам 1-3 нормы вектора. Эти наиболее распространенные на практике нормы называют соответственно октаэдрической, сферической и кубической по названию поверхностей, определяемых уравнением xr =1, xr R3 .

Пространства E с различными введенными нормами, вообще говоря, являются различными линейными нормированными пространствами. Вектор x из E называется нормированным, если его норма равна единице, т.е. xr =1.

Пусть А - квадратная n × n матрица. Тогда для любого xr E вектор y = Ax тоже принадлежит E , и его норма равна yr = Axr .

Норму

A

матрицы А определим следующим образом

A

= supr

Ar

xr

,

(2.6)

xr

Axr

xr≠0

x

`

где

и

- введённые выше нормы

векторов x и Axr . Норма матрицы, оп-

ределённой равенством (2.6), называется согласованной с нормой вектора. Из определения (2.6) имеем

Axr

Axr

r

≠ 0 .

A

= sup

r

≥

r

,

x

xr

≠0

x

x

Отсюда

Axr

A

xr

(2.7)

≤

причём это соотношение справедливо и при xr = 0 .

На основании (2.7) для любого x из Е выполняются неравенства

|| BAxr ||≤|| B(Ax) ||≤|| B || || Ax ||≤|| B || || A || || xv || ,

отсюда

|| A2 x ||≤|| A ||2 || xr||,...,|| Ak xr||≤|| Ak || || xr||,

т.е.

|| Ak ||≤|| A ||k .

Вид согласованной нормы матрицы зависит от избранной нормы вектора в Е. Так, для x1 имеем

25

n

n

n n

n

n

n

n

r

1

∑

∑

ij

j

∑∑

ij

j

∑

j

∑

ij

ij

∑

j

Ax

=

a

x

≤

a

x

=

x

a

1≤ j≤n ∑

a

x

.

≤ max

i=1

j=1

i=1 j=1

j=1

i=1

i=1

j=1

Можно найти вектор x из E такой, что это неравенство выполняется как верное равенство, т.е. будем иметь для найденного x :

		n
r			a				r
Ax		= max					x
	1	1≤ j≤n ∑		ij						1
		i=1

Отсюда и из определения (2.6) получаем, что если норма вектора определяется равенством (2.2), то согласованная норма матрицы А даётся формулой

n

,

(2.8)

A

a

= max

1

1≤ j≤n ∑

ij

т.е. норма

i=1

A

1 есть максимальная сумма

модулей элементов

столбцов

матрицы А.

Для нормы xr 2 , определяемой скалярным произведением, согласованная норма A 2 матрицы А имеет вид

1

n

n

2

A

2

=

∑∑aij2

.

(2.9)

Для нормы

xr

i=1 j=1

3 согласованная норма матрицы А даётся формулой

n

A

= max

a

,

(2.10)

3

1≤i≤n ∑

ij

j =1

т.е. норма A 3 есть максимальная сумма модулей элементов строк матрицы А.

Системы линейных алгебраических уравнений (СЛАУ). Система m

линейных алгебраических уравнений с n неизвестными записывается в виде:

a x		+ a x			+K+ a x			= b ,
11 1		12		2		1n n			1
a21x1 + a22 x2 +K+ a2n xn = b2 ,										(2.11)
LLLLLLLLLLLL
a	x	+ a	m2	x		+K+ a	x		= b .
	m1 1			2			mn n		m

Здесь aij – заданные числа – коэффициенты системы; b1, b2 ,..., bm – тоже заданные числа – правые части системы; x1 , x2 ,..., xn – неизвестные, которые

следует найти. Введя матрицу A коэффициентов и векторы-столбцы неизвестных и правых частей (свободных членов), систему (2.11) можно коротко запи-

сать в векторно-матричном виде:

Решением	системы	Axr = b.			(2.12)
		(2.11)	называется	набор	чисел

x1 =α1 , x2 =α2 ,..., xn =αn , подстановка которых в каждое из равенств (2.11) превращают всех их в тождества. Система (2.11) называется совместной, если

26

она имеет хотя бы одно решение. Если система не имеет решений, она называ-

ется несовместной (противоречивой). Система (2.11) называется определенной,

если она имеет единственной решение, и неопределенной, если решений больше одного.

Расширенной матрицей системы (2.11) называется матрица B размерами n ×( n +1), полученная из матрицы A добавлением столбца свободных членов и содержащая всю информацию о системе:

a11

a12

....

a1n

b1

r

a

a

....

a

b

22

2n

2

B = [A

b ]=

11

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

....

.

an2

....

ann

bn

an1

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

Вопрос о совместности системы линейных алгебраических уравнений (СЛАУ) полностью решается теоремой Кронекера-Капелли, которая утвержда-

ет, что система линейных алгебраических уравнений (2.11) совместна тогда и только тогда, когда ранг расширенной системы В равен рангу матрицы A : r(А)= r(B). Совместная система имеет единственное решение, если ранг системы равен числу неизвестных n, и бесконечно много решений, если r(A)<n.

Если матрица A - квадратная и ее определитель A ≠ 0 , то такая система

всегда совместна и имеет единственной решение, которое может быть найдено по формулам Крамера,

	xi =	det Ai	, i =1,..., n	(2.13)
		det A

где det Ai - определитель матрицы, в которой столбец коэффициентов при вычисляемом неизвестном заменен столбцов свободных членов. Эти формулы применимы на практике только при небольших n (n = 2,3).

При решении совместной системы линейных алгебраических уравнений с прямоугольной матрицей вначале устанавливают ранг матрицы - r(A)= r, а затем переставляют уравнения и неизвестные таким образом, чтобы отличный от нуля минор матрицы, имеющий порядок r, занял положение в верхнем левом углу матрицы. Если r=n и r≤m, система имеет единственное решение, которое можно найти, решая первые n уравнений преобразованной системы. Если r<n и r≤m, то систему решают относительно первых r неизвестных x1 , x2 ,..., xr , а неиз-

вестным xr+1 , xr+2 ,..., xn придают произвольные значения. Различных наборов

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

Таким образом, при решении СЛАУ общего вида необходимо рассмотреть следующее вопросы:

- совместна или несовместна система;

27

-определена или не определена совместная система;

-как описать совокупность решений неопределенной системы;

-как разработать метод получения единственного решения определенной системы.

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

Методы решения СЛАУ делятся на две группы – прямые и итерационные. Прямые методы используют конечные соотношения или формулы для вычисления неизвестных. К прямым (или точным) методам решения СЛАУ относятся алгоритмы, которые в предположении, что вычисления ведутся без округлений, позволяют получить точное решение системы за конечное число арифметических действий Q(A). В прямых методах число необходимых для решения задачи арифметических операций зависит только от вида вычислительной схемы и от

порядка матрицы n. Считают, что для произвольной невырожденной матрицы прямые методы с оценкой числа операций Q(A)=O(n3) являются оптимальными.

Итерационные методы – это методы последовательного улучшения решения и методы последовательных приближений. В них задается начальное приближенное решение. Потом с помощью некоторого алгоритма проводится один цикл вычислений, уточняющих начальное приближенное решение. Если уточненное значение не удовлетворяет требуемой точности, то, приняв это новое значение за начальное, снова с помощью алгоритма улучшают его. Итерации проводятся до тех пор, пока не будет получено решение с требуемой точностью. Алгоритмы решения линейных систем с использованием итерационных методов обычно более сложные по сравнению с прямыми методами и объем вычислений заранее трудно определить.

2.2. Прямые методы решения систем линейных алгебраических уравнений (СЛАУ)

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

В основе рассматриваемых ниже методов лежит следующее утверждение

Пусть все главные миноры матрицы A отличны от нуля, т.е.

a11 ≠ 0,	a11	a12	≠ 0, ..., det A ≠ 0.	(2.14)
	a21	a22

Тогда существуют такие нижняя L и верхняя U треугольные матрицы, что матрица A представима в виде произведения этих матриц:

A = LU.

(2.15)

Если элементы диагонали одной из треугольных матриц фиксированы (ненулевые), то такое разложение единственно.

28

Например,

a11	a12		....	a1n		l11
a	a	22	....	a		l
11			....		2n	= 11
.... ....				....		....
	an2		....
an1				ann		ln1

	0	....	0		1	U12	....	U1n
l	22	....	0		0	1	....	U
						....	....		2n
.... ....			.... ....					....
ln2		....			0	0	....	1
			lnn

Для решения исходной задачи не обязательно искать в явном виде именно разложение (2.15). Принципиальное значение этого утверждения заключается в обосновании возможности сведения общей задачи решения системы (2.12) к двум специальным:

Ly =b,

(2.16)

Ux = y.

Таким образом, решение системы прямыми методами может быть осуществлено поэтапно: на первом этапе систему (2.12) преобразуют к одной из систем (2.16), на втором - решают упрощенную систему и получают значения неизвестных.

К числу наиболее известных прямых методов относится метод Гаусса и его модификации. Рассмотрим вначале вариант, который иногда называют схе-

мой единственного деления.

Метод Гаусса (схема единственного деления). Основная идея метода состоит в том, что система (2.11) приводится к эквивалентной ей системе с верхней треугольной матрицей и единичной диагональю (прямой ход исключения). Из полученной системы последовательно, начиная с последнего уравнения, находятся неизвестные xn , xn −1,..., x1 (обратный ход исключения).

Рассмотрим применение метода на примере системы третьего порядка:

a			x +a			x +a			x =b ,
	11		1	12		2	13			3	1
a21x1+a22 x2+a23 x3=b2 ,												(2.17)
a		31	x +a		32	x +a		33		x =b .
			1			2				3	3

Предположим, что a11 ≠ 0 , разделим первое уравнение системы на коэффициент a11 . Для исключения x1 из второго уравнения вычтем из него первое уравнение, умноженное на a21 . Затем, умножив первое уравнение на a31 , вычтем его из третьему уравнению, исключая тем самым из него x1 . Получим равно-

	сильную систему уравнений вида
	′	′	′
	x1+a12 x2+a13 x3=b1,
	′	′	′	,		(2.18)
	a22 x2+a23 x3=b2
	′	′	′	,
	a32 x2+a33 x3=b3
	где
	a1′j = a1 j / a11 , j = 1,2,3,	b1′ = b1 / a11 ,				(2.19)
	aij′ = aij −ai1a1′j ,	bi′	= bi −ai1b1′,		i, j = 2,3

29

Аналогично исключим x2 из третьего уравнения системы (2.17). Пусть

′	≠ 0 ,								′	. Для исклю-
a22		разделим второе уравнение системы на коэффициент a22
чения x2		из третьего уравнения вычтем из него первое уравнение, умноженное
	′
на a32 . В результате получим систему
				′		′		′
			x1 + a12 x2			+ a13 x3	= b1,
						′′		′′		(2.20)
					x2 + a23 x3		= b2 ,
						′′		′′
где						a33 x3	= b3 ,
		′′	′	′	,,	′′	′	′	,
		a23	= a23 / a22			b2 = b2		/ a22
		′′	′	′	′′	′′	′	′	′′
		a33 = a33		− a32a23		, b3 = b3 − a32b2.

Обратный ход начинается с решения третьего уравнения системы (2.20). Находим

x3 =	b3′′	.
	′′
	a33

Используя найденноеx3 , из второго уравнения системы (2.20) найдем x2 :

x2 = b2′′ − a23′′ x3 ,

затем – x1 :

x1 = b1′ − a12′ x2 − a13′ x3.

Обобщим этот алгоритм на случай линейной системы с произвольным числом уравнений. Для удобства дальнейших преобразований обозначим эле-

менты правой части в расширенной матрице как ai(n+1)

a

a

a

....

a

a

11

12

13

1n

1(n+1)

a21

a22

a23 ....

a2n

a2(n+1)

(2.21)

B =

.

... ...

... .... ...

...

an1

an2

an3 ....

ann

an(n+1)

На первом шаге разделим первое уравнение системы на a11 . Обозначим

коэффициенты

полученного

приведенного

уравнения

U (1) = a

/ a ,

( j = 2,..., n +1) , домножим его на коэффициент a

21

и вычтем

1 j

1 j

11

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

из второго уравнения

(обнуляя

коэффициент a21 матрицы).

Поступим

аналогично с

остальными

уравнениями и получим новую систему, матрица которой в первом столбце, кроме первого элемента, содержит только нули, т.е.

30

	1	U12(1)	U13(1)	....	U1(1)n		U1((1)n+1)
	0	a22(1)	a23(1)	.... a2(1)n			a2((1)n+1)
... ...			... .... ...				...	.	(2.22)
	0	a(1)	a(1)	....	a(1)		a(1)
		n2	n3		nn
							n(n+1)

Первое уравнение в дальнейших преобразования не участвует. Описанный выше процесс исключения неизвестных применим к матрице [aij(1) ] разме-

рами (n-1) ×n, которую называют активной. После k аналогичных шагов получим k приведенных уравнений с коэффициентами

Ulj(l) = alj(l−1) / all(l−1) , (l =1,...,k, j = k +1,..., n +1)

(2.23)

и активную матрицу [aij(k ) ] размерами (n - k) × (n - k+1), элементы которой вы-

числяются по формулам

aij(k ) = aij(k −1) − aik(k −1) ×U kj(k )

, (i = k +1,..., n, j = k +1,..., n +1) .

(2.24)

Прямой ход метода Гаусса заканчивается после n шагов определением Un(n(n)+1) .

Обратный ход метода Гаусса заключается в последовательном определении компонент решения, начиная с хn и заканчивая x1 , по следующим формулам:

n
xn =Un(n(n)+1) , xk =Uk(k(n)+1) − ∑Uki(k ) xi , k = n −1,...,1.	(2.25)

i=k+1

Элементы akk(k −1) , на которые осуществляется деление, называются ведущими или главными элементами. Для возможности выполнения расчетов по описанному алгоритму предполагалось, что все они отличны от нуля ( akk(k −1) ≠ 0 ).

			Решение системы алгебраических уравнений методом							Гаусса требует
	1	(n3	+ 3n2 + n) операций умножения и деления и					1	(2n3 + 3n2	−5n) операций
3							6
сложения. Таким образом, оценка					числа необходимых для решения задачи
арифметических операций Q( A) ≈				2	n	3	и метод Гаусса относится к числу опти-
				3

мальных. Для реализации метода требуется хранить n2 + n +1 машинных слов.

Метод Гаусса с выбором ведущего элемента. Возможность прекраще-

ния вычислений из-за деления на нуль - не единственная неприятность, связанная с реализацией алгоритма метода Гаусса по формулам (2.23) – (2.25). Так как реальные вычисления производятся с машинными числами, то неизбежны ошибки округления и их накопление. Если ведущий элемент очень мал, то при

31

делении и дальнейшем вычитании приведенного уравнения из последующих возникает большая вычислительная погрешность. Для устранения обрыва вычислительного процесса и уменьшения этой погрешности рекомендуется осуществлять перестановку уравнений. Если на каком-либо шаге m ведущий элемент равен нулю (или меньше некоторого малого числа ε > 0 близкого к минимально допустимому на ЭВМ числу), то нужно переставить строки расширенной матрицы местами, сделав ведущей ту строку, у которой коэффициент при xm не равен нулю. Это всегда возможно, т.к. по предположению система одно-

значно разрешима и все коэффициенты при xm не могут одновременно рав-

няться нулю.

Метод с выбором ведущего (разрешающего или главного) элемента является одной из модификаций метода Гаусса, позволяющей исключить деление на нуль и повысить точность вычислений. Идея метода заключается в том, чтобы при прямом ходе в алгоритме метода Гаусса на каждом шаге исключения производить выбор наибольшего по модулю элемента в качестве ведущего. Этого достигают перестановкой строк или столбцов матрицы коэффициентов. Существуют различные стратегии выбора ведущего элемента: выбор в строке, столбце или по всей активной матрице. Наибольший эффект достигается при выборе наибольшего по модулю элемента во всей активной матрице, однако такая модификация метода Гаусса встречается довольно редко, также как и модификация с выбором элемента в строке. Это связано с сильным усложнением алгоритма, т.к. перестановка столбцов местами равносильна переобозначению переменных. Возникает необходимость запоминать номера переставляемых столбцов и выполнять обратное переобозначение переменных после окончания вычислений. Наиболее распространённой в вычислительной практике является стратегия выбора главного элемента столбца - нахождение максимального по

модулю элемента k-го столбца матрицы [aij(k ) ] и использование его в качестве

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

Эффективность использования максимального ведущего элемента во многих случаях зависит от величины элементов строк и столбцов матрицы А. Желательно, чтобы матрица была уравновешена, т.е. выполнялись условия

n					n
ri ≈ ci , где ri = ∑		aij		,	ck = ∑		a jk		.
j =1					j =1

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

32

симального элемента в строке – неявное масштабирование, при котором в каче-

стве ведущего выбирают элемент a(jjj−1) , для которого			a(jjj−1)
			rj

				~( j−1)
				aij
=	max
				ri
	i= j,...,n

⎟.

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

Специфической особенностью решения СЛАУ методом Гаусса является то, что в процессе приведении матрицы системы к треугольному виду преобразованиями прямого хода правая часть системы участвует только пассивно – в виде дополнительного столбца расширенной матрицы. Это означает, что можно одновременно решать несколько систем с одинаковой матрицей коэффициентов и различными правыми частями. Для этого достаточно составить расширенную матрицу вида

a	a	a	....	a		a	...	a
11	12	13		1n		1(n+1)		1(n+k )
B = a21	a22	a23	....	a2n		a2(n+1) ... a2(n+k )				(2.26)
...	...	... .... ...				...	...	...	,
	an2	an3	....	ann		an(n+1)	...
an1								an(n+k )

в которой для удобства дальнейших преобразований элементы m-й правой части обозначены как ai(n+m) , и выполнять все операции метода Гаусса над этой рас-

ширенной матрицей.

Метод оптимального исключения. В этой модификации метода Гаусса операции прямого и обратного хода метода выполняются попеременно. На первом шаге после приведения первого уравнения исключается неизвестное x1 из второго уравнения, а затем с помощью приведенного второго уравнения - неизвестное x2 из первого. После k таких шагов матрица системы имеет вид

	1
...
	0
a((kk+)1)1
...
	a( k )
	n1

... 0

... ...

... 1

... a((kk+)1) k

... ...

... ank( k))

a1((kk)+1)

...

ak( k( k) +1) a((kk+)1)(k +1)

...

an( k(k)+1)

... a1(nk )

... ...

... akn( k )

... a((kk+)1)n

... ...

... ann( k )

( k )
a1(n+1)
...		(2.26)
ak( k(n)+1)
a((kk+)1)(n+1)
...
a( k )
n(n+1)

33

На k+1-м шаге, используя первые k уравнений, исключают неизвестные x1 , x2 ,..., xk из (k+1)-го уравнения. Затем посредством этого уравнения исключа-

ется неизвестное xk+1 из первых k уравнений и т.д.

Таким образом, в результате прямого хода матрица системы приводится к диагональному виду с единицами на главной диагонали. При этом отпадает необходимость обратного хода, поскольку столбец правой части приведенной матрицы [ai((nn)+1) ]и является вектором решения.

Метод Гаусса-Жордана. В этой модификации метода Гаусса операции исключения переменных для каждого приводимого уравнения осуществляют не только ниже, но и выше главной диагонали. Операции с первым уравнением системы полностью аналогичны стандартной схеме. Второе уравнение системы после приведения и домножения на соответствующие коэффициенты вычитают не только из третьего и последующих уравнений, но и из первого. В результате k таких шагов получаем матрицу

	1	...	0	a(k)
				1(k+1)
... ... ...				...
	0	...	1	a(k )
				k (k+1)
	0	...	0	a(k)
				(k+1)(k+1)
... ... ...				...
	0	...	0	a(k )
				n(k+1)

...	a(k)
	1n
...	...
...	a(k)
	kn
...	a(k)
	(k+1)n
...	...
...	ann(k)

(k)
a1(n+1)
...
(k )
ak (n+1)	.	(2.27)
(k)
a(k+1)(n+1)
...
(k)
an(n+1)

На k+1-м шаге ведущим является элемент	a(k)	. Разделив k+1 урав-
	(k +1)(k+1)

нение системы на ведущий элемент, получают приведенное уравнение, посредством которого исключают неизвестное xk+1 из остальных уравнений системы.

Как и в методе оптимального исключения, матрица системы приводится к диагональному виду и вектором решения является столбец [ai((nn)+1) ].

LU – разложение или компактная схема Гаусса. Как отмечалось выше,

матрицу A можно представить в виде произведения нижней треугольной матрицы L (lower) и верхней треугольной матрицы U ( upper). Это разложение единственно с точностью до выбора диагональных элементов матриц. Рассмотрим разложение, которое соответствует полученному в схеме метода единственного деления метода Гаусса, где U - верхняя треугольная матрица с единичной диагональю:

a11	a12		....	a1n		l11		0	....	0		1	U12	....	U1n
a	a	22	....	a		l	l	22	....	0		0	1	....	U
11			....		2n	= 11			....				....	....		2n
.... ....				....		.... ....				.... ....					....
	an2		....				ln2		....			0	0	....	1
an1				ann		ln1				lnn

34

Это равенство равносильно следующей вестных элементов lij и Uij матриц L и U :

min(i, j)

aij = ∑lik ×U kj (i =1,..., n

k =1

системе относительно n2 неиз-

j =1,..., n) .

(2.28)

Специфика этой системы позволяет находить неизвестные одно за другим посредством алгоритма, который называется компактной схемой метода Гаусса или схемой Халецкого. Элементы lim и Umi (m=1,…,n) вычисляются по формулам (2.29) в следующем порядке: вначале элементы m-го столбца матрицы L, затем m -й строки матрицы U .

li1 = ai1,	U1 j	= a1 j / l11, i =1,...,n,		j = 2,...,n,
	m−1				(2.29)
lim = aim − ∑lik ×Ukm ,			m = 2,...,n,	i = m,...,n,
	k =1
	m−1
Umj = (amj − ∑lmk ×Ukj ) / lmm ,				j = m +1,...,n.
	k =1	Ax=b	сводится к последовательному решению
Тогда решение системы

двух систем:

Ly = b,

(2.30)

Ux = y.

В отличие от схемы единственного деления компактная схема менее удобна для усовершенствования с цель уменьшения влияния вычислительных погрешностей путем выбора ведущих элементов. Достоинством ее можно считать то, что LU – разложение играет роль обратной матрицы и может использоваться при решении серии систем с одной и той же матрицей и разными правыми частями.

Нахождение определителя и обратной матрицы. Как отмечалось выше,

особенностью решения СЛАУ методом последовательного исключения неизвестных является применение только элементарных преобразований системы– вычитания строк и деления строки на ведущий элемент в процессе приведении матрицы системы к треугольному виду. Это позволяет применить прямой ход метода Гаусса для нахождения определителя матрицы:

n

det A = (−1)p ∑akk(k−1) ,

k=1

где akk(k−1) - ведущий элемент на k – м шаге, а р – число перестановок строк мат-

рицы при выполнении прямого хода. Если ведущий элемент на каком-то шаге окажется равным нулю, то матрица - вырожденная, т.е. det A = 0.

Метод Гаусса можно использовать и для нахождения обратной матрицы для

невырожденной матрицы A. По определению обратная матица A−1 удовлетворяет условию:

35

	1	0	... ...		0
	0	1	0 ...		0
						.
AA−1 =E=...		...	... ...		...
	0	...	0	1	0
	0	0	...	0	1

Рассмотрим обратную матрицу как n столбцов:

A−1 =[Y ,Y ,...,Y ] ,	где Y - j-ый столбец A−1 .
1 2 n	j

Тогда

AA−1 =[AY1, AY2,..., AYn ]=E .

Таким образом, каждый столбец обратной матрицы является решением системы с одной и той же матрицей коэффициентов:

		1					0					0
		0		,	AY		1		,…,			0	.
AY =						=				AY =
1		...			2		...			n		...
		0					0					1

Как уже отмечалось выше, метод Гаусса позволяет одновременно решать несколько систем с одинаковой матрицей коэффициентов и различными правыми частями. Это связано с тем, что в преобразованиях прямого хода правая часть системы участвует только пассивно – в виде дополнительного столбца расширенной матрицы. Отсюда очевиден алгоритм построения обратной матрицы:

необходимо решить n систем с основной матрицей A и правыми частями, в которых один элемент равен 1, а все остальные равны 0. Для этого достаточно составить расширенную матрицу вида

a	...	a	1	0	...	0
11		1n
[A,E]= ... ... ...			...	...	...	...
	...	a	0	0	...	1
a
1n		nn

и выполнять все операции метода Гаусса над этой расширенной матрицей. На месте единичной матрицы будет получена обратная матрица A−1 .

Метод прогонки решения систем с трехдиагональными матрицами коэффициентов. Этот метод является модификацией метода Гауcса для систем специального вида с трехдиагональной матрицей A, т.е. с матрицей, все элементы которой, не лежащие на главной и двух побочных диагоналях, равны нулю (aij = 0 при j > i +1 и j < i −1).

36

Рассмотрим систему линейных алгебраических уравнений следующего

вида
b1 x1 + c1 x2									= d1 ,
a2 x1 + b2 x2 + c2 x3									= d2 ,
a3 x2 + b3 x3 + c3 x4									= d3 ,
LLLLLLLLLLLLLLLLLLLLLLLLLLLL(2.31)
			an−1 xn−2 + bn−1 xn−1 + cn−1 xn = dn−1 ,
						an xn−1 + bn xn			= dn .
Матрица А системы (2.31) имеет вид:
b			c	0		0 ...	0	0	0
	1		1	c		0 ...	0	0	0
a			b
		2	2		2	c3 ...
	0		a3	b3			0	0	0
A =									...	,
... ... ... ... ... ... ...
	0		0	0		0 ...	an−1	bn−1	cn−1
	0		0	0		0 ...	0	an	bn

но для записи такой системы не используют матричные обозначения для коэффициентов при неизвестных. На главной диагонали её расположены элементы b1,b2 ,...,bn , над ней - элементы c1,c2 ,...,cn−1, под ней - элементы

a2 ,a3,...,an .

Системы типа (2.31) получается при моделировании некоторых инже-

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

Метод прогонки состоит из двух этапов - прямой прогонки (аналог прямого хода метода Гаусса) и обратной прогонки (аналога обратного хода метода

Гаусса). Прямая прогонка состоит в том, что каждое неизвестное xi выражается

через xi+1 с помощью так называемых прогоночных коэффициентов Ai ,							Bi ра-
венствами
xi = Ai xi+1 + Bi ,				i =		.	(2.32)
					1, n −1
Из первого уравнения системы (2.31) имеем
x = − c1		x	2	+ d1 ,
1	b1			b1
а из соотношения (2.32) получаем, что
x1 = A1x2 + B1.							(2.33)

37

Отсюда и из предыдущего равенства находим, что
A = −			c1		,
1			b1
					.	(2.34)
B	=	d1
1		b1

Далее, из второго уравнения системы (2.31) и соотношения (2.33) получаем равенство

a2 (A1x2 + B1)+ b2 x2 + c2 x3 = d2 ,

Отсюда выразим x2 через x3 :

x2 =	− c2 x3 + d2 − a2B1							,
а из (2.32) имеем			a2 A1 + b2
x2 = A2 x3 + B2 .
Из этих двух равенств находим коэффициенты A2 и B2 :
A2		= −		c2			,
				a2 A1	+b2 .
B		=	d2 −a2 B1
2			a2 A1 + b2

Продолжая этот процесс далее, можно вычислить прогоночные коэффициенты Ai и Bi для любого номера i:

	A = −	ci	, B	=	di − ai Bi−1		,i =			(2.35)
								2, n −1.
	i	ai Ai−1 + bi	i		ai Ai−1	+ bi
	Таким образом, коэффициенты			Ai и Bi		уравнений (2.32),				связывающие

соседние значения xi и xi+1, i =1, n −1, можно определить из рекуррентных со-

отношений (2.35), поскольку начальные коэффициенты				A1 и B1, согласно
(2.34), известны.
Обратная прогонка состоит в последовательном вычислении неизвестных
xi . Сначала определяется xn .
Из (2.32) при i = n −1 получаем равенство
xn−1 = An−1xn + Bn−1 .
Последнее уравнение системы (2.31) имеет вид
an xn−1 + bn xn = dn .
Исключив из этих двух уравненийxn−1, найдём, что
x	=	dn − an Bn−1	.	(2.36)
n		bn + an An−1

Далее из (2.32) при i = n − 2 и последнего уравнение системы (2.31) получаем

38

xn−2 = An−2 xn−1 + Bn−2 ,

an−1 xn−2 + bn−1 xn−1 + cn−1 xn = dn−1 ,

где xn уже известно. Исключив из этих двух уравнений xn−2 , найдёмxn−1. Продолжая процесс далее, последовательно вычислим значение xn−1 ,..., x2 , x1 .

Заметим, что при нахождении прогоночных коэффициентов Ai и Bi по

формулам (2.35) надо учитывать возможность деления на нуль. Можно пока-

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

условия bi ≥ ai + ci , причём хотя бы для одного i имеет место строгое нера-

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

коэффициентов Ai не превосходят единицы.

Условие преобладания диагональных элементов обеспечивает также устойчивость метода прогонки относительно погрешностей округлений. Это позволяет использовать метод прогонки для решения больших систем уравнений. Отметим, что в ряде случаев для хорошо обусловленных систем вида (2.31) метод прогонки оказывается устойчивым и при нарушении условия преобладания диагональных элементов.

Число арифметических действий, выполняемых при решении системы (2.31) методом прогонки по формулам (2.32), (2.35) и (2.36), зависит линейно от числа неизвестных:

Q = 6(n −2) +5 + 2(n −1) =8n −9

При решении же произвольной системы линейных алгебраических уравнений методом Гаусcа число действий пропорционально кубу числа неизвестных.

Исследование погрешностей решения СЛАУ прямыми методами. Ошибки округления, совершаемыe в процессе вычислений, почти всегда приводят к то-

му, что вычисленное решение xr* в определенной степени отличается от теоре-
r	−1
тического решения x = A	b . Оценить погрешность вычисленного решения
xr* можно при помощи двух векторов: ошибки ∆xr = xr− xr* и невязки rr = Axr* −b .

Из теории матриц известно, что если одна из этих величин равна нулю, то и другая также должна быть равна нулю. Но они не обязаны одновременно быть “малыми”. Доказано, что метод Гаусса с частичным выбором ведущего элемента гарантированно дает малые невязки, но отсюда, к сожалению, нельзя сделать вывод, что и ошибка мала.

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

процесс решения исходной системы Axr = b с учетом погрешностей округления эквивалентен точному решению некоторой возмущенной системы

39

(A + ∆A)(xr+ ∆xr) = (br+ ∆br) . И каждому вычислительному алгоритму соответст-

вует своя матрица эквивалентных возмущений ∆A.

Рассмотрим системуr уравнений, правая часть которой получила приращение («возмущение»)r ∆b . Реакцией решения на возмущение правой части будет вектор поправок ∆x , т.е. решением «возмущенной» системы

A(xr + ∆xr) = (b + ∆b)

будет вектор «возмущенный» вектор x + ∆x .

Абсолютную погрешность решения будем оценивать нормой разности между точным и приближенным решениями ∆xr = xr − xr* , а относительную –

отношением абсолютной погрешности к норме решения (точного или приближенного).

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

r

∆b

∆rx

≤ (?)

r

x

b

с неизвестным пока коэффициентом.

Для этого подставим исходнуюr систему в «возмущенную» и получим уравнение для поправки ∆x :

A∆xr = ∆b ,

из которого найдем её явное выражение

r

−1

∆x = A

∆b . Нормируя это выражение

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

br

≤

A

xr

,

∆xr

≤

A−1

∆b

.

Эти неравенства одинаковогоr смысла можно перемножить

b ∆xr ≤ Axr A−1 ∆b

и разделить на b xr . Получим искомую оценку связи между величиной «возмущения» правой части системы и величиной реакции системы:

∆xr

A

A−1

∆b

(2.37)

r

≤

r

.

x

b

Коэффициент этой связи называется числом (мерой) обусловленности матрицы A и обозначают cond(A ) или M A

cond(A) = A A−1 .

Это же самое число является коэффициентом роста относительных погрешностей при неточном задании элементов матрицы исходной системы, т.е. для решения «возмущенной» системы

(A + ∆A)(xr + ∆xr) = b .

Справедливо неравенство

40

∆xr

∆A

.

(2.38)

r

≤ cond(A)

A + ∆A

x

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

∆x

≤

cond(A)

∆A

+

∆b

(2.39)

.

x

∆A

A

b

1

−cond(A)

A

Отсюда видно, что на точность решения влияют два фактора: число обусловленности матрицы и эквивалентные возмущения. Если cond( A) велико, то ма-

лым относительным погрешностям системы могут соответствовать большие относительные погрешности решения. Другими словами, неравенства (2.37) – (2.39) показывают, что чем больше число обусловленности матрицы A, тем сильнее сказывается на решении линейной системы ошибка в исходных данных. Системы с большим cond(A) называются плохо обусловленными. Плохо обу-

словленные системы рекомендуется решать, увеличивая, насколько можно, точность задания коэффициентов матрицы и правой части, а также проводить промежуточные вычисления с повышенной точностью.

Назвать «пороговое» значение, начиная с которого систему можно считать хорошо или плохо обусловленной, трудно. Здесь может играть роль размерность задачи, точность, с которое должно быть найдено ее решение, точность представления чисел в ЭВМ и т.д. Нижней границей числа обусловленности является 1, т.е. cond(A) ≥1 всегда. Можно также указать верхнюю границу

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

Установлено, что относительная невязка и относительные эквивалентные возмущения системы имеют величину, сравнимую с ошибкой округления, независимо от того, насколько плохо обусловлена матрица:

b − Ax

= O(β−t ),

δA

= O(n β−t ).

Ax

A

Здесь β − основание системы представления чисел с плавающей точкой, а t − число цифр мантиссы, используемых при представлении чисел в ЭВМ.

Таким образом, для относительной погрешности решения справедлива оследующая оценка

∆x

−t

(2.40)

x

= O(cond(A) n β

)

и из нее видно, что нельзя рассчитывать на получение результата с точностью

ε, если cond(A) ≥ ε n β −t .

41

2.3. Итерационные методы решения СЛАУ

Для решения итерационным методом исходная система линейных алгеб-
r		r
раических уравнений Ax	= b должна быть каким-либо способом (таких спо-
собов существует бесконечное множество) приведена к виду
	r	r	r	(2.41)
	x	= Gx	+ f ,
где G - некоторая матрица,	r	- преобразованный вектор свободных членов.
	f

Систему (2.41) можно рассматривать как задачу о неподвижной точке линейного отображения G и определить последовательность {xr(k ) }приближений к не-

подвижной точке в пространстве Rn рекуррентным равенством
xr(k ) = Gxr(k −1) + fr.	(2.42)

Для начала вычислений задается начальное приближение - произвольный век-
тор	xr(0)	- и строится	рекуррентная последовательность	векторов
xr(0) , xr(1) ,......, xr(k ) ,...... по формуле (2.42). Этот итерационный процесс называ-
ется методом простых итераций (МПИ). Если последовательность				{xr(k ) }схо-
дится к некоторому вектору xr			, то, согласно (2.41), этот вектор и есть решение
	r	r
системы Ax		= b .

Для выяснения условий применимости метода простой итерации надо найти ответы на следующие вопросы:

-Какие требования нужно предъявить к G, f , xr(0) , чтобы последовательность {xr(k ) }сходилась к решению исходной задачи?

-С какой скоростью сходится этот процесс?

-Сколько нужно сделать итераций, чтобы при заданном начальном приближении получить решение задачи с заданной точностью?r

	Ответы на эти вопросы содержатся в следующих теоремах.
	Теорема 2.1. Если || G ||<1, то система уравнений (2.41) имеет единст-
венное решение xr * и последовательные приближения				xr(k ) сходится к реше-
нию x* со скоростью геометрической прогрессии.
ство	∆ Допустим, что xr* - решение системы (2.41), то есть имеет место равен-
	r	r	r
				(2.43)
	x*	= Gx	* + f .
Отсюда, используя неравенство треугольника для нормы вектора, получаем
	|| xr* ||=|| Gxr* + fr ||≤|| Gxr* || + ||		fr ||≤|| G || || xr* || + || fr ||, т. е.

|| xr* ||≤|| G || || xr* || + || fr || . Отсюда

42

r*		|| f ||
|| x	||≤		.		(2.44)
		1−|| G ||
Из этого неравенства следует, что система x = Gx					имеет только нулевое
			r	r	r
решение. Следовательно, неоднородная система x				= Gx	+ f имеет единствен-

ное решение при любом свободном f .

Пусть теперь rrk = xr* −xrk , где xr* - решение системы (2.41). Вычитая из равенства (2.43) равенство (2.41), будем иметь

rrk = xr* −xrk = Gxr* −Gxrk −1 = G(xr* −xrk −1) = Grrk −1, т.е.

rrk = Grrk −1 = G2rrk −2 =... = Gk rr0 ,

где rr0 = xr* − xr0 - погрешность начального приближения. Отсюда и из (2.41) получаем неравенство || rrk ||=|| Gk rr0 ||≤|| G ||k || rr0 || .

Так как || G ||<1, то из полученного неравенства следует, что rrk стремит-

ся к нулю при k → ∞ не медленнее геометрической прогрессии со знаменате-

лем q =|| G ||<1.□

Из теоремы 1 вытекает

Следствие: Если для матрицы G = [gij ]выполнено одно из условий:

	n
1)	| G ||1 ={ ∑| gij |} <1;
	1≤ j ≤n
	n n	12
2)	|| G ||2 ={∑∑gij2}		<1;
	i =1 j =1
	n
3)	|| G ||3 ={∑| gij |} <1,
	i =1

то последовательность ( x ), построенная по рекуррентной формуле (2.42), сходится к единственному решению системы (2.41).

Теперь рассмотрим проблему оценки погрешности итераций.

Пусть || G ||<1, тогда, согласно теореме 1, система (2.41) имеет единст-

венное решение xr ,	для которого справедливо равенство (2.42). Из равенств
(2.41) и (2.43) имеем:	= G(xr* − xr(k −1) ) xr* = xr(k ) +G(xr* − xr(k −1) ) .
xr* − xr(k )		(2.45)

Вычтя xr(k −1) из обеих частей последнего равенства, будем иметь

x* − xr(k −1) = xr(k ) − xr(k −1) +G(xr* − xr(k −1) ) .

Отсюда, используя свойства нормы, получаем

|| x* − xr(k −1) ||≤|| xr(k ) − xr(k −1) || + || G || || xr* − xr(k −1) ||

43

или

||≤|| xr(k ) − xr(k −1) || .

(1−|| G ||) || x* − xr(k −1)

Так как || G ||<1, то из последнего неравенства следует, что

*

r(k −1)

1

r(k )

r(k −1)

|| x

− x

||≤

||x

−x

|| .

(2.46)

1−|| G ||

Из (2.45) находим

|| xr* − xr(k ) ||≤|| G || || xr* − xr(k −1) ||.

Отсюда и из неравенства (2.46) будем иметь

|| x* − xr(k )

||≤

G

||xr(k ) −xr(k −1)

||.

(2.47)

1−|| G ||

Неравенство (2.47) позволяет оценить норму погрешности k -го приближения через норму разности двух последовательных приближений и норму матрицы G. Эта оценка называется апостериорной и широко используется на практике. Если применить это неравенство k раз, то получим априорную оценку погрешности

|| x* − xr(k ) ||≤

G

k

||xr(1) −xr(0) ||

(2.48)

1−|| G ||

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

Теорема 2.2. Для сходимости метода простой итерации (МПИ) при любом

начальном приближении xr(0) к решению системы (2.41) необходимо и достаточно, чтобы все собственные значения матрицы G были по абсолютной величине меньше единицы.

На практике это трудно проверить, и обычно пользуются достаточными условиями сходимости – итерации сходятся, если какая-нибудь норма матрицы меньше единицы, т.е.

G

n

G

n

1

= max ∑

gij

= α < 1

или

2

= max ∑

gij

= β <1.

1≤i≤n j=1

1≤ j≤n i =1

Чем меньше норма матрицы G, тем быстрее сходится итерационный процесс. Рассмотрим один из простейших способов преобразования системы к ви-

ду (2.41) - решим каждое i-е уравнение исходной системы относительно xi :

1			n
xi = −			∑aij x j
	a

ii j =1, j ≠i

−bi . (2.49)

Рассмотрим два классических метода, использующих преобразованную таким образом систему - метод Якоби и метод Зейделя.

Метод Якоби использует следующий алгоритм построения приближений:

44

		1		n		.
	xi(k ) = −			∑aij x(jk −1)	−bi		(2.50)
		ii
		a	j =1, j ≠i

n

Если A - матрица с доминирующей диагональю, т.е. aii > ∑ aij , то метод

j ≠i

Якоби сходится при любом начальном приближении x( 0 ) .

Модификация метода Якоби, в которой при k-й итерации для вычисления следующего xi(k) используются уже вычисленные на этом шаге x1(k),x2(k),.,xi-1(k),

называется методом Зейделя:

xi(k ) = − 1 ∑i−1 aij x(jk ) aii j =1

n
+ ∑aij x(jk −1)	−bi	,	(2.51)
j =i +1

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

вое значение xi(k) запоминается на месте старого.

Приведенные выше методы Якоби и Зейделя относятся к одношаговым итерационным методам, когда для нахождения x(k+1) требуется помнить только одну предыдущую итерацию x(k). Для исследования сходимости удобнее записывать итерационные методы не в координатной, а в матричной форме, придерживаясь стандартной формы записи итерационных методов.

Канонической формой одношагового итерационного метода ре-

шения СЛАУ называется его запись в виде

	Bk +1	x(k +1) − x(k )	+ Ax(k ) = b ,	(2.52)
		τk +1

где Bk+1 - матрица, задающая тот или иной итерационный метод, τk+1 - итерационный параметр. Числовые параметры τk вводят для ускорения сходимости. Способ выбора итерационных параметров определяется при исследовании сходимости метода, когда выясняется при каких значениях параметров метод сходится и когда сходимость будет наиболее быстрой (соответствующие параметры называются оптимальными).

Итерационный метод называют явным, если Bk+1 - единичная матрица. Неявные итерационные методы имеет смысл применять лишь в том случае, когда решение системы уравнений с матрицей Bk требует меньше машинной памяти или времени или алгоритмически проще, чем решение исходной системы.

Итерационный метод называется стационарным, если Bk+1=B и τk+1=τ, (т.е. не зависят от номера итерации), и нестационарным в противоположном случае. Согласно этой классификации метод Якоби является одношаговым явным стационарным методом с диагональной матрицей B=DA (bii.=aii) и τ =1, метод Зейделя - неявным методом с нижней треугольной матрицей B= DA+LA.

45

Методом простой итерации называют явный метод с постоянным параметром

x(k +1) − x(k )	+ Ax(k ) = b , или	x(k +1) = x(k ) −τ r(k ) ,
τ

где r(k) = Ax(k)-b - вектор невязки. Метод сходится для симметричных поло-

жительно определенных матриц при 0 < τ < 2A .

Для окончания итерационного процесса на практике используют три спо-

соба. При первом определяют величину стабилизации и прекращают вычисления, если она меньше ε, т.е.

x(k +1) − x(k )

max{x(k ) , x(k +1) }≤ ε .

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

При втором способе вычисляют нормы невязки до начала итераций и на каждой итерации. Итерации прекращают при выполнении неравенства

Ax(k ) −b

Ax(0) −b ≤ ε .

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

x(k ) − x ≤ qk x(0) − x ,

где q (0,1), то говорят, что метод сходится со скоростью геометрической про-

грессии со знаменателем q. Можно определить, потребовав, чтобы qn < ε, число итераций n, достаточное для того, чтобы начальная погрешность уменьшилась в заданное число раз:

n ≥ n0	(ε) =	ln(1	/ ε)	.
		ln(1
			/ q)

Целая часть числа n0(ε) называется минимальным числом итераций, необходимым для получения заданной точности ε.

Величину ln(1/q) называют скоростью сходимости итерационного

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

46

точности. Качество различных итерационных методов сравнивают обычно по их скорости сходимости: чем выше скорость сходимости, тем лучше метод.

2.3 Вычисление собственных значений и собственных векторов

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

a − λ a

... a

11

12

1n

ϕ(λ) =

A − λE

=

a21

a22

− λ

...

a2 n

=

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

(2.53)

a

a

...

a

− λ

n1

n 2

nn

= (−1)n (λn − p λn−1

− ... − p

n−1

λ

− p

n

= 0

1

Определение компонент собственного вектора требует решения системы п однородных уравнений с п неизвестными; для вычисления всех собственных векторов матрицы, вообще говоря, требуется решить п систем вида:

(A −λi E)Xi = 0 ,

(2.54)

где Xi = (x1i , x2i ,...,xni ) — собственный вектор матрицы А, принадлежащий соб-

ственному значению λ i .

Коэффициенты pi характеристического полинома являются, с точностью

до знака, суммами всех миноров определителя матрицы А порядка i, опирающихся на главную диагональ. Непосредственное вычисление коэффициентов pi

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

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

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

47

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

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

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

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

Метод А. М. Данилевского.

Элегантный и весьма эффективный метод вычисления коэффициентов характеристического полинома предложен А.М. Данилевским. Геометрический

смысл этого метода состоит в следующем. Данная матрица A рассматривается как матрица оператора в базисе:

e1 = (1,0...0)′, e2 = (0,1...0)′,..., en = (0,0...1)′

Предполагается, что векторы e1 , Ae1 , A n −1 e1 линейно-независимы. Тогда

An e1 = p1 A n −1 e1 + p2 A n −2 e1 + pn e1

Ясно, что коэффициенты p1, p2 ,...pn - искомые коэффициенты характеристического полинома.

В базисе e1, Ae1 ,A n −1 e1 рассматриваемый оператор, очевидно, будет иметь так называемую матрицу Фробениуса

48

		0	0 ...	p n
			0 ...	p n −1
	1
P =		0	1 ...	p n − 2		,
	.......... ..........
		0	0 ...	p1

содержащую в явном виде искомые коэффициенты характеристического полинома.

Переход от базиса e1 , e2 ,..., en к базису e1, Ae1 ,A n −1 e1 осуществляется по-

степенно	в п-1 шагов.	Каждый шаг состоит в	переходе от базиса
e1 , Ae1 ,	A k −1e1 , ek +1 ,..., en	к базису e1 , Ae1 , A k −1e1 ,	A k e1 , ek +2 ,..., en .

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

Обозначим через A(k ) матрицу, полученную при к-1-м шаге процесса, так что A = A(1) , P = A(n) . Столбцами матрицы A(k ) являются координаты век-

торов e1 , A2 e1 , ... , A k e1 , Aek +1 ,..., Aen	в базисе e1 , Ae1 , A k −1e1 , ek +1 ,..., en .
Поэтому первые к -1 столбцов матрицы	A(k ) будут совпадать с одноименны-

ми столбцами матрицы Фробениуса P . Имеем A(k +1) =Sk−1 A(k ) Sk гдеS k - соответствующая матрица преобразования координат. Ясно что

			1 ... s1			, k +1 ... 0
S		=		0 ... s	2	, k +1 ... 0
	k		.......... ........
				0 ... s			... 1
					n	, k +1

где s1,k +1 , s2,k +1 ,...sn ,k +1 координаты вектора A(k )e1 в базисе e1 , Ae1 , A k −1e1 , ek +1 ,..., en . Эти координаты, как мы видели выше, есть не что

иное, как элементы aik(k ) k-ого столбца матрицы A(k ) Имеем далее:

49

B(k )

A(k )

		1 ...	−		s 1 , k + 1		...			0
				s k + 1 , k + 1
					s 2 , k	+ 1
		0 ...	−						...	0
					s k + 1 , k + 1
S k− 1 =	..........			..........		.....
				1
		0 ...						...		0
					s k + 1 , k + 1
	..........			..........		.....
					s n , k
		0 ...	−			+ 1			...	1
					s k + 1 , k + 1
Вычисление матрицы	A(k +1)		целесообразно проводить в два приема. В на-

чале вычисляется вспомогательная матрица B(k ) = Sk−1 A(k ) . Это действие, в си-

лу строения	матрицы S k−1 , состоит в том, что (k+1)-я строка матрицы A(k ) ум-
ножается на			1	, а от каждой из остальных строк вычитается полученная
		a	(k )
			k +1,k

(k+1)-я строка матрицы B(k ) , умноженная на соответствующий элемент k-го

столбца матрицы . В результате этих действий первые (k-1) столбцов не изменятся, в k-ом же столбце на (k+1)-м месте окажется единица, а все осталь-

ные элементы станут нулями. Вычисления остальных элементов матрицы будут двухчленными, напоминающими вычисления метода Гаусса. Полученная

матрица B(k ) умножается затем справа наS k . При этом изменяется только один, именно (k+1) столбец. Его элементами будут:

n

∑b1(kj )a(jkk ) j =1

n

, ∑b2(kj)a(jkk ) j =1

n

, ... , ∑bnj(k )a(jkk ) j =1

Иначе говоря, (k+1)-й столбец матрицы A(k +1) есть результат итерации k-го столбца матрицы A(k ) матрицей B(k ) .

Таким образом, переход от матрицы A(k ) к матрице A(k +1) происходит по формулам:

50

b(k )		=		1	a(k )		,
k +1, j					k +1, j
			ak(k+1,) k
b(k )	= a		(k )	− a(k )b(k )				,	(i ≠ k +1)
i, j		i, j			i,k	k +1, j
		a(k +1)		= b(k ) ,					( j ≠ k +1)
		i, j			i, j

n

ai(,kk++1)1 = ∑bi(,kj) a(jk,k) . j=1

Число операций, нужных для вычисления по методу А.М. Данилевского, значительно меньше, чем по другим методам – общее число умножений и деле-

ний равно n3 − n2 .

Таким образом, после нахождения матрицы A(k +1) , характеристический многочлен выписывается по виду последнего столбца.

Т.к A(k +1) y = λy и A(k +1) =Sk−1 A(k ) Sk ,

то Sk−1 A(k ) Sk = λy .

Умножая это равенство слева на матрицу Sk , получим:

A(k ) Sk y = λSk y .

Таким образом, собственные векторы исходной матрицы A легко находятся по соответствующим собственным векторам ее канонической формы Фробениуса. Если λ - известное собственное значение матрицы Фробениуса, то

A(k +1) y = λy

p1 y1 + p2 x2 +... + pn yn = λ y1 y1 = λ y2

y2 = λ y3

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

yn−1 = λ yn

Положим yn =1 и найдем все остальные значения: xi = Sk yi = M n−1M n−2 ...M1 yi

где матрицы M i (i =1,...n −1) отличаются только одной строкой от единичной матрицы.

Метод А. Н. Крылова.

Идея А. Н. Крылова заключалась в предварительном преобразовании уравнения (1) в эквивалентное ему уравнение вида

51

b11 − λ		b12	... b1n
b	− λ2	b	... b
21		22	2 n
D(λ) = ................................. = 0 ,					(2.55)
	− λn	b	... b
b
n1		n 2	nn
развертывание которого по степеням λ осуществляется, очевидно, значительно
проще, при помощи разложения определителя по минорам 1-ого столбца.
Равенство нулю определителя (2.63) есть необходимое и					достаточное ус-
ловие для того, чтобы система однородных уравнений

λx1 = a11 x1 + a12 x2 +... + a1n xn ,
λx		= a x + a x +...		+ a x ,
	2	21 1	22 2	2n n
.................................................					(2.56)
λx		= a x + a x +...		+ a x
	n	n1 1	n2 2	nn n
имела решение х1, х2, ...,	хп, отличное от нулевого.

Преобразуем систему (2.56) следующим образом. Умножим первое уравнение на λ и заменим λ х1, ..., λ хп их выражениями (2.56) через х1, ... , хп. Это дает

где		λ2 x1 = b21 x1 +b22 x2 +... + b2n xn ,						(2.57)
				n
		b2k = ∑a1s ask .						(2.58)
				s=1
Умножим далее уравнение (2.57) на λ и снова заменим λ х1, ...,								λ хп их выраже-
ниями через х1, ... , хп :	λ3 x1 = b31 x1 +b32 x2 +...+b3n xn .
Повторяя этот процесс (n -1) раз, мы перейдем от системы (2.56) к системе
λx1		= b11 x1 + b12 x2 +... + b1n xn ,
λx		= b	x +b		x +... + b		x ,
	2	21	1	22	2	2n	n
.................................................								(2.59)
λx		= b	x +b		x +	... + b	x
	n	n1	1	n2	2	nn	n

коэффициенты которой bik будут определяться по рекуррентным формулам

52

b1k	= a1k ,
	n	(2.60)
bik	= ∑bi−1,s ask , i = 2,..., n; k = 1,..., n

s=1

Очевидно, что определитель системы (2.59) будет иметь вид (2).

Система (2.59) имеет ненулевое решение для всех значений λ , удовлетворяющих уравнению (1). Таким образом, D(λ) обращается в нуль при всех λ , являющихся корнями уравнения (1). Покажем, что

				1 0	...	0
D(λ)				= b11	b12	... b1n		= N
ϕ	(	λ	)	.................................
					bn−1,2 ...
				bn−1,1			bn−1,n

т. е. при N ≠ 0		D(λ)	отличается от искомого характеристического полинома
только численным множителем.
Введем в рассмотрение векторы Bi c компонентами bi1						bi 2 ,..., bin .
	n
Равенства bik = ∑b i−1,s a sk				показывают, что Bi = A′Bi −1		(i=2,..,n),
s =1					B = (1,0,...,0)′.
		B = A′B , где
	i		1	0	0
′		B0	(i =1,..., n) .
Тогда Bi = ( A )

Итерационный степенной метод.

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

Пусть дана матрица A и пусть ее собственные значения, соответствующие собственным векторам x(1) , x(2) , x(3) ...x(n) , упорядочены по абсолютным вели-

чинам

λ

>

λ

2

≥

λ

3

≥ ... ≥

λ

n

. Тогда, выбрав некоторый вектор y(0)

, напри-

1

мер, вектор, компоненты которого равны единице y(0) = (1,1,...1)T , можно для определения λ1 построить следующий итерационный процесс:

53

y(1) = Ay(0)	λ(1)	=		y(1)j
	1		y(j0)
y(2) = Ay(1)	λ(2)	=		y(2)j
	1			y(j1)				,
....................................
y(k) = Ay(k−1)	λ(k) =				y(jk)
	1				y	(jk		−1)

где y(jk −1) , y(jk ) - соответствующие компоненты векторов yr(k −1) , yr(k ) . При этом в качестве номера j может использоваться любое число из диапазона j = 1,…,n. В связи с тем, что вектор y(k ) на k -ой итерации может быть представлен в виде y(k ) = Ay(k −1) = Ak y(0) , рассматриваемый итерационный процесс носит название

"степенной метод".

Покажем, что последовательность y(k ) сходится к собственному вектору x(1) . Поскольку {x(1) , x(2) ,..., x(n) } - базис пространства Rn, то имеет место представление:

y(0) = α1 x(1)

+α2 x(2) +... +αn x(n)

Предположим, что α ≠ 0. Согласно определению

Ax(i)

= λ

x(i) , i = 1,…,n.

1

i

Тогда A(k) x(i) = λ k x(i)

и имеет место разложение:

i

y(k ) = Ak y(0) = α1λ1k x(1) +α2 λ2k x(2)

+... +αn λnk x(n) .

(2.61)

Представим его в виде

λ2 )k

λn

λ −k y(k )

= α

x(1) +

α

(

x(2) +... +α

n

(

)k x(n)

1

1

2

λ

λ

и перейдем к пределу при k → ∞ .

1

1

Поскольку

λ2 k → 0,...,

λn

k

→ 0,

k → ∞

, то

λ

−k y(k )

→ α

1

x(1)

. Это зна-

λ1

1

λ1

чит что последовательность { y(k ) } сходится к собственному вектору x(1) по направлению, т.е. вектор y(k ) можно считать приближением к вектору α1λ1k x(1) ,

пропорциональному собственному вектору x(1) .

Построим числовую последовательность, сходящуюся к собственному значению λ1 на основе информации об y(k) . Учитывая разложение (2.61), обра-

зуем скалярные произведения (считаем, что система собственных векторов ортонормирована):

54

( y ( k ) , y ( k ) ) = α1

2 λ1

2 k + α2

2 λ2

2k

+ ... + αn

2 λn

2 k

=

λ

2 k (α

1

2 + δ

1

(k ))

,

δ

1

(k ) = α

2

2 (

λ2

) 2k + ... +

α

2

(

λn

)2 k ;

1

λ

1

n

λ

1

( y ( k ) , y

( k +1) ) = α

1

2 λ

2 k +1

+ α

2

2 λ

2 k +1

+ ... + α

2

λ

2 k +1 =

1

2

n

n

λ

2 k (α

1

2 λ +

δ

2

(k ))

,δ

2

(k ) = α

2 λ

2

( λ2

) 2 k + ... + α

n

2 λ

n

(

λn ) 2k .

1

1

2

λ

λ

1

1

Отметим, что при k → ∞,

δ1 (k) → 0 ,

δ2 (k) → 0 .

Положим

µ k

=

( y ( k ) , y ( k +1) )

( y

( k )

, y

( k )

)

Согласно предыдущим представления:

µ

=

α 2 λ

+δ

(k)

k

1

2

1

2

α1

+δ1 (k)

Следовательно, имеет место сходимости µk → λ1, k → ∞ , т.е. величина µk может служить приближением к искомому собственному значению λ1 . Отметим, что ρ(A) = λ1 . Следовательно, величина µk для больших k является приближенным значением для спектрального радиуса матрицы А.

Пусть r(x) = (x, Ax) , x ≠ 0 - отношение Релея. Поскольку y(k +1) = Ay(k ) , то

(x, x)

µk = ( y(k) , Ay(k) ) = r( y(k) ) .

( y(k) , y(k) )

В практических вычислениях степенной метод сопровождается процедурой евклидовой нормировки последовательных приближений y(k ) . Тогда модифицированная схема степенного метода имеет вид:

y(0) Rn , y(k) =		Ay(k−1)		, µk = ( y(k) , Ay(k) ), k =1,2,K
		Ay(k−1)

Пусть A - симметричная, невырожденная матрица. Поставим задачу отыскания минимального по модулю собственного значения λn , в предположении его

единственности. Рассмотрим обратную матрицу A−1 . Как известно, обратная

величина 1 есть максимальное по модулю собственное значение матрицы

λn

A−1 с соответствующим собственным вектором x(n) . Поэтому для отыскания

55

			1
	собственной пары			, x	(n)	можно использовать степенной метод примени-
			λn

тельно к матрице A−1 :

y(0) Rn , y(k) = A−1 y(k−1) , k =1,2,K

y(k) ≈ αxn ,	µk ≈	1	.
		λn

Отметим, что вектор y(k ) является решением линейной системы Ay = y(k −1) . Таким образом, степенной метод для отыскания собственной пары ( λn , x(n) ) имеет вид (обратные итерации):

y(0) Rn , Ay(k) = y(k −1) , k =1,2,K.

Рассмотрим также случай, когда матрица A имеет два наибольших по модулю собственных значения и эти значения вещественны и противоположны по знаку. Будем считать, что собственные значения матрицы A распределены по модулю следующим образом λ1 = λ2 > λ3 ≥ ... ≥ λn и λ1 = −λ2 .

Тогда, в силу формулы (2.61), получим:

( y(k) , y(k) ) =α				2λ	2k +α		2	λ			2k + ... +α				2	λ	n	2k =[λ	= −λ			2	]
			1	1			2			2					n			1
λ	2k (α 2	+α	2	+ δ	1	(k))	,δ		1		(k) =α	2	(	λ3 )2k + ... +α					2	(	λn )2k ;
1	1		2									3		λ					n		λ
															1							1

( y(k) , y(k +1) ) =α					2	λ 2k +1 +α				2	2λ		2	2k +1 + ... +α				2λ	2k
				1			1											n	n
λ	2k (α	2λ	−α	2	λ		+δ	2	(k)) ,δ			2		(k) =α	3	2λ	(	λ3 )2k
1	1	1		2	1											3		λ
																		1

+1 =[λ1 = −λ2 ]

+ ... +αn2λn (λn )2k .

λ1

Отсюда видно, что скалярные произведения ( y(k) , y(k) ) и ( y(k) , y(k +1) ) нельзя одновременно использовать для определения λ1 , ибо у них различные главные

части. Однако, в этом случае можно определить λ																	2 .
Рассмотрим:																1
			2		2k +α			2λ		2k +... +α				2		2k
( y(k) , y(k) ) =α				λ			2								λ		= [λ	= −λ			2	]
			1	1						2				n		n	1
λ 2k (α	2 +α	2	+δ		1	(k))	,	δ	1	(k) =α	2	(	λ3 )2k +... +α					2	(	λn )2k ;
1	1	2									3		λ					n		λ
														1						1

( y(k) , y(k+2) ) =α					2	λ	2k+2		+α		2	λ	2k+2	+... +α					2	λ	2k+2	= [λ
					1	1				2			2						n		n	1
λ 2k+2	(α	2	+α	2	+δ		2	(k))	,δ		2	(k) =α		2	λ	3	(	λ3 )2k			+2 +... +α		n
1		1		2										3				λ
																		1

= −λ2 ]

2λn (λn )2k .

λ1

56

Отметим, что при k → ∞, δ1 (k) → 0 , δ2 (k) → 0 . Следовательно:

( y(k) , y(k+2) )		λ 2k+2 (α			2	+α		2	+δ		1	(k))	= λ 2
	=	1			1			2
( y(k) , y(k) )			2k (α		2 +α			2 +δ			(k))
		λ		1			2			1			1
		1

Для нахождения собственных векторов, принадлежащих λ1 и λ2 целесообразно построить:

y(k ) + λ1

y(k −1)

= α

1

λ k x(1)

−α

λ k x(2)

+... +α

n

λ

k x(n) +

1

2 1

n

k ),

+α λ k x(1)

+α

2

λ k x(2)

+... +α

λ λ

k−1x(n) = λ k (2α

1

x(1)

+ O(

λ3

1

1

1

n 1

n

1

λ1

y(k ) - λ1

y(k −1)

= α

λ k x(1)

−α

λ k x(2)

+... +α

n

λ

k x(n) +

1 1

2 1

n

k ).

+α

λ k x(1)

−α

2

λ k x(2)

+... +

α

λ λ

n

k−1x(n) = λ

k (2α

2

x(2)

+ O(

λ3

1 1

1

n 1

2

λ1

Отсюда видно, что с точностью до постоянного множителя в качестве собствен-

ного вектора,		отвечающему λ1 , можно приближенно взять	вектор
y(k) + λ y(k −1)	,	а в качестве собственного вектора, отвечающему λ2 ,	можно
1

приближенно взять вектор y(k) - λ1 y(k−1) .

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

λ1 = λ2 > λ3 ≥ ... ≥ λn ,

λ1 = reiθ ,λ2 = re−iθ .

В силу формулы (2.61) имеем:

y(k ) = α λ k x(1)	+α	λ k x(2)	+... +α	λ k x(n)	(2.62)
1 1		2 2		n n

Пусть α1x(1) и α2 x(2) отличны от нуля. Поскольку матрица A и начальный вектор y(0) вещественны, то вещественным будет и вектор y(k ) = Ak y(0) . Значит, в

формуле (2.62) величины α1x(1) и α2 x(2) должны быть комплексносопряженными. Пусть α1x(1) = Reiµ и α2 x(2) = Re−iµ . Теперь из формулы

(2.62) получим:

y(k ) = 2Rrk cos(kθ + µ) +... +αn λnk x(n)

(2.63)

57

При достаточно большом k в силу условия

λ2

>

λ3

из формулы (2.63) полу-

чим:

k )

y(k ) = 2Rrk cos(kθ + µ) + O(

λ3

(2.64)

Наряду с равенством (2.64) будут иметь место следующие аналогичные равенства:

y(k +1)

= 2Rrk +1 cos((k + 2)θ + µ) + O(

λ3

k +1 ) ,

(2.65)

y(k +2)

= 2Rrk +2 cos((k + 2)θ + µ) + O(

λ3

k +2 ) .

(2.66)

Введем в рассмотрение определитель:

I (k ) =

y(k )

y(k+1)

y(k+1)

y(k+2)

и вычислим его значения, используя формулы (2.64)-(2.66). Получим:

I (k) = 4R2 r 2k+2{cos((k + 2)θ + µ) cos(kθ + µ) − cos2 ((k +1)θ + µ)} + Mr k O(

λ3

k ) =

= −4R2 r 2k+2 sin 2 θ + Mr k O(

λ3

k ),

где M – некоторая константа,

sin2 θ ≠ 0 , ибо λ1 и λ2

по предположению ком-

плексные числа. Аналогичным путем получим:

I (k −1) = −4R2 r2k sin2 θ + Mrk −1O( λ3 k )

Значит, модуль комплексного числа λ1 можно приближенно вычислить по формуле:

r	2	≈	I (k )	≈	y(k ) y(k +2) − y(k +1) y(k +1)		.
			I (k −1)		y(k −1)	y(k +1) − y(k ) y(k )

После того, как мы найдем r , аргумент комплексного числа можно находить приближенно по формуле:

y(k +2) + r 2 y(k )

cosθ ≈ 2ry(k +1) ,

т.к.

58

y(k +2) + r 2 y(k) = 2Rr k +2{cos((k + 2)θ + µ) + cos(kθ + µ)} + MO(

λ

k

3

= 2Rr k +2 cos((k +1)θ + µ) cosθ + MO(

λ

3

k

= 2ry(k +1) cosθ + MO(

λ

3

k .

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

λ1 = r(cosθ + i sinθ) , λ2 = r(cosθ −i sinθ).

59