Автор cайта: Владимир Потемкин fortran-90@yandex.ru |
|
Собственные значения матриц
LA60. Приведение матрицы A к верхней почти треугольной форме Гессенберга с помощью ортогонального преобразования.
LA61. Вычисление собственных значений матрицы в форме Гессенберга с помощью QR-алгоритма Френсиса.
LA64. Вычисление собственных значений симметричной матрицы методом вращений (Якоби).
LA60
Программа LA60 производит преобразование матрицы А порядка N к верхней почти треугольной форме с помощью ортогонального преобразования. При этом собственные значения матрицы не меняются.
Алгоритм преобразования заключается в последовательном умножении исходной матрицы на ортогональные матрицы Q(k)
As+1 = QsAsQs
для s = 1, 2, ..., N-2. Здесь ортогональная матрица Q представляет собой матрицу отражения Qs = E - wswsT/γ,
где wsT = (as+1,s+α, as+2,s, ..., an,s) - часть столбца матрицы А, преобразуемая на s-ом шаге; α = Sign(as+1,s)•SQRT(as+1,s2 + as+2,s2 + ... + an,s2); γ = α(as+1,s + α). При таком преобразовании в исходной матрице последовательно обнуляется нижняя часть столбцов с элементами as+2,s, ..., an,s, а элемент as+1,s становится равным -α и этим самым матрица А за N-2 шагов приводится к верхней почти треугольной форме Гессенберга.
Каждый s-ый шаг алгоритма разделяется на два этапа. Сначала исходная матрица As умножается на матрицу Qs слева
as+1,s = as+1,s + α;
n
vj = ∑ aisaij;
i=s+1
ai,j = ai,j - ai,svj/γ (i = s+1, s+2, ..., n);
для всех j от s+1 до n. На втором этапе получившаяся промежуточная матрица умножается на матрицу Qs справа
n
vi = ∑ aijajs;
j=s+1
ai,j = ai,j - ai,svi/γ (j = s+1, s+2, ..., n);
для всех i от 1 до n. Вычисления повторяются для s=1, 2, ..., n-2.
Вызов программы
CALL LA60(A, H)
Параметры программы
A - входной параметр;
H - выходной параметр.
Real A(1:N, 1:N) - массив, содержащий исходную матрицу;
Real H(1:N, 1:N) - массив, содержащий матрицу в форме Гессенберга. При вызове матрицы А и H можно совмещать, тогда матрица H будет записана вместо матрицы А.
Вернуться к оглавлению Скачать LA60
LA61
Программа LA61 вычисляет собственные значения верхней почти треугольной матрицы Гессенберга, используя двойной QR алгоритм Френсиса с дополнительным сдвигом.
Произвольную квадратную матрицу А можно представить в виде произведения ортогональной Q (Q-1 = АT) и правой треугольной матрицы R
A = QR. (1)
Отсюда имеем, что R = QTA. Домножая на матрицу Q справа, получаем, что RQ = QTAQ = A1, где A1 - матрица, подобная матрице А. Строим теперь последовательность матриц As по следующему правилу. Матрицу As разлагаем в произведение ортогональной и правой треугольной матриц As = QsRs и полагаем As+1 = RsQs = QsTAsQs, причем матрицы из последовательности As подобны между собой и подобны исходной матрице A.
Построенная последовательность матриц при определенных условиях сходится к правой диагональной матрице, диагональные элементы которой являются собственными значениями исходной матрицы. Исключение в смысле сходимости представляют отдельные поддиагональные элементы, каждый из которых связан с матрицей второго порядка на диагонали, которая имеет комплексно-сопряженную пару собственных значений.
Ускорения сходимости можно достигнуть, если применить вариант QR-алгоритма со сдвигом, который заключается в следующем. Последовательность матриц строится исходя из соотношений
A0 - k0E = Q0R0,
A1 = R0Q0 + k0E = Q0T(A0 - k0E)Q0 + k0E = Q0TA0Q0;
. . . . . . . . . . . . . . . . . . . . . . . . .
As - ksE = QsRs,
As+1 = RsQs + ksE = QsTAsQs, (2)
где ks - сдвиг, который для увеличения скорости сходимости должен быть близок к одному из собственных значений матрицы А (и, следовательно, всех Аs).
Однако даже когда А - вещественная матрица, некоторые из ее собственных значений могут быть комплексными. Тогда в преобразовании (2) необходимо использовать комплексное значение ks и матрица Аs+1 становится комплексной. Эту трудность можно преодолеть, выполняя одновременно два шага вида (2) со сдвигами ks и ks*. Тогда при точном вычислении матрица Аp+2 остается вещественной. Выполняя два шага в соответствии с (2) получаем следующие соотношения
Ap+2 = Qs+1TQsTAsQsQs+1. (3)
Поскольку QsRs = As - ksE и Qs+1Rs+1 = As+1 - ks+1E = QsTAsQs - ks+1E, получаем, что
QsQs+1Rs+1QsT = As - ks+1E и
QsQs+1Rs+1Rs = (As - ks+1E)(As - ksE).
Приняв QsQs+1 = Q, Rs+1Rs = R, (As - ks+1E)(As - ksE) = M, получим
QR = M, (4)
Ap+2 = QTAsQ. (5)
Отсюда следует, что R = QTM, т.е. что матрица QT приводит произведение матриц, обозначенное через M, к треугольному виду. Это произведение является вещественной матрицей, если ks и ks+1 вещественны или комплексно сопряжены. Следовательно, QT и Q - вещественные матрицы. Любую вещественную матрицу можно привести к треугольному виду, последовательно умножая слева на вещественные симметричные и ортогональные матрицы отражения U1, U2, ..., Un-1
Us = E - 2wswsT, (6)
где ws - единичный вектор, у которого первые s-1 компонент равны нулю. Поэтому QT = Un-1TUn-2T ... U2TU1T и Q = Un-1Un-2 ... U2U1.
Первый столбец матрицы Q совпадает с первым столбцом матрицы U1, следовательно, остается определить матрицу Q, удовлетворяющую соотношению (4) и имеющую своим первым столбцом первый столбец матрицы U1. После этого матрицу Ap+2 легко найти из соотношения (5).
Для нахождения матрицы Q может быть использовано следующее свойство: если матрица B неособая и выполнено условие
H = QTBQ,
где Q - ортогональная матрица, Q - правая почти треугольная матрица Гессенберга, то полные матрицы Q и H определяются первым столбцом матрицы Q и это представление единственно, если матрица H имеет положительные поддиагональные элементы.
Теперь предположим, что перед применением очередного двойного шага QR-алгоритма матрица Аs каким-либо образом преобразована к правому почти треугольному виду Гессенберга. В этом случае для вычисления матрицы Аp+2 нет необходимости производить вычисления матрицы М. В место этого для получения Аp+2 достаточно знать лишь первый столбец матрицы М, в котором отличны от нуля лишь первые три элемента, равные
│ p1 = a112 + a12a21 - a11(ks + ks+1) + ksks+1
│ q1 = a21(a11 + a22 - ks - ks+1) (7)
│ r1 = a21a32
где aij - элементы матрицы As. Далее вычисляем матрицу U1T, которая обнуляет q1 и r1, причем в ее первом столбце (и первой строке) также отличны от нуля лишь первые три элемента. Тогда матрица U1TAsU1 для n=8 имеет вид
│ × × × × × × × × │
│ × × × × × × × × │
│ + × × × × × × × │
│ + + × × × × × × │
│ × × × × × │
│ × × × × │
│ × × × │
│ × × │
т.е. является матрицей Гессенберга с тремя дополнительными элементами, обозначенными знаком "+".
Любая вещественная матрица K при помощи ортогонального подобного преобразования может быть сведена к матрице Гессенберга
H = Un-1T ... U2TKU2 ... Un-1,
где Us определяются из соотношения (6). В частности, можно считать матрицу K равной U1TAsU1. Однако первый столбец матрицы U1U2 ... Un-1 есть первый столбец матрицы U1 и тем самым искомая матрица Q равна
Q = U1U2 ... Un-1,
а матрица Ap+2 равна матрице H.
На каждом шаге преобразования Хаусхольдера (приведения матрицы к форме Гессенберга) текущая матрица имеет вид
│ × × × × × × × × │
│ × × × × × × × × │
│ × × × × × × × │
│ × × × × × × │
│ + × × × × × │
│ + + × × × × │
│ × × × │
│ × × │
Эта матрица возникает перед преобразованием U4 для матрицы 8-го порядка и является матрицей Гессенберга с дополнительными элементами, обозначенными знаком "+". Таким образом матрица Us имеет вектор ws вида
(0, 0, ..., 0, ps, qs, rs, 0, ..., 0),
где ненулевые элементы ps, qs и rs определяются элементами (s, s-1), (s+1, S-1), (s+2, S-1) текущей матрицы. Очевидно, что и первый шаг, состоящий в левом и правом умножении на матрицу U1, имеет такой же вид, как и при преобразовании Хаусхльдера, поскольку матрица U1 определяется также только тремя элементами из соотношений (7).
Работа программы.
Преобразование Хаусхольдера на каждом шаге определяется тремя элементами ps, qs и rs. На первом шаге они вычисляются согласно соотношению (7), на всех последующих используют элементы текущей матрицы с номерами (s, s-1), (s+1, S-1), (s+2, S-1). Во избежание переполнения разрядной сетки величины ps, qs и rs нормируются путем деления их на величину (|ps| + |qs| + |rs|), т.к. преобразование Хаусхольдера зависит лишь от отношения величин ps, qs, rs.
Число арифметических операций при вычислениях с матрицей E - 2wswsT можно значительно уменьшить, если выражение 2wswsT представить в виде usvsT, где векторы us и vsT параллельны ws, а их элементы, не равные нулю, определяются из соотношений
(ps ± as), qs/as, rs/as или 1, qs/(ps ± as), rs/(ps ± as), где as2 = ps2 + qs2 + rs2.
Перед каждой итерацией матрица Hs проверяется на малость поддиагональных элементов. Если последний малый элемент есть (l, l-1), то это позволяет расщепить исходную матрицу на две матрицы более низкого норядка. В этом случае итерации продолжаются с матрицей, расположенной в нижнем правом углу с номерами строк и столбцов от l до n. Если ни один поддиагональный элемент не является малым, то l принимается равным 1. Критерий малости элементов, используемый в программе
|hl,l-1| <= MashEpsilon*(|hl-1,l-1 + hl,l|),
который сравнивает на малость элемент hl,l-1 с ближайшими диагональными элементами.
Сдвиги в каждой итерации зависят от собственных значений матрицы 2×2, расположенной в правом нижнем углу текущей матрицы Hs и их можно вычислить из системы уравнений
│ks + ks+1 = hn-1,n-1(s) + hnn(s)
│ksks+1 = hn-1,n-1(s)hn-1,n-1(s)hnn(s) - hn-1,n(s)hn,n-1(s)
При использовании указанных сдвигов после некоторого количества итераций возникает ситуация, когда либо элемент hn-1,n-1(s), либо элемент hn-1,n-2(s) становится принебрежимо малым. Если выполняется первое условие, то элемент hnn(s) можно считать собственным значением и порядок матрицы можно уменьшить на 1 вычеркиванием последней строки и последнего столбца. Если выполнено второе условие, т.е. элемент hn-1,n-2(s) пренебрежимо мал, то можно выделить пару собственных значений из матрицы 2×2, расположенной в правом нижнем углу, а порядок исходной матрицы снижается на 2 вычеркиванием двух последних строк и двух последних столбцов.
Итерационный процесс может прерваться, если достигнуто максимальное количество итераций или если один из двух последних поддиагональных элементов не изменился после проведенной итерации. В этом случае наименьший из двух последних поддиагональных элементов принимается равным нулю и выделяется одно или пара собственных значений, в зависимости от того, какой из элементов оказался наименьшим.
Вызов программы
CALL LA61(H, WR, WI, Cnt, Error)
Параметры программы
H - входной и выходной параметр;
WR, WI, Cnt, Error - выходные параметры.
Real H(1:N, 1:N) - массив, содержащий исходную матрицу в форме Гессенберга. В процессе вычислений массив видоизменяется;
Real WR(1:N), WI(1:N) - массивы, содержащие вещественную и мнимую части собственных чисел;
Integer Cnt(1:N) - массив, показывающий, сколько потребовалось итераций для вычисления каждого собственного значения. Если одновременно найдено два собственных значения, то число итераций выдается со знаком плюс для первого собственного значения и со знаком минус для второго.
Integer Error - индикатор ошибки. Error=0, если ошибок нет, и равно 65, если после 30 итераций какое-либо из собственных значений не было найдено.
Пример
!Вычисление собственных значений вещественной матрицы.
program TestLA61
use NML
implicit none
integer, parameter:: n=4
real:: A(n,n), H(n,n), WR(n), WI(n)
integer:: Cnt(n), Er
data A/1.0, 1.1, 1.2, 1.4, 1.1, 1.1, 1.2, 1.3, 1.2, 1.2, 1.2, 1.3, 1.4, 1.3, 1.3, 1.3/
integer i, j
!begin
H=A
call LA60(A, H)
print 10
print 11, ((A(i,j), j=1,n), i=1,n)
print 12
print 13, ((H(i,j), j=1,n), i=1,n)
call LA61(H, WR, WI, Cnt, Er)
if (Er==0) then
print 14
print 15, (WR(i), WI(i), Cnt(i), i=1,n)
else
print 16
end if
10 format(/' ?aa®¤ i ¬ aa?? :'/)
11 format(8x,<n>F6.1)
12 format(/' ? aa?? ? a®a¬? ??aa???a? :'/)
13 format(1x,<n>F13.6)
14 format(/' ‘®?aa??e? § c??i:'/)
15 format(8x,2F15.6,I6)
16 format(/' ?a®?a ¬¬ ®aa ®?«? .'/)
end program TestLA61
Исходная матрица:
1.0 1.1 1.2 1.4
1.1 1.1 1.2 1.3
1.2 1.2 1.2 1.3
1.4 1.3 1.3 1.3
Матрица в форме Гессенберга:
1.000000 -2.147091 0.000000 0.000000
-2.147091 3.719523 -0.261293 0.000000
0.000000 -0.261293 -0.083925 -0.012079
0.000000 0.000000 -0.012079 -0.035598
Собственные значения:
4.911706 0.000000 0
-0.271466 0.000000 0
-0.001959 0.000000 3
-0.038279 0.000000 -3
Вернуться к оглавлению Скачать LA61
LA64
Программа LA64 вычисляет собственные значения симметричной матрицы методом вращений (методом Якоби).
Произвольную симметричную матрицу A можно представить в виде произведения QTAQ = D, где Q - ортогональная (QT = Q-1), а D - диагональная матрицы, причем матрица D подобна матрице A, т.е. имеет те же собственные значения, что и матрица A. Поскольку собственными значениями диагональной матрицы являются ее диагональные элементы, то зная матрицу Q можно найти все собственные значения матрицы A.
Матрица Q получается как предел последовательности произведений матриц простых поворотов Q = V1...Vp..., где Vp имеет вид (для k<l)
│ 1 │
│ . │
│ 1 │
│ cosθ . . sinθ │ k
│ 1 │
│ . │
│ . │
│ 1 │
│ -sinθ . . cosθ │ l
│ 1 │
│ . │
│ 1 │
k l
где все невыписанные элементы - нулевые, а θ - угол поворота. Матрица D получается как предел последовательности матриц Ap, удовлетворяющих соотношению
Ap = VpTAp-1Vp (p = 0, 1, 2, ...), где A0 = A.
При преобразовании матрицы Ap-1 с помощью матрицы простого поворота на каждом шаге аннулируется внедиагональный элемент akl(p-1) (k≠l). Из условия akl(p) = 0, где akl(p) - элемент преобразованной матрицы, получаем соотношение
tg2θ = 2akl(p-1)/(all(p-1) - akk(p-1)), где θ всегда берется в пределах |2θ| ≤ π/2.
При таком выборе θ имеем
Матрица Ap-1 преобразуется по формулам
aij(p) = aij(p-1) (i,j ≠ k,l);
aik(p) = aik(p-1)cosθ - ail(p-1)sinθ (i ≠ k,l);
ail(p) = aik(p-1)sinθ + ail(p-1)cosθ (i ≠ k,l);
akk(p) = akk(p-1)cos2θ + all(p-1)sin2θ - 2akl(p-1)cosθsinθ;
all(p) = akk(p-1)sin2θ + all(p-1)cos2θ + 2akl(p-1)cosθsinθ;
akl(p) = (akk(p-1) - all(p-1))cosθsinθ + akl(p-1)(cos2θ - sin2θ).
Внедиагональные элементы аннулируются циклически по столбцам в соответствии со следующей нумерацией пар (k,l)
(1,2), (1,3), (2,3), (1,4), (2,4), (3,4), ..., (1,n), ..., (n-1,n).
Для того, чтобы не аннулировать малые внедиагональные элементы, в то время как присутствуют элементы, большие по значению, осуществляется проверка условия
|akl(p-1)| < σ/n,
где σ вначале выбирается раной σ0 = (2 ∑ akl2)½ для k < l и ликвидируются все элементы, большие σ0, после чего значение σ0 уменьшается в n раз. Процесс продолжается до тех пор, когда все внедиагональные элементы по модулю не станут меньше, чем величина σ1 = εσ0/n, где ε - заданная погрешность решения.
Вызов программы
CALL LA64(A, Eps)
Параметры программы
A - входной и выходной параметр;
Eps - входной параметр;
Real A(1:N, 1:N) - массив, содержащий исходную симметричную матрицу. На выходе - на главной диагонали располагаются собственные значения матрицы A. Элементы матрицы выше главной диагонали используются при вычислениях, элементы ниже главной диагонали не используются и не изменяются;
Real Eps - погрешность решения.
Пример
!Вычисление собственных значений вещественной симметричной матрицы.
program TestLA64
use NML_LA64
implicit none
integer, parameter:: n=4
real:: A(n,n)
real, parameter:: Eps=1.E-8
data A/1.0, 1.1, 1.2, 1.4, 1.1, 1.1, 1.2, 1.3, 1.2, 1.2, 1.2, 1.3, 1.4, 1.3, 1.3, 1.3/
integer i, j
!begin
print 10
print 11, ((A(i,j), j=1,n), i=1,n)
call LA64(A,Eps)
print 12
print 13, (A(i,i), i=1,n)
10 format(/' ?aa®¤ i ¬ aa?? :'/)
11 format(8x,<n>F6.1)
12 format(/' ‘®?aa??e? § c??i:'/)
13 format(1x,<n>F13.6)
end
Исходная матрица:
1.0 1.1 1.2 1.4
1.1 1.1 1.2 1.3
1.2 1.2 1.2 1.3
1.4 1.3 1.3 1.3
Собственные значения:
-0.001959 4.911704 -0.038279 -0.271466
Вернуться к оглавлению Скачать LA64
|