Собственные значения матриц. Алгоритмы Френсиса и Якоби. Библиотека численных методов на языке Fortran 90

LA60

Программа LA60 производит преобразование матрицы А порядка N к верхней почти треугольной форме с помощью ортогонального преобразования. При этом собственные значения матрицы не меняются.
Алгоритм преобразования заключается в последовательном умножении исходной матрицы на ортогональные матрицы Q^(k)
A_s+1 = Q_sA_sQ_s
для s = 1, 2, ..., N-2. Здесь ортогональная матрица Q представляет собой матрицу отражения
Q_s = E - w_sw_s^T/γ,
где w_s^T = (a_s+1,s+α, a_s+2,s, ..., a_n,s) - часть столбца матрицы А, преобразуемая на s-ом шаге; α = Sign(a_s+1,s)•SQRT(a_s+1,s² + a_s+2,s² + ... + a_n,s²); γ = α(a_s+1,s + α). При таком преобразовании в исходной матрице последовательно обнуляется нижняя часть столбцов с элементами a_s+2,s, ..., a_n,s, а элемент a_s+1,s становится равным -α и этим самым матрица А за N-2 шагов приводится к верхней почти треугольной форме Гессенберга.
Каждый s-ый шаг алгоритма разделяется на два этапа. Сначала исходная матрица A_s умножается на матрицу Q_s слева
a_s+1,s = a_s+1,s + α;
_n
v_j = ∑ a_isa_ij;
^i=s+1
a_i,j = a_i,j - a_i,sv_j/γ (i = s+1, s+2, ..., n);
для всех j от s+1 до n. На втором этапе получившаяся промежуточная матрица умножается на матрицу Q_s справа
_n
v_i = ∑ a_ija_js;
^j=s+1
a_i,j = a_i,j - a_i,sv_i/γ (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 = Q^TA. Домножая на матрицу Q справа, получаем, что RQ = Q^TAQ = A₁, где A₁ - матрица, подобная матрице А. Строим теперь последовательность матриц A_s по следующему правилу. Матрицу A_s разлагаем в произведение ортогональной и правой треугольной матриц A_s = Q_sR_s и полагаем A_s+1 = R_sQ_s = Q_s^TA_sQ_s, причем матрицы из последовательности A_s подобны между собой и подобны исходной матрице A.
   Построенная последовательность матриц при определенных условиях сходится к правой диагональной матрице, диагональные элементы которой являются собственными значениями исходной матрицы. Исключение в смысле сходимости представляют отдельные поддиагональные элементы, каждый из которых связан с матрицей второго порядка на диагонали, которая имеет комплексно-сопряженную пару собственных значений.
   Ускорения сходимости можно достигнуть, если применить вариант QR-алгоритма со сдвигом, который заключается в следующем. Последовательность матриц строится исходя из соотношений
A₀ - k₀E = Q₀R₀, A₁ = R₀Q₀ + k₀E = Q₀^T(A₀ - k₀E)Q₀ + k₀E = Q₀^TA₀Q₀;
. . . . . . . . . . . . . . . . . . . . . . . . .
A_s - k_sE = Q_sR_s, A_s+1 = R_sQ_s + k_sE = Q_s^TA_sQ_s, (2)
где k_s - сдвиг, который для увеличения скорости сходимости должен быть близок к одному из собственных значений матрицы А (и, следовательно, всех А_s).
   Однако даже когда А - вещественная матрица, некоторые из ее собственных значений могут быть комплексными. Тогда в преобразовании (2) необходимо использовать комплексное значение k_s и матрица А_s+1 становится комплексной. Эту трудность можно преодолеть, выполняя одновременно два шага вида (2) со сдвигами k_s и k_s^*. Тогда при точном вычислении матрица А_p+2 остается вещественной. Выполняя два шага в соответствии с (2) получаем следующие соотношения
A_p+2 = Q_s+1^TQ_s^TA_sQ_sQ_s+1. (3)
Поскольку Q_sR_s = A_s - k_sE и Q_s+1R_s+1 = A_s+1 - k_s+1E = Q_s^TA_sQ_s - k_s+1E, получаем, что
Q_sQ_s+1R_s+1Q_s^T = A_s - k_s+1E и
Q_sQ_s+1R_s+1R_s = (A_s - k_s+1E)(A_s - k_sE).
Приняв Q_sQ_s+1 = Q, R_s+1R_s = R, (A_s - k_s+1E)(A_s - k_sE) = M, получим
QR = M, (4)
A_p+2 = Q^TA_sQ. (5)
Отсюда следует, что R = Q^TM, т.е. что матрица Q^T приводит произведение матриц, обозначенное через M, к треугольному виду. Это произведение является вещественной матрицей, если k_s и k_s+1 вещественны или комплексно сопряжены. Следовательно, Q^T и Q - вещественные матрицы. Любую вещественную матрицу можно привести к треугольному виду, последовательно умножая слева на вещественные симметричные и ортогональные матрицы отражения U₁, U₂, ..., U_n-1
U_s = E - 2w_sw_s^T, (6)
где w_s - единичный вектор, у которого первые s-1 компонент равны нулю. Поэтому Q^T = U_n-1^TU_n-2^T ... U₂^TU₁^T и Q = U_n-1U_n-2 ... U₂U₁.
Первый столбец матрицы Q совпадает с первым столбцом матрицы U₁, следовательно, остается определить матрицу Q, удовлетворяющую соотношению (4) и имеющую своим первым столбцом первый столбец матрицы U₁. После этого матрицу A_p+2 легко найти из соотношения (5).
Для нахождения матрицы Q может быть использовано следующее свойство: если матрица B неособая и выполнено условие
H = Q^TBQ,
где Q - ортогональная матрица, Q - правая почти треугольная матрица Гессенберга, то полные матрицы Q и H определяются первым столбцом матрицы Q и это представление единственно, если матрица H имеет положительные поддиагональные элементы.
Теперь предположим, что перед применением очередного двойного шага QR-алгоритма матрица А_s каким-либо образом преобразована к правому почти треугольному виду Гессенберга. В этом случае для вычисления матрицы А_p+2 нет необходимости производить вычисления матрицы М. В место этого для получения А_p+2 достаточно знать лишь первый столбец матрицы М, в котором отличны от нуля лишь первые три элемента, равные
│ p₁ = a₁₁² + a₁₂a₂₁ - a₁₁(k_s + k_s+1) + k_sk_s+1
│ q₁ = a₂₁(a₁₁ + a₂₂ - k_s - k_s+1) (7)
│ r₁ = a₂₁a₃₂
где a_ij - элементы матрицы A_s. Далее вычисляем матрицу U₁^T, которая обнуляет q₁ и r₁, причем в ее первом столбце (и первой строке) также отличны от нуля лишь первые три элемента. Тогда матрица U₁^TA_sU₁ для n=8 имеет вид
│ × × × × × × × × │
│ × × × × × × × × │
│ + × × × × × × × │
│ + + × × × × × × │
│ × × × × × │
│ × × × × │
│ × × × │
│ × × │
т.е. является матрицей Гессенберга с тремя дополнительными элементами, обозначенными знаком "+".
Любая вещественная матрица K при помощи ортогонального подобного преобразования может быть сведена к матрице Гессенберга
H = U_n-1^T ... U₂^TKU₂ ... U_n-1,
где U_s определяются из соотношения (6). В частности, можно считать матрицу K равной U₁^TA_sU₁. Однако первый столбец матрицы U₁U₂ ... U_n-1 есть первый столбец матрицы U₁ и тем самым искомая матрица Q равна
Q = U₁U₂ ... U_n-1,
а матрица A_p+2 равна матрице H.
На каждом шаге преобразования Хаусхольдера (приведения матрицы к форме Гессенберга) текущая матрица имеет вид
│ × × × × × × × × │
│ × × × × × × × × │
│ × × × × × × × │
│ × × × × × × │
│ + × × × × × │
│ + + × × × × │
│ × × × │
│ × × │
Эта матрица возникает перед преобразованием U₄ для матрицы 8-го порядка и является матрицей Гессенберга с дополнительными элементами, обозначенными знаком "+". Таким образом матрица U_s имеет вектор w_s вида
(0, 0, ..., 0, p_s, q_s, r_s, 0, ..., 0),
где ненулевые элементы p_s, q_s и r_s определяются элементами (s, s-1), (s+1, S-1), (s+2, S-1) текущей матрицы. Очевидно, что и первый шаг, состоящий в левом и правом умножении на матрицу U₁, имеет такой же вид, как и при преобразовании Хаусхльдера, поскольку матрица U₁ определяется также только тремя элементами из соотношений (7).
Работа программы.
Преобразование Хаусхольдера на каждом шаге определяется тремя элементами p_s, q_s и r_s. На первом шаге они вычисляются согласно соотношению (7), на всех последующих используют элементы текущей матрицы с номерами (s, s-1), (s+1, S-1), (s+2, S-1). Во избежание переполнения разрядной сетки величины p_s, q_s и r_s нормируются путем деления их на величину (|p_s| + |q_s| + |r_s|), т.к. преобразование Хаусхольдера зависит лишь от отношения величин p_s, q_s, r_s.
Число арифметических операций при вычислениях с матрицей E - 2w_sw_s^T можно значительно уменьшить, если выражение 2w_sw_s^T представить в виде u_sv_s^T, где векторы u_s и v_s^T параллельны w_s, а их элементы, не равные нулю, определяются из соотношений
(p_s ± a_s), q_s/a_s, r_s/a_s или 1, q_s/(p_s ± a_s), r_s/(p_s ± a_s),
где a_s² = p_s² + q_s² + r_s².
Перед каждой итерацией матрица H_s проверяется на малость поддиагональных элементов. Если последний малый элемент есть (l, l-1), то это позволяет расщепить исходную матрицу на две матрицы более низкого норядка. В этом случае итерации продолжаются с матрицей, расположенной в нижнем правом углу с номерами строк и столбцов от l до n. Если ни один поддиагональный элемент не является малым, то l принимается равным 1. Критерий малости элементов, используемый в программе
|h_l,l-1| <= MashEpsilon*(|h_l-1,l-1 + h_l,l|),
который сравнивает на малость элемент h_l,l-1 с ближайшими диагональными элементами.
Сдвиги в каждой итерации зависят от собственных значений матрицы 2×2, расположенной в правом нижнем углу текущей матрицы H_s и их можно вычислить из системы уравнений
│k_s + k_s+1 = h_n-1,n-1^(s) + h_nn^(s)
│k_sk_s+1 = h_n-1,n-1^(s)h_n-1,n-1^(s)h_nn^(s) - h_n-1,n^(s)h_n,n-1^(s)
При использовании указанных сдвигов после некоторого количества итераций возникает ситуация, когда либо элемент h_n-1,n-1^(s), либо элемент h_n-1,n-2^(s) становится принебрежимо малым. Если выполняется первое условие, то элемент h_nn^(s) можно считать собственным значением и порядок матрицы можно уменьшить на 1 вычеркиванием последней строки и последнего столбца. Если выполнено второе условие, т.е. элемент h_n-1,n-2^(s) пренебрежимо мал, то можно выделить пару собственных значений из матрицы 2×2, расположенной в правом нижнем углу, а порядок исходной матрицы снижается на 2 вычеркиванием двух последних строк и двух последних столбцов.
Итерационный процесс может прерваться, если достигнуто максимальное количество итераций или если один из двух последних поддиагональных элементов не изменился после проведенной итерации. В этом случае наименьший из двух последних поддиагональных элементов принимается равным нулю и выделяется одно или пара собственных значений, в зависимости от того, какой из элементов оказался наименьшим.
Перед вызовом программы LA61 нужно выполнить программу LA60.

Вызов программы

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 можно представить в виде произведения
Q^TAQ = D,
где Q - ортогональная (Q^T = Q^-1), а D - диагональная матрицы, причем матрица D подобна матрице A, т.е. имеет те же собственные значения, что и матрица A. Поскольку собственными значениями диагональной матрицы являются ее диагональные элементы, то зная матрицу Q можно найти все собственные значения матрицы A.
Матрица Q получается как предел последовательности произведений матриц простых поворотов Q = V₁...V_p..., где V_p имеет вид (для k<l)
│ 1 │
│ . │
│ 1 │
│ cosθ . . sinθ │    k
│ 1 │
│ . │
│ . │
│ 1 │
│ -sinθ . . cosθ │    l
│ 1 │
│ . │
│ 1 │
k l
где все невыписанные элементы - нулевые, а θ - угол поворота. Матрица D получается как предел последовательности матриц A_p, удовлетворяющих соотношению
A_p = V_p^TA_p-1V_p    (p = 0, 1, 2, ...),
где A₀ = A.
При преобразовании матрицы A_p-1 с помощью матрицы простого поворота на каждом шаге аннулируется внедиагональный элемент a_kl^(p-1) (k≠l). Из условия a_kl^(p) = 0, где a_kl^(p) - элемент преобразованной матрицы, получаем соотношение
tg2θ = 2a_kl^(p-1)/(a_ll^(p-1) - a_kk^(p-1)),
где θ всегда берется в пределах |2θ| ≤ π/2.
При таком выборе θ имеем

formula
Матрица A_p-1 преобразуется по формулам
a_ij^(p) = a_ij^(p-1)    (i,j ≠ k,l);
a_ik^(p) = a_ik^(p-1)cosθ - a_il^(p-1)sinθ    (i ≠ k,l);
a_il^(p) = a_ik^(p-1)sinθ + a_il^(p-1)cosθ    (i ≠ k,l);
a_kk^(p) = a_kk^(p-1)cos²θ + a_ll^(p-1)sin²θ - 2a_kl^(p-1)cosθsinθ;
a_ll^(p) = a_kk^(p-1)sin²θ + a_ll^(p-1)cos²θ + 2a_kl^(p-1)cosθsinθ;
a_kl^(p) = (a_kk^(p-1) - a_ll^(p-1))cosθsinθ + a_kl^(p-1)(cos²θ - sin²θ).
Внедиагональные элементы аннулируются циклически по столбцам в соответствии со следующей нумерацией пар (k,l)
(1,2), (1,3), (2,3), (1,4), (2,4), (3,4), ..., (1,n), ..., (n-1,n).
Для того, чтобы не аннулировать малые внедиагональные элементы, в то время как присутствуют элементы, большие по значению, осуществляется проверка условия
|a_kl^(p-1)| < σ/n,
где σ вначале выбирается раной σ₀ = (2 ∑ a_kl²)^½ для k < l и ликвидируются все элементы, большие σ₀, после чего значение σ₀ уменьшается в n раз. Процесс продолжается до тех пор, когда все внедиагональные элементы по модулю не станут меньше, чем величина σ₁ = εσ₀/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