Участник:Алексей Будюк/Метод Якоби вычисления собственных значений симметричной матрицы

Основные авторы описания: А.М. Будюк, В.А. Сальников

Свойства и структура алгоритма

Общее описание алгоритма

Целый ряд инженерных задач сводится к рассмотрению систем уравнений, имеющих единственное решение лишь в том случае, если известно значение некоторого входящего в них параметра. Этот особый параметр называется характеристическим, или собственным значением системы. С задачами на собственные значения инженер сталкивается в различных ситуациях. Так, для тензоров напряжений собственные значения определяют главные нормальные напряжения, а собственными векторами задаются направления, связанные с этими значениями. При динамическом анализе механических систем собственные значения соответствуют собственным частотам колебаний, а собственные вектора характеризуют моды этих колебаний. При расчете конструкций собственные значения позволяют определять критические нагрузки, превышение которых приводит к потере устойчивости.
Таким образом, задача нахождения собственных значений и собственных векторов матриц(СЗВМ) является одной из основных задач для многих разделов физики. Одним из методов решения является метод Якоби. Этот метод позволяет привести матрицу к диагональному виду, последовательно, исключая все элементы, стоящие вне главной диагонали. В итоге, накопление в процессе преобразований произведения трансформационных матриц дает матрицу собственных векторов, в то время как диагональные элементы являются собственными значениями. К сожалению, приведение к строго диагональному виду, в общем случае, требует бесконечно большого числа шагов, так как образование нового нулевого элемента на месте одного из элементов матрицы часто ведет к появлению ненулевого элемента там, где ранее был нуль, так как при преобразовании затрагиваются и другие элементы матрицы. На практике метод Якоби рассматривают, как итерационную процедуру, которая в принципе позволяет с заданной точностью подойти к диагональной форме.
В случае симметричной матрицы [math]\displaystyle{ A }[/math] действительных чисел преобразование выполняется с помощью ортогональных матриц [math]\displaystyle{ P }[/math], полученных в результате вращении в действительной плоскости. Вычисления осуществляются следующим образом. Из исходной матрицы [math]\displaystyle{ A }[/math] образуют матрицу [math]\displaystyle{ A_1 }[/math], при этом ортогональная матрица плоского вращения [math]\displaystyle{ P_1 }[/math] выбирается так, чтобы в матрице [math]\displaystyle{ A_1 }[/math] появился нулевой элемент, стоящий вне главной диагонали (обычно выбирается наибольший по модулю элемент). Затем из [math]\displaystyle{ A_1 }[/math] с помощью второй преобразующей матрицы [math]\displaystyle{ P_2, }[/math] образуют новую матрицу [math]\displaystyle{ A_2 }[/math]. При этом [math]\displaystyle{ P_2 }[/math], выбирают так, чтобы в [math]\displaystyle{ A_2 }[/math] появился еще один нулевой внедиагональный элемент. Эту процедуру продолжают, стремясь, чтобы на каждом шаге в нуль обращался наибольший внедиагональный элемент. В итоге получаем некоторый бесконечный итерационный процесс, на каждом шаге которого преобразующаяся матрица становится всё более и более близкой к диагональной форме (в смысле уменьшения суммы квадратов всех своих недиагональных элементов). На практике, при решении задач, выбирается некоторая величина [math]\displaystyle{ \epsilon }[/math], и итерационный процесс останавливается, когда наибольший по модулю внедиагональный элемент становится меньше [math]\displaystyle{ \epsilon }[/math].

Математическое описание алгоритма

Пусть [math]\displaystyle{ A=(a_{ij}) }[/math] симметричная матрица, а [math]\displaystyle{ P(i,j,\theta) }[/math] ортогональная матрица, называемая вращением Якоби, зависящая от параметра [math]\displaystyle{ \theta }[/math], тогда итерационная последовательность строится следующим образом:

[math]\displaystyle{ A_k = P_k^TA_{k-1}P_k }[/math]

Матрицы [math]\displaystyle{ A_i }[/math] образуют итерационную поледовательность симетричных матриц, сходящуюся к диагональной матрице из собственных значений [math]\displaystyle{ D }[/math]. А так как

[math]\displaystyle{ A_k = P_k^TA_{k-1}P_k = P_k^TP_{k-1}^TA_{k-2}P_kP_{k-1} = ... = P^TA_1P }[/math],

то матрицы вращений образуют матрицу из собственных векторов [math]\displaystyle{ P=P_kP_{k-1}..P_1 }[/math]. Матрица вращения [math]\displaystyle{ P(i,j,\theta) }[/math] имеет следующий вид:

[math]\displaystyle{ P(i,j,\theta)=\begin{pmatrix} 1 & & & & & & 0 \\ &\ddots & & & & & \\ & & cos(\theta) & & sin(\theta) & & \\ & & & \ddots & & & \\ & & -sin(\theta) & & cos(\theta) & & \\ & & & & &\ddots & \\ 0 & & & & & & 1 \end{pmatrix} }[/math]

Где синусы и косинусы стоят только на [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math] позициях.
Если же мы произведём матричное умножение [math]\displaystyle{ P_k^TA_kP_k }[/math], то в результате первого умножения [math]\displaystyle{ P_k^TA_k }[/math] изменит только строки [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math] матрицы [math]\displaystyle{ A_k }[/math], в то время как [math]\displaystyle{ A_kP_k }[/math] меняет только столбцы [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math]. Таким образом изменённые элементы расположены только в строках и столбцах [math]\displaystyle{ i }[/math] и [math]\displaystyle{ j }[/math] матрицы [math]\displaystyle{ A_k }[/math]. И в результате этого перемножения мы получаем формулы для элементов изменённой матрицы:

[math]\displaystyle{ \begin{align} a_{ii}^{k+1}&=c^2a_{ii}^k-2sca_{ij}^k+s^2a_{jj}^k \\ a_{jj}^{k+1}&=s^2a_{ii}^k+2sca_{ij}^k+c^2a_{jj}^k \\ a_{ij}^{k+1}&=a_{ji}^{k+1}=(c^2-s^2)a_{ij}^k+sc(a_{ii}^k-a_{jj}^k)\\ a_{il}^{k+1}&=a_{li}^{k+1}=ca_{il}^k-sa_{jl}^k\\ a_{jl}^{k+1}&=a_{lj}^{k+1}=sa_{il}^k+ca_{jl}^k\\ a_{lm}^{k+1}&=a_{ml}^k,\\ \end{align} }[/math]

[math]\displaystyle{ l,m\neq i,j }[/math]

Приравнивая внедиагональный элемент к нулю, получаем формулы для параметра [math]\displaystyle{ \theta }[/math]: [math]\displaystyle{ r=ctg(2\theta)=\frac{c^2-s^2}{2sc}=\frac{a_{jj}-a_{ii}}{2a_{ij}} }[/math]

[math]\displaystyle{ tg(\theta)=t=\frac{sign(r)}{|r|+\sqrt{1+r^2}} }[/math]

[math]\displaystyle{ c=\frac{1}{\sqrt{1+t^2}} }[/math]

[math]\displaystyle{ s=ct }[/math]

Макроструктура алгоритма

кек