Вычислительная математика лекции

			1		0		0 a11						a12					a13
1	1			0	1		0		0 a22				l21a12									U
L2	L1		A														a23 l21a13
				0	l		1		0 a				l a				a l a
					32						32			31	12		33		31	13
		a11			a12									a13
U			0 a		l a						a	23		l a						.
				22	21	12								21		13
			0		0			a l a					(a		23	l a )l			23
									33	31	13						21	13

Следовательно, A=LU, где L=L1L2. Можно показать, что

	1	0	0
L l		1	0
	21
		l32	1
l31

Существенным условием возможности такой записи L является нижняя треугольная форма с единичной главной диагональю для матриц

L11 и L21. Кроме того каждая из этих матриц должна иметь единственный столбец с ненулевыми элементами ниже главной диагонали.

Решение системы Ax=b осуществляют в следующей последовательности:

1.Производят факторизацию A=LU. Эта операция требует приблизительно 2n3/3 операций.

2.Решают систему с треугольной матрицей Ly=b. Необходимое количество операций приблизительно n2/2.

3.Решают систему с треугольной матрицей Ux=y. Необходимое количество операций приблизительно n2/2.

Общее число операций как и прежде имеет порядок 2n3/3. 2.Метод Гаусса с выбором главного элемента по столбцу (схема

частичного выбора).

В процессе расчета коэффициентов lij=aij/ajj, i=j+1, …, n, может оказаться, что ведущий элемент ajj = 0. В этом случае алгоритм становится некорректным, а при условии ajj aij вычислительно неустойчивым.

	10 5	1 x			1		с точным решением
Рассмотрим систему		1
	2	1 x2			0

x1=-0.499995…, x2=0.999995…. . Пусть длина слова десятичного

131

компьютeра соответствует четырем знакам. Иными словами числа представлены в форме 0.****×10p. Формально мы можем воспользоваться схемой единственного деления, так как a11≠0. Множитель

l21= 0.2 *101 0.2 *106 вычисляется точно. Значение нового

0.1*10 4

a22 =0.1*101-(-0.2*106)(0.1*101)=0.1*101+0.2*106=0.2000*106. На самом деле a22 =0.200001*106, но машинное слово вмещает только четыре цифры. При вычислениях нового b2 и компонентов x1, x2 в процедуре обратного хода ошибки округления не возникают:

b2 =-(-0.2*106)(0.1*101)=0.2*106

x			0.2	*106		0.1*101
				*106
2			0.2
x		0.1*101			0.1*101		0
1				0.1*10 4

Ошибка округления в шестом знаке после запятой при расчете а22

привела к полному искажению результата. Отметим, что cond(A)=2, задача хорошо обусловлена.

Для исключения подобных ситуаций осуществляют следующую модификацию алгоритма. Во время обработки очередного j-ого столбца просматриваются все его элементы, находящиеся на главной диагонали и

ниже ее. Выбирается элемент aij из условия | aij | max | aij | , после чего

j i n

уравнения с номерами i и j меняются местами. Дальнейшие вычисления идут по прежнему сценарию.

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

	2		1 x			0
	10	5	1			1
			1 x2

Решение дает следующие результаты:

l21 0.1*10 4 0.5*10 5 0.2 *101

132

Новое а21=0.1*101-(-0.5*10-5)(1)=0.1*101 (ошибка округления)

Новое b2=0.1*101-(-0.5*10-5)(0)=0.1*101

	x		0.1*101	0.1*101
		0.1*101
	2
	x	(0.1*101)(1)			0.5
	1	0.2 *101

Перестановка i и j строк эквивалентна умножению преобразуемой матрицы слева на матрицу перестановки Pij. Pij получается из единичной матрицы перестановкой ее i и j-ой строк . Для матриц перестановки справедливо соотношение Pij= Pij-1, Pij Pij=E.

Процесс формирования LU разложения модифицируется следующим образом. При обработке очередного k-ого столбца сначала необходимо осуществить перестановку выбранных строк, т. е. сформировать матрицу PkAk. Для этой матрицы найти Lk-1, а затем Ak= Lk-1PkAk-1. Например,

цепочка преобразований при n=4 может выглядеть следующим образом

A→A1→A2→A3=U.

L3-1P3L2-1P2L1-1P1A=U; A1=L1-1P1A; A2=L2-1P2A1; A3=L3-1P3A2.

Так как матрицы P3L2-1 и P2L1-1 не являются нижними треугольными с единичной диагональю, получить разложение A=LU из последнего соотношения не удается. Однако это соотношение можно записать в

следующей форме

L3-1 (P3L2-1 P 3 ) (P3 P2L1-1 P2P3) ( P3P2 P1A)=U . Обозначим

1						1
L1	(P P L 1P P ),					L2	(P L		1P ) , ( P3P2 P1A)=PA.
	3	2	1	2	3		3	2	3

1

Рассмотрим матрицу L1 (P3P2L1 1P2P3 ) . Здесь L1-1 нижняя треугольная матрица с единичной главной диагональю и с единственным, а

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

133

			1	0	0	0
1			l	1	0	0
			21
L1
			l31	0	1	0
			l41	0	0	1
Умножение L1-1 слева на P2		преобразует эту матрицу в

недиагональную. Вторая строка меняется местами с некоторой,

расположенной ниже. Умножение преобразованной матрицы справа на ту же P2 переставляет столбцы с теми же номерами, что и строки.

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

1

преобразование с помощью матрицы P3. Таким образом, L1 имеет требуемую форму и отличается от L1-1 только перестановкой элементов в

1

первом столбце. То же самое относится к L2 . Разложение имеет вид

PA=LU, где L= L1 L2 L3. Как и прежде не требуется инвертировать и

1 1

перемножать матрицы. Элементы из столбцов L1 , L2 и L31 просто переносятся с противоположными знаками в те же столбцы матрицы L.

Обобщим сказанное на случай матрицы размером n.

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

Ln-1-1 Pn-1 Ln-2-1 …Pi+1 Li-1 ….P3 L2-1P2L1-1P1A=U.

В результате определяется набор и последовательность перестановок. Далее, встраивая в нужные места соответствующие сомножители, приводим это разложение к виду

1	1	1	1	1
Ln 1	Ln 2 … Li		… L1	Pn-1…P1А=U, где Ln 1 = Ln11 ,

1

Li = Pn 1...Pi 1Li 1Pi 1...Pn 1 , (i=1,2,…,n-2).

134

Если при обработке очередного k – ого столбца перестановки не требуются, полагаем Pk=E.

В процессе разложения требуется определить три матрицы P, L и U.

Если эти матрицы определены, решение системы разбивается на следующие этапы:

1.Рассчитывают Pb.

2.Решают систему Ly=Pb.

3.Решают систему Ux=y.

Упражнения. Разработать программу на языке Matlab,

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

11.4.2.Метод Холецкого (метод квадратного корня).

Метод можно использовать в случае симметрической положительно определенной матрицы А. симметрическая матрица называется положительно определенной, если для любого вектора x выполнено условие xTAx>0. Симметрическая матрица имеет вещественные собственные числа λ. Она положительно определена, если все её собственные числа λ>0.

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

Разложение Холецкого имеет вид

135

					l11	0		0
					l21	l22		0
A LLT ; L										;
					ln1	ln 2			lnn
				l11	l21		ln1
					l22
T			0				ln2
L
			0	0			lnn

Если разложение выполнено, решение системы Ax=b осуществляется

вдва этапа:

1.Ly=b;

2.LTx=y.

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

Для нахождения элементов матрицы L элементы произведения LLT

приравнивают соответствующим элементам матрицы А. Матрица А – симметрическая, поэтому вычисляют только элементы на главной

диагонали и ниже её ( aij ;i j, n; j 1,n) . В результате получают n(n 1) 2

уравнений для того же количества неизвестных lij.

l	2	a
11		11
li1l11 ai1; i=2, n; элемeнты 1 -ого столбца
l	2	l 2	a	22
	21	22
li1l21 li 2l22				ai 2		; i=3, n; элемeнты 2 -ого столбца
.....................................
l	2	l2	...	l2		a
	k1	k 2			kk	kk
li1lk1 li 2lk 2					... lik lkk aik ; i=k 1, n элемeнты k -ого столбца
..............................................
l	2	l2	...	l2		a ; элемeнты n -ого столбца
	n1	n2			nn	nn

Эта система решается последовательно: сначала вычисляют l11 затем

li1, l22 и т. д.

136

l11
	a11
li1 ai1 / l11; i=2, n
l		a		22		l 2
22						21
li 2 (ai 2 li1l21 ) / l22								; i=3, n
.....................................
l a						l2	l 2	...	l 2	;k 1, n
kk			kk			k1	k 2		k ,k 1
lik (aik li1lk1 li 2lk 2									... li,k 1lk ,k 1 ) / lkk ; i=k 1, n
..............................................
l a				nn		l 2	l 2	...	l 2
nn						n1	n2		n,n 1

Подкоренные выражения положительны, если матрица А положительно определена.

Компактная форма записи тех же расчетных формул: 1. Для k 1, n проделать

	k 1
1.1 lkk	akk lki2 ;
	i 1

1.2 Для i=k 1, n проделать

k 1

aik lijlkj

	lik	j 1	.
		lkk

2. Вывод матрицы L.

Необходимое число операций Q≈n3/3, что вдвое меньше требуемого в методе Гаусса. Другим положительным качеством метода является высокая вычислительная устойчивость алгоритма.

Упражнения. Разработать программу на языке Matlab,

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

137

11.4.3.QR разложение матрицы и решение систем линейных алгебраических уравнений.

QR разложение означает факторизацию вида A=QR, где Q –

ортогональная матрица, R – верхняя треугольная матрица.

Матрица Q является ортогональной, если выполнено хотя бы одно

( а, следовательно, все) из следующих условий:

1.Q-1 = QT;

2.Cтолбцы Q образуют ортонормированную систему (т. е. (qi, qj)= 0

и(qi, qi ) = 1);

3.Cтроки Q, рассматриваемые как векторы, образуют

ортонормированную систему;

4. ||Qx|| = ||x|| для любого x.

Из второго и третьего условия следует, что если Q ортогональная матрица, то и QT ортогональная.

Из условия (4) следует, что ||Q||p=1, а из (1) и (4) – cond(Q)=1 (cond(A)=||A|| ||A-1||).

Отметим следующее свойство: произведение ортогональных матриц есть матрица ортогональная.

Если вычислено A =QR, система Ax = b решается следующим образом.

1.Решают систему Qy = b, y =QTb. При этом отпадает необходимость обращения матрицы Q.

2.Решают систему с верхней треугольной матрицей Rx = y.

1. Метод вращения.

Матрица вращения Tkl получается из единичной матрицы добавлением четырех элементов, расположенных на пересечении k-ой и l-ой строк с k-ым и l-ым столбцами (k<l)

138

		1
			1
			ckl	skl
Tkl
			skl	ckl
				1

Все элементы, кроме отмеченных, равны нулю. Матрица вращения ортогональна, т. е. ckl2+skl2=1, что эквивалентно skl=sin(φkl), ckl=cos(φkl).

Вектор y=Tklx отличается от вектора x элементами в k-ой и l-ой строках yk=cklxk+sklxl, yl=cklxl-sklxk. Действие Tkl на вектор x эквивалентно его повороту на угол φkl вокруг оси, перпендикулярной гиперплоскости Оxk xl . Для обнуления элемента xl требуется выполнение следующих условий: yl=cklxl-sklxk=0, ckl2+skl2=1, откуда следует

ckl

xk

, skl

xl

.

x2

x2

x2

x2

k

l

k

l

Преобразование исходной матрицы A к верхней треугольной R

осуществляется в результате последовательного обнуления элементов столбцов, расположенных ниже главной диагонали, начиная с первого столбца и заканчивая (n-1)-ым: A →A1 →…→Ak →…→An-1 =R . Первые k

столбцов матрицы Ak имеют нулевые компоненты ниже главной диагонали. На первом шаге вычисляют A1 , обнуляя (n-1) элементов ai1

(i=2,…, n) первого столбца матрицы А. Для исключения элемента a21

формируют матрицу вращения T12, полагая xk=a11, xl=a21, и вычисляют T12A=(T12a1,T12a2,…T12ai,…,T12an). Здесь ai i-ый столбец матрицы А. Для исключения очередного элемента ai1 первого столбца необходимо определить матрицу вращения T1i, положив xk=a11, xl=ai1. В результате получим A1=T1nT1,n-1…T13T12A=T1A. Матрица T1 ортогональна, так как является произведением ортогональных матриц. На k-ом шаге при

139

переходе Ak-1 →Ak после обнуления элементов k-ого столбца получают

Ak= Tk,nTk,n-1…Tk,k+2Tk,k+1Ak-1=TkAk-1. Выполнив (n-1) шагов (k=1,2,…,n-1),

получим Tn-1Tn-2…T1A=R или A=QR, где Q= T1T T2T…Tn-1T .

При формировании соответствующих матриц вращения учитывают номер обнуляемого элемента. Если, например, обнуляется элемент aij (i>j), то элементы матрицы Tji равны (xk=ajj, xl=aij).

c ji

a jj

, s ji

aij

.

a2

a2

a2

a2

jj

ij

jj

ij

Количество требуемых операций Q≈2n3. Метод обладает повышенной вычислительной устойчивостью в сравнении с методом Гаусса.

2.Метод отражения.

Матрица отражения (матрица Хаусхольдера) V=E-2wwT, w Rn,

	V Rn×n, ||w||=
				wT w			w12
	1 2w2					2w w
		1			1		2
		2w w			1 2w2
	V=
		2	1		2
		2wn w1				2wn w2

w22 ... wn2 =1.

2w1wn 2w2wn

1 2wn2

Матрица отражения симметрическая и ортогональная.

Действительно V=VT и VTV=(E-2wwT)(E-2wwT)=E-4wwT+4wwTwwT=E. То есть VT=V-1.

В основе использования матрицы отражения для QR разложения лежит следующая теорема

Теорема. Пусть заданы любые ненулевые векторы a и b. Существует единственный вектор w единичной длины такой, что определяемая им матрица отражения переводит вектор a в b (Va=b).

Действие матрицы V на вектор a есть зеркальное отражение вектора a относительно гиперплоскости, проходящей через начало координат и

перпендикулярной вектору w.

140