100 FORMAT(1Х,'ТНЕ EIGENVALUES ARE'')
101 FORMAT(1X,E15.8)
STOP
END
Результат работы программы получаем в виде:
Собственные значения равны
0.33709179E 08
0.19149061E 08
0.71417603E 07
Метод Гивенса для симметричных матриц
Метод Гивенса основан на преобразовании подобия, аналогичном применяемому в методе Якоби. Однако в этом случае алгоритм построен таким образом, что вновь образованные нулевые элементы при всех последующих преобразованиях сохраняются. Поэтому метод Гивенса требует выполнения конечного числа преобразований и по сравнению с методом Якоби связан с меньшими затратами машинного времени. Его единственный недостаток состоит в том, что симметричная матрица приводится не к диагональному, а к трехдиагональному виду. Ниже будет показано, что такая форма матрицы может быть весьма полезной и оправдывает усилия, затраченные на ее получение.
В случае матрицы размерности п х п метод Гивенса требует п — 2 основных шагов, на каждом из которых выполняется ряд преобразований, число которых зависит от числа нулей, которое хотят получить в данном столбце или строке. На k-м шаге обращают в нули элементы, стоящие вне трех диагоналей k-й строки и k -го столбца, сохраняя в то же время нулевые элементы, полученные на предыдущих шагах. Таким образом, перед началом k -го шага преобразованная матрица является трехдиагональной, если ограничиться рассмотрением ее первых k — 1 строк и столбцов. По мере преобразований симметричная матрица размерности 5х5 приобретает следующие формы:
* | * | * | * | * | ||
* | * | * | * | * | ||
A0= | * | * | * | * | * | исходная матрица, |
* | * | * | * | * | ||
* | * | * | * | * |
* | * | 0 | 0 | 0 | ||
* | * | * | * | * | ||
A1= | 0 | * | * | * | * | после первого основного шага, |
0 | * | * | * | * | состоящего из трех преобразований, | |
0 | * | * | * | * |
* | * | 0 | 0 | 0 | ||
* | * | * | 0 | 0 | ||
A2= | 0 | * | * | * | * | после второго основного шага, |
0 | 0 | * | * | * | состоящего из двух преобразований, | |
0 | 0 | * | * | * |
* | * | 0 | 0 | 0 | ||
* | * | * | 0 | 0 | после третьего основного шага, | |
A3= | 0 | * | * | * | 0 | состоящего из одного преобразования. |
0 | 0 | * | * | * | Теперь матрица имеет трехдиагональный вид. | |
0 | 0 | 0 | * | * |
На каждом основном шаге изменяются лишь те элементы матрицы аij, которые расположены в ее правой нижней (заштрихованной) части. Таким образом на k-м шаге преобразуется только матрица порядка (п —k + 1), занимающая правый нижний угол исходной матрицы. Ясно, что на каждой следующей стадии выполняется меньшее число преобразований, чем на предыдущей. Всего для приведения матрицы к трехдиагональному виду требуется выполнить (n2 — Зп + 2)/2 преобразований.
Наш опыт применения метода Гивенса показывает, что можно при выполнении одного шага преобразований обратить в нуль сразу все элементы целой строки и столбца, стоящие вне трех диагоналей матрицы. Метод, позволяющий выполнить такое преобразование, предложил Хаусхолдер .
Метод Хаусхолдера для симметричных матриц
Метод Хаусхолдера позволяет привести матрицу к трехдиагональному виду, выполнив почти вдвое меньше вычислений по сравнению с другими методами. Это обусловлено тем, что при его применении становятся нулевыми сразу все элементы строк и столбцов, стоящие вне трех диагоналей матрицы. Метод Хаусхолдера позволяет получить требуемый результат быстрее, чем метод Гивенса, так как связан с выполнением меньшего числа, хотя и более сложных преобразований. Это его свойство особенно ярко проявляется применительно к большим матрицам. Хотя в методе Хаусхолдера вместо плоских вращении используются эрмитовыортогональные преобразования матриц, трехдиагональная форма матрицы, которую получают этим методом, имеет те же собственные значения, что и трехдиагональная матрица, получаемая методом Гивенса. При использовании метода Хаусхолдера на п — 2 основных шагах выполняются следующие преобразования:
Аk = РkAk-1Рk, k=1, 2, ...,п-2,
где Aо == А.
Каждая преобразующая матрица имеет вид
uk ukT
Pk = E - -------------- ,
2Kk2
где
ui,k = 0 при i = 1, 2, …, k,
ui,k = ak,i при i = k+2, …, n,
uk+1,k = ak,k+1 ± Sk.
Здесь
n 1/2Sk = Sa2k,i
i=k+1
2K2k = S2k ± ak, k+1 Sk.
В этих уравнениях берется знак, соответствующий элементу ak,k+1. Это позволяет сделать значение иk+1,k максимальным. Отметим, что методами Гивенса и Хаусхолдера можно пользоваться и в случае несимметричных матриц, приводя их, правда, не к трехдиагональному, а другому частному виду треугольной матрицы известной как матрица Гессенберга:
* | * | 0 | 0 | 0 | 0 |
* | * | * | 0 | 0 | 0 |
* | * | * | * | 0 | 0 |
* | * | * | * | * | 0 |
* | * | * | * | * | * |
* | * | * | * | * | * |
5. ОПРЕДЕЛЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИЙСИММЕТРИЧНОЙ ТРЕХДИАГОНАЛЬНОЙ МАТРИЦЫ
Приведя симметричную матрицу к трехдиагональному виду методом Гивенса или Хаусхолдера, необходимо найти ее собственные значения. Чтобы ясней были достоинства трехдиагональной формы, сформулируем задачу о собственных значениях в виде
dеt(А—lE) = 0,
где А — симметричная трехдиагональная матрица. Раcкрыв выражение в скобках, получим
a1 - l | b2 | 0 |
b1 | a2 - l | = 0 |
bn | ||
0 | bn | an - l |
Произвольный определитель порядка п можно выразить через п миноров порядка п — 1, каждый из которых в свою очередь выражается через п — 1 миноров порядка п — 2. Удобство трехдиагональной формы в том, что на каждом шаге все миноры, кроме двух, оказываются равными нулю. В результате исходный определитель представляется последовательностью полиномов
fm(l) = (am - l) fm-1 (l) – b2mfm-2(l).
Приняв
f0 (l) = 1 и f1 (l) = a1 - l при r = 2, .... п,
получим совокупность полиномов, известную как последовательность Штурма и обладающую тем свойством, что корни полинома fj (l) располагаются между корнями полинома fj+1(l). Поэтому для f1 (l) = a1—l можно утверждать, что значение lК = а1 заключено между корнями полинома f2 (l) == (a2 — l) (a1 — l) —b22. Это облегчает итерационное определение корней полинома, так как если известны границы интервалов, в которых лежат значения корней полинома, то их можно найти методом половинного деления. Так последовательно находят корни всех полиномов, и последний из них fn (l) дает все искомые п собственные значения. Эту процедуру можно проиллюстрировать графически (см. рис. 3).
Последовательность Штурма обладает еще и таким свойством: для любого значения b, при котором fn (b) <> 0, число собственных значений матрицы A, больших b, равно числу изменений знака последовательности
1, f1 (b), f2 (b), … , (1)n fn (b).
Если целое число, равное числу изменений знака, обозначить через V(b), то число собственных значений в интервале действительных чисел [b, с] будет равно V(b)—V(c).
Корни многочленаf2(l) f1 (b) |
Корни многочленаf3(l) f1 (b) |
Корни многочленаfn(l) f1 (b) |