Дискуссионный математический форумМатематический форум

Математический форум Math Help Planet

Обсуждение и решение задач по математике, физике, химии, экономике

Теоретический раздел
Часовой пояс: UTC + 4 часа [ Летнее время ]
MathHelpPlanet.com RSS-лента Математического форума

Часовой пояс: UTC + 4 часа [ Летнее время ]


Методы решения задач о собственных значениях и векторах матрицы

Методы решения задач о собственных
значениях и векторах матриц


Постановка задачи


Пусть [math]A[/math] — действительная числовая квадратная матрица размера [math](n\times n)[/math]. Ненулевой вектор [math]X= \bigl(x_1,\ldots,x_n\bigr)^T[/math] размера [math](n\times1)[/math], удовлетворяющий условию


[math]A\cdot X= \lambda\cdot X,\qquad \mathsf{(2.1)}[/math]

называется собственным вектором матрицы [math]A[/math]. Число [math]\lambda[/math] в равенстве (2.1) называется собственным значением. Говорят, что собственный вектор [math]X[/math] соответствует (принадлежит) собственному значению [math]\lambda[/math].


Равенство (2.1) равносильно однородной относительно [math]X[/math] системе:


[math](A-\lambda E)\cdot X=0\quad (X\ne 0).\qquad \mathsf{(2.2)}[/math]

Система (2.2) имеет ненулевое решение для вектора [math]X[/math] (при известном [math]\lambda[/math]) при условии [math]|A-\lambda E|=0[/math]. Это равенство есть характеристическое уравнение:


[math]|A-\lambda E|= P_n(\lambda)=0,[/math]
(2.3)

где [math]P_n(\lambda)[/math] — характеристический многочлен n-й степени. Корни [math]\lambda_1, \lambda_2,\ldots,\lambda_n[/math] характеристического уравнения (2.3) являются собственными (характеристическими) значениями матрицы [math]A[/math], а соответствующие каждому собственному значению [math]\lambda_i,~ i=1,\ldots,n[/math], ненулевые векторы [math]X^i[/math], удовлетворяющие системе


[math]AX^i=\lambda_iX^i\quad \text{or}\quad (A-\lambda_i E)X^i=0,~~ i=1,2,\ldots,n,[/math]
(2.4)

являются собственными векторами.


Требуется найти собственные значения и собственные векторы заданной матрицы. Поставленная задача часто именуется второй задачей линейной алгебры.


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


Различают полную и частичную проблему собственных значений, когда необходимо найти весь спектр (все собственные значения) и собственные векторы либо часть спектра, например: [math]\rho(A)= \max_{i}|\lambda_i(A)|[/math] и [math]\min_{i}|\lambda_i(A)|[/math]. Величина [math]\rho(A)[/math] называется спектральным радиусом.


Замечания


1. Если для собственного значения [math]\lambda_i[/math] — найден собственный вектор [math]X^i[/math], то вектор [math]\mu X^i[/math], где [math]\mu[/math] — произвольное число, также является собственным вектором, соответствующим этому же собственному значению [math]\lambda_i[/math].


2. Попарно различным собственным значениям соответствуют линейно независимые собственные векторы; k-кратному корню характеристического уравнения соответствует не более [math]k[/math] линейно независимых собственных векторов.


3. Симметрическая матрица имеет полный спектр [math]\lambda_i,~ i=\overline{1,n}[/math], действительных собственных значений; k-кратному корню характеристического уравнения симметрической матрицы соответствует ровно [math]k[/math] линейно независимых собственных векторов.


4. Положительно определенная симметрическая матрица имеет полный спектр действительных положительных собственных значений.




Метод непосредственного развертывания


Полную проблему собственных значений для матриц невысокого порядка [math](n\leqslant10)[/math] можно решить методом непосредственного развертывания. В этом случае будем иметь


[math]|A-\lambda E|= \begin{vmatrix}a_{11}-\lambda& a_{12}& a_{13}& \cdots& a_{1n}\\ a_{21}& a_{22}-\lambda& a_{23}& \cdots& a_{2n}\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ a_{n1}& a_{n2}& a_{n3}& \cdots& a_{nn}-\lambda \end{vmatrix}= P_n(\lambda)=0.[/math]
(2.5)

Уравнение [math]P_n(\lambda)=0[/math] является нелинейным (методы его решения изложены в следующем разделе). Его решение дает [math]n[/math], вообще говоря, комплексных собственных значений [math]\lambda_1,\lambda_2,\ldots,\lambda_n[/math], при которых [math]P_n(\lambda_i)=0~ (i=\overline{1,n})[/math]. Для каждого [math]\lambda_i[/math] может быть найдено решение однородной системы [math](A-\lambda_iE)X^i=0,~ i=\overline{1,n}[/math]. Эти решения [math]X^i[/math], определенные с точностью до произвольной константы, образуют систему [math]n[/math], вообще говоря, различных векторов n-мерного пространства. В некоторых задачах несколько этих векторов (или все) могут совпадать.


Алгоритм метода непосредственного развертывания


1. Для заданной матрицы [math]A[/math] составить характеристическое уравнение (2.5): [math]|A-\lambda E|=0[/math]. Для развертывания детерминанта [math]|A-\lambda E|[/math] можно использовать различные методы, например метод Крылова, метод Данилевского или другие методы.


2. Решить характеристическое уравнение и найти собственные значения [math]\lambda_1, \lambda_2, \ldots,\lambda_n[/math]. Для этого можно применить методы, изложенные далее.


3. Для каждого собственного значения составить систему (2.4):


[math](A-\lambda_iE)\cdot X^i=0,\quad i=1,2,\ldots,n[/math]

и найти собственные векторы [math]X^i[/math].


Замечание. Каждому собственному значению соответствует один или несколько векторов. Поскольку определитель [math]|A-\lambda_iE|[/math] системы равен нулю, то ранг матрицы системы меньше числа неизвестных: [math]\operatorname{rang}(A-\lambda_iE)=r<n[/math] и в системе имеется ровно [math]r[/math] независимых уравнений, а [math](n-r)[/math] уравнений являются зависимыми. Для нахождения решения системы следует выбрать [math]r[/math] уравнений с [math]r[/math] неизвестными так, чтобы определитель составленной системы был отличен от нуля. Остальные [math](n-r)[/math] неизвестных следует перенести в правую часть и считать параметрами. Придавая параметрам различные значения, можно получить различные решения системы. Для простоты, как правило, попеременно полагают значение одного параметра равным 1, а остальные равными 0.


Пример 2.1. Найти собственные значения и собственные векторы матрицы [math]A\in \mathbb{R}^{2\times 2}[/math], где [math]A=\begin{pmatrix}3&-2\\-4&1\end{pmatrix}[/math].


▼ Решение

Воспользуемся методикой.


1. Запишем уравнение (2.5): [math]|A-\lambda E|= \begin{vmatrix}3-\lambda&-2\\-4& 1-\lambda \end{vmatrix}= \lambda^2-4 \lambda-5=0[/math], отсюда получаем характеристическое уравнение [math]P_2(\lambda)\equiv \lambda^2-4 \lambda-5=0[/math].


2. Находим его корни (собственные значения): [math]\lambda_1=5,~ \lambda_2=-1[/math].


3. Составим систему [math](A-\lambda_iE)X^i=0,~ i=1,2[/math], для каждого собственного
значения и найдем собственные векторы:


[math]\begin{pmatrix}3-\lambda_1&-2\\-4& 1-\lambda_1 \end{pmatrix}\! \cdot\! \begin{pmatrix}x_1^1\\x_2^1\end{pmatrix}=0[/math] или [math]\begin{cases}-2x_1^1-2x_2^1=0,\\-4x_1^1-4x_2^1=0.\end{cases}[/math]

Отсюда [math]x_1^1=-x_2^1[/math]. Если [math]x_2^1=\mu[/math], то [math]x_1^1=-\mu[/math]. В результате получаем [math]X^1= \bigl\{x_1^1, x_2^1\bigr\}^T= \bigl\{\mu(-1;1)\bigr\}^T[/math].


Для [math]\lambda_2=-1[/math] имеем


[math]\begin{pmatrix}3-\lambda_2&-2\\-4& 1-\lambda_2 \end{pmatrix}\! \cdot\! \begin{pmatrix}x_1^2\\x_2^2\end{pmatrix}=0[/math] или [math]\begin{cases}4x_1^2-2x_2^2=0,\\-4x_1^2+2x_2^2=0.\end{cases}[/math]

Отсюда [math]x_2^2=2x_1^2[/math]. Если [math]x_1^2=\mu[/math], то [math]x_2^2=2\mu[/math]. В результате получаем [math]X^2= \bigl\{x_1^2, x_2^2\bigr\}^T= \bigl\{\mu(1;2)\bigr\}^T[/math], где [math]\mu[/math] — произвольное действительное число.


Пример 2.2. Найти собственные значения и собственные векторы матрицы [math]A= \begin{pmatrix}2&-1&1\\-1&2&-1\\0&0&1\end{pmatrix}[/math].


▼ Решение

Воспользуемся методикой.


1. Запишем характеристическое уравнение (2.5):


[math]\begin{vmatrix}2-\lambda&-1&1\\-1&2-\lambda&-1\\0&0&1-\lambda \end{vmatrix}=0[/math] или [math](1-\lambda)\bigl[(2-\lambda)^2-1\bigr]=0[/math].

2. Корни характеристического уравнения: [math]\lambda_{1,2}=1[/math] (кратный корень), [math]\lambda_3=3[/math] — собственные значения матрицы.


3. Найдем собственные векторы.


Для [math]\lambda_{1,2}=1[/math] запишем систему [math](A-\lambda_{1,2}E)\cdot X^{1,2}=0\colon[/math]


[math]\begin{pmatrix}1&-1&1\\-1&1&-1\\ 0&0&0\end{pmatrix}\! \cdot\! \begin{pmatrix} x_1^{1,2}\\ \x_2^{1,2}\\ x_3^{1,2}\end{pmatrix}=0.[/math]

Поскольку [math]\operatorname{rang}(A-\lambda_{1,2}E)=1[/math], в системе имеется одно независимое уравнение


[math]x_1^{1,2}-x_2^{1,2}+x_3^{1,2}=0[/math] или [math]x_1^{1,2}=x_2^{1,2}-x_3^{1,2}[/math].

Полагая [math]x_2^{1,2}=1,~ x_3^{1,2}=3[/math], получаем [math]x_1^{1,2}=1[/math] и собственный вектор [math]X^1= \begin{pmatrix}1&1&0\end{pmatrix}^T[/math].


Полагая [math]x_2^{1,2}=0,~ x_3^{1,2}=1[/math], получаем [math]x_1^{1,2}=-1[/math] и другой собственный вектор [math]X^2= \begin{pmatrix}-1&0&1\end{pmatrix}^T[/math]. Заметим, что оба собственных вектора линейно независимы.


Для собственного значения [math]\lambda_3=3[/math] запишем систему [math](A-\lambda_3E)\cdot X^3=0\colon[/math]


[math]\begin{pmatrix}-1&-1&1\\-1&-1&-1\\ 0&0&-2\end{pmatrix}\!\cdot\! \begin{pmatrix} x_1^3\\ x_2^3\\ x_3^3\end{pmatrix}=0.[/math]

Поскольку [math]\operatorname{rang}(A-\lambda_3E)=2[/math], то выбираем два уравнения:


[math]-x_1^3-x_2^3+x_3^3=0,\qquad-2x_3^3=0.[/math]

Отсюда [math]x_3^3=0,~ x_1^3=-x_2^3[/math]. Полагая [math]x_2^3=1[/math], получаем [math]x_1^3=-1[/math] и собственный вектор [math]X^3=\begin{pmatrix}-1&1&0 \end{pmatrix}^T[/math].




Метод итераций для нахождения собственных значений и векторов


Для решения частичной проблемы собственных значений и собственных векторов в практических расчетах часто используется метод итераций (степенной метод). На его основе можно определить приближенно собственные значения матрицы [math]A[/math] и спектральный радиус [math]\rho(A)= \max_{i}\bigl|\lambda_i(A)\bigr|[/math].


Пусть матрица [math]A[/math] имеет [math]n[/math] линейно независимых собственных векторов [math]X^i,~ i=1,\ldots,n[/math], и собственные значения матрицы [math]A[/math] таковы, что


[math]\rho(A)= \bigl|\lambda_1(A)\bigr|> \bigl|\lambda_2(A)\bigr|\geqslant \ldots\geqslant \bigl|\lambda_n(A)\bigr|.[/math]

Алгоритм метода итераций


1. Выбрать произвольное начальное (нулевое) приближение собственного вектора [math]X^{1(0)}[/math] (второй индекс в скобках здесь и ниже указывает номер приближения, а первый индекс без скобок соответствует номеру собственного значения). Положить [math]k=0[/math].


2. Найти [math]X^{1(1)}=AX^{1(0)},~ \lambda_1^{(1)}= \frac{x_i^{1(1)}}{x_i^{1(0)}}[/math], где [math]i[/math] — любой номер [math]1\leqslant i\leqslant n[/math], и положить [math]k=1[/math].


3. Вычислить [math]X^{1(k+1)}=A\cdot X^{1(k)}[/math].


4. Найти [math]\lambda_1^{(k+1)}= \frac{x_i^{1(k+1)}}{x_i^{1(k)}}[/math], где [math]x_i^{1(k+1)}, x_i^{1(k)}[/math] — соответствующие координаты векторов [math]X^{1(k+1)}[/math] и [math]X^{1(k)}[/math]. При этом может быть использована любая координата с номером [math]i,~ 1\leqslant i\leqslant n[/math].


5. Если [math]\Delta= \bigl|\lambda_1^{(k+1)}- \lambda_1^{(k)}\bigr|\leqslant \varepsilon[/math], процесс завершить и положить [math]\lambda_1\cong \lambda_1^{k+1}[/math]. Если [math]\Delta>\varepsilon[/math], положить [math]k=k+1[/math] и перейти к пункту 3.




Замечания


1. Процесс последовательных приближений


[math]\begin{aligned}&X^{1(1)}= AX^{1(0)},\quad X^{1(2)}= AX^{1(1)}= A^{2}X^{1(0)},\quad \ldots,\\ &X^{1(k)}= AX^{1(k-1)}= AA^{k-1}X^{1(0)}= A^kX^{1(0)},\quad \ldots \end{aligned}[/math]

сходится, т.е. при [math]x\to\infty[/math] вектор [math]X^{1(k)}[/math] стремится к собственному вектору [math]X^1[/math]. Действительно, разложим [math]X^{1(0)}[/math] по всем собственным векторам: [math]\textstyle{X^{1(0)}= \sum\limits_{i=1}^{n} c_iX^i}[/math]. Так как, согласно (2.4), [math]AX^i= \lambda_iX^i[/math], то


[math]\begin{aligned}& AX^{1(0)}= X^{1(1)}= \sum\limits_{i=1}^{n} c_i \lambda_iX^i,\quad AX^{1(1)}= A^2X^{1(0)}= X^{1(2)}= \sum\limits_{i=1}^{n} c_i \lambda_i^2X^i,\quad \ldots\\ &A^kX^{1(0)}= X^{1(k)}= \sum\limits_{i=1}^{n} c_i \lambda_i^kX^i= \lambda_1^k\! \left[c_1X^1+ c_2{\left(\frac{\lambda_2}{\lambda_1}\right)\!}^k X^2+ \ldots+ c_n{\left(\frac{\lambda_n}{\lambda_1}\right)\!}^k X^n\right]\!. \end{aligned}[/math]

При большом [math]k[/math] дроби [math]{\left(\frac{\lambda_2}{\lambda_1}\right)\!}^k, \ldots, {\left(\frac{\lambda_n}{\lambda_1}\right)\!}^k[/math] малы и поэтому [math]A^kX^{1(0)}= c_1\lambda_1^kX^1[/math], то есть [math]X^{1(k)}\to X^1[/math] при [math]k\to\infty[/math]. Одновременно [math]\lambda_1= \lim\limits_{k\to\infty} \frac{x_{i}^{1(k+1)}}{x_{i}^{1(k)}}[/math].


2. Вместо применяемой в пункте 4 алгоритма формулы для [math]\lambda_1^{(k+1)}[/math] можно взять среднее арифметическое соответствующих отношений для разных координат.


3. Метод может использоваться и в случае, если наибольшее по модулю собственное значение матрицы [math]A[/math] является кратным, т.е.


[math]\lambda_1= \lambda_2= \ldots= \lambda_s[/math] и [math]\bigl|\lambda_1\bigr|> \bigl|\lambda_k\bigr|[/math] при [math]k>s[/math].

4. При неудачном выборе начального приближения [math]X^{1(0)}[/math] предел отношения [math]\frac{x_i^{1(k+1)}}{x_i^{1(k)}}[/math] может не существовать. В этом случае следует задать другое начальное приближение.


5. Рассмотренный итерационный процесс для [math]\lambda_1[/math] сходится линейно, с параметром [math]c=\frac{\lambda_2}{\lambda_1}[/math] и может быть очень медленным. Для его ускорения используется алгоритм Эйткена.


6. Если [math]A=A^T[/math] (матрица [math]A[/math] симметрическая), то сходимость процесса при определении [math]\rho(A)[/math] может быть ускорена.


7. Используя [math]\lambda_1[/math], можно определить следующее значение [math]\lambda_2[/math] по формуле [math]\lambda_2= \frac{x_i^{1(k+1)}- \lambda_1 x_i^{1(k)}}{x_i^{1(k)}- \lambda_1 x_i^{1(k-1)}}~ (i=1,2,\ldots,n)[/math]. Эта формула дает грубые значения для [math]\lambda_2[/math], так как значение [math]\lambda_1[/math] является приближенным. Если модули всех собственных значений различны, то на основе последней формулы можно вычислять и остальные [math]\lambda_j~(j=3,4,\ldots,n)[/math].


8. После проведения нескольких итераций рекомендуется "гасить" растущие компоненты получающегося собственного вектора. Это осуществляется нормировкой вектора, например, по формуле [math]\frac{X^{1(k)}}{\|X^{1(k)}\|_1}[/math].


Пример 2.3. Для матрицы [math]A=\begin{pmatrix}5&1&2\\ 1&4&1\\ 2&1&3 \end{pmatrix}[/math] найти спектральный радиус степенным методом с точностью [math]\varepsilon=0,\,1[/math].


▼ Решение

1. Выбирается начальное приближение собственного вектора [math]X^{(0)}= \begin{pmatrix} 1&1&1 \end{pmatrix}^T[/math]. Положим [math]k=0[/math].


2. Найдем [math]X^{1(0)}= AX^{1(0)}= \begin{pmatrix}5&1&2\\ 1&4&1\\ 2&1&3\end{pmatrix}\!\cdot\! \begin{pmatrix}1\\1\\1\end{pmatrix}= \begin{pmatrix}8\\6\\6\end{pmatrix}[/math]; [math]\lambda_1^{(1)}= \frac{x_1^{1(1)}}{x_1^{1(0)}}= \frac{8}{1}=8[/math], положим [math]k=1[/math].


3. Вычислим [math]X^{1(2)}= AX^{1(1)}= \begin{pmatrix}5&1&2\\ 1&4&1\\ 2&1&3 \end{pmatrix}\!\cdot\! \begin{pmatrix}8\\6\\6\end{pmatrix}= \begin{pmatrix} 58\\38\\40 \end{pmatrix}[/math].


4. Найдем [math]\lambda_1^{(2)}= \frac{x_1^{1(2)}}{x_1^{1(1)}}= \frac{58}{8}=7,\!25[/math].


5. Так как [math]\bigl|\lambda_1^{(2)}- \lambda_1^{(1)}\bigr|= 0,\!75> \varepsilon[/math], то процесс необходимо продолжить. Результаты вычислений удобно представить в виде табл. 10.10.


[math]\begin{array}{|c|c|c|c|c|c|}\hline \phantom{\begin{matrix}|\\.\end{matrix}} k \phantom{\begin{matrix}|\\.\end{matrix}} & x_{1}^{1(k)} & x_{2}^{1(k)} & x_{3}^{1(k)} & \lambda_1^{(k)} & \bigl|\lambda_1^{(k)}- \lambda_1^{(k-1)}\bigr| \\ \hline 0& 1 & 1 & 1 &- &- \\ \hline 1& 8& 6& 6 & 8 &- \\ \hline 2& 58& 38& 40& 7,\!25& 0,\!75\\ \hline 3& 408& 250& 274& 7,\!034& 0,\!116\\ \hline 4& 2838& 1682& 1888& 6,\!9559& 0,\!078< \varepsilon\\ \hline \end{array}[/math]

Точность по достигнута на четвертой итерации. Таким образом, в качестве приближенного значения [math]\lambda_1[/math] берется 6,9559, а в качестве собственного вектора принимается [math]X^1= \begin{pmatrix} 2838& 1682& 1888\end{pmatrix}^T[/math].


Так как собственный вектор определяется с точностью до постоянного множителя, то [math]X^1[/math] лучше пронормировать, т.е. поделить все его компоненты на величину нормы. Для рассматриваемого примера получим


[math]X^1= \frac{1}{2838}\cdot\! \begin{pmatrix}2838\\ 1682\\ 1888\end{pmatrix}= \begin{pmatrix}1,\!000\\ 0,\!5927\\ 0,\!6652 \end{pmatrix}\!.[/math]

Согласно замечаниям, в качестве собственного значения [math]\lambda_1[/math] матрицы можно взять не только отношение


[math]\frac{x_1^{1(4)}}{x_1^{1(3)}}= \frac{2838}{408}\approx 6,\!9559[/math], но и [math]\frac{x_2^{1(4)}}{x_2^{1(3)}}= \frac{1682}{250}\approx 6,\!7280;~~ \frac{x_3^{1(4)}}{x_3^{1(3)}}= \frac{1888}{274}\approx 6,\!8905[/math],

а также их среднее арифметическое [math]\frac{6,\!9559+6,\!728+6,\!8905}{3}\approx 6,\!8581[/math].


Пример 2.4. Найти максимальное по модулю собственное значение матрицы [math]A=\begin{pmatrix}2&-1&1\\ -1&2&-1\\ 0&0&3 \end{pmatrix}[/math] и соответствующий собственный вектор.


▼ Решение

1. Зададим начальное приближение [math]X^{1(0)}= \begin{pmatrix}1&-1&1 \end{pmatrix}^T[/math] и [math]\varepsilon=0,\!0001[/math].


Выполним расчеты согласно методике (табл. 10.11).


[math]\begin{array}{|c|c|c|c|c|c|}\hline \phantom{\begin{matrix}|\\.\end{matrix}} k \phantom{\begin{matrix}|\\.\end{matrix}} & x_{1}^{1(k)} & x_{2}^{1(k)} & x_{3}^{1(k)} & \lambda_1^{(k)} & \bigl|\lambda_1^{(k)}-\lambda_1^{(k-1)}\bigr| \\ \hline 0 & 1&-1& 1&-&-\\ \hline 1 & 4&-4& 1& 4&-\\ \hline 2 & 13&-13& 1& 3,\!25& 0,\!75 \\ \hline 3 & 40&-40& 1& 3,\!0769& 0,\!17307\\ \hline 4 & 121& -121& 1& 3,\!025& 0,\!0519\\ \hline 5 & 364&-364& 1& 3,\!00826& 0,\!01673\\ \hline 6 & 1093&-1093& 1& 3,\!002747& 0,\!005512\\ \hline 7 & 3280&-3280 & 1& 3,\!000914& 0,\!00183\\ \hline 8 & 9841& -9841& 1& 3,\!000304& 0,\!000609\\ \hline 9 & 29524& -29524& 1& 3,\!000101& 0,\!000202\\ \hline 10 & 88573& -88573& 1& 3,\!000034& 0,\!000067\\ \hline \end{array}[/math]

В результате получено собственное значение [math]\lambda_1\cong 3,\!00003[/math] и собственный вектор [math]X^1= \begin{pmatrix} 88573&-88573&1\end{pmatrix}^T[/math] или после нормировки


[math]X^1= \frac{1}{88573}\cdot \begin{pmatrix} 88573&-88573&1 \end{pmatrix}^T = \begin{pmatrix} 1&-1&0,\!0000113\end{pmatrix}^T.[/math]



Метод вращений для нахождения собственных значений


Метод используется для решения полной проблемы собственных значений симметрической матрицы и основан на преобразовании подобия исходной матрицы [math]A\in\mathbb{R}^{n\times n}[/math] с помощью ортогональной матрицы [math]H[/math].


Напомним, что две матрицы [math]A[/math] и [math]A^{(i)}[/math] называются подобными ([math]A\sim A^{(i)}[/math] или [math]A^{(i)}\sim A[/math]), если [math]A^{(i)}=H^{-1}AH[/math] или [math]A=HA^{(i)}H^{-1}[/math], где [math]H[/math] — невырожденная матрица.


В методе вращений в качестве [math]H[/math] берется ортогональная матрица, такая, что [math]HH^{T}=H^{T}H=E[/math], т.е. [math]H^{T}=H^{-1}[/math]. В силу свойства ортогонального преобразования евклидова норма исходной матрицы [math]A[/math] не меняется. Для преобразованной матрицы [math]A^{(i)}[/math] сохраняется ее след и собственные значения [math]\lambda_i\colon[/math]


[math]\operatorname{tr}A= \sum_{i=1}^{n}a_{ii}= \sum_{i=1}^{n} \lambda_i(A)= \operatorname{tr}A^{(i)}.[/math]

При реализации метода вращений преобразование подобия применяется к исходной матрице [math]A[/math] многократно:


[math]A^{(k+1)}= \bigl(H^{(k)}\bigr)^{-1}\cdot A^{(k)}\cdot H^{(k)}= \bigl(H^{(k)}\bigr)^{T}\cdot A^{(k)}\cdot H^{(k)},\quad k=0,1,2,\ldots[/math]
(2.6)

Формула (2.6) определяет итерационный процесс, где начальное приближение [math]A^{(0)}=A[/math]. На k-й итерации для некоторого выбираемого при решении задачи недиагонального элемента [math]a_{ij}^{(k)},~ i\ne j[/math], определяется ортогональная матрица [math]H^{(k)}[/math], приводящая этот элемент [math]a_{ij}^{(k+1)}[/math] (а также и [math]a_{ji}^{(k+1)}[/math]) к нулю. При этом на каждой итерации в качестве [math]a_{ij}^{(k+1)}[/math] выбирается наибольший по модулю. Матрица [math]H^{(k)}[/math] называемая матрицей вращения Якоби, зависит от угла [math]\varphi^{(k)}[/math] и имеет вид


Матрица вращения Якоби

В данной ортогональной матрице элементы на главной диагонали единичные, кроме [math]h_{ii}^{(k)}= \cos\varphi^{(k)}[/math] и [math]h_{jj}^{(k)}=\cos\varphi^{(k)}[/math], а остальные элементы нулевые, за исключением [math]h_{ij}^{(k)}=-\sin\varphi^{(k)}[/math], [math]h_{ji}^{(k)}=\sin\varphi^{(k)}[/math] ([math]h_{ij}[/math] -элементы матрицы [math]H[/math]).


Угол поворота [math]\varphi^{(k)}[/math] определяется по формуле


[math]\operatorname{tg}2\varphi^{(k)}= \frac{2a_{ij}^{(k)}}{a_{ii}^{(k)}-a_{jj}^{(k)}}= \overline{P}_k\,;\qquad \varphi^{(k)}= \frac{1}{2}\operatorname{arctg}\overline{P}_k\,,[/math]
(2.7)

где [math]|2\varphi^{(k)}|\leqslant \frac{\pi}{2},~ i<j[/math] ([math]a_{ij}[/math] выбирается в верхней треугольной наддиагональной части матрицы [math]A[/math]).


В процессе итераций сумма квадратов всех недиагональных элементов [math]\sigms(A^{(k)})[/math] при возрастании [math]k[/math] уменьшается, так что [math]\sigms(A^{(k+1)})< \sigma(A^{(k)})[/math]. Однако элементы [math]a_{ij}^{(k)}[/math] приведенные к нулю на k-й итерации, на последующей итерации немного возрастают. При [math]k\to\infty[/math] получается монотонно убывающая ограниченная снизу нулем последовательность [math]\sigma(A^{(1)})> \sigma(A^{(2)})> \ldots> \sigma(A^{(k)})>\ldots[/math]. Поэтому [math]\sigma(A^{(k)})\to0[/math] при [math]k\to\infty[/math]. Это и означает сходимость метода. При этом [math]A^{(k)}\to \Lambda= \operatorname{diag}(\lambda_1,\ldots,\lambda_n)[/math].


Замечание. В двумерном пространстве с введенной в нем системой координат [math]Oxy[/math] с ортонормированным базисом [math]\{\vec{i},\vec{j}\}[/math] матрица вращения легко получается из рис. 2.1, где система координат [math]Ox'y'[/math] повернута на угол [math]\varphi\colon[/math]


[math]\begin{cases}\vec{i}\,'= \vec{i}\cos \varphi+ \vec{j}\sin \varphi,\\[2pt] \vec{j}\,'=-\vec{i}\sin \varphi+ \vec{j}\cos \varphi. \end{cases}[/math]

Таким образом, для компонент [math]\vec{i}\,',\, \vec{j}\,'[/math] будем иметь [math]\bigl(\vec{i}\,',\vec{j}\,'\bigr)= \bigl(\vec{i},\vec{j}\bigr)\cdot\! \begin{pmatrix} \cos \varphi&-\sin \varphi\\ \sin \varphi& \cos \varphi\end{pmatrix}[/math]. Отсюда следует, что в двумерном пространстве матрица вращения имеет вид [math]H= \begin{pmatrix} \cos \varphi&-\sin \varphi\\ \sin \varphi& \cos \varphi\end{pmatrix}[/math]. Отметим, что при [math]n=2[/math] для решения задачи требуется одна итерация.


Графическая интерпретация матрицы вращения Якоби



Алгоритм метода вращений


1. Положить [math]k=0,~ A^{(0)}=A[/math] и задать [math]\varepsilon>0[/math].


2. Выделить в верхней треугольной наддиагональной части матрицы [math]A^{(k)}[/math] максимальный по модулю элемент [math]a_{ij}^{(k)},~ i<j[/math].


Если [math]|a_{ij}^{(k)}|\leqslant \varepsilon[/math] для всех [math]i\ne j[/math], процесс завершить. Собственные значения определяются по формуле [math]\lambda_i(A^{(k)})=a_{ii}^{(k)},~ i=1,\ldots,n[/math].


Собственные векторы [math]X^i[/math] находятся как i-e столбцы матрицы, получающейся в результате перемножения:


[math]\nu_k= H^{(0)}\cdot H^{(1)}\cdot H^{(2)}\cdot \ldots\cdot H^{(k-1)}= \bigl( X^1,X^2,X^3,\ldots,X^n\bigr).[/math]

Если [math]\bigl|a_{ij}^{(k)}\bigr|>\varepsilon[/math], процесс продолжается.


3. Найти угол поворота по формуле [math]\varphi^{(k)}= \frac{1}{2} \operatorname{arctg} \frac{2a_{ij}^{(k)}}{a_{ii}^{(k)}-a_{jj}^{(k)}}[/math].


4. Составить матрицу вращения [math]H^{(k)}[/math].


5. Вычислить очередное приближение [math]A^{(k+1)}= \bigl(H^{(k)}\bigr)^T A^{(k)} H^{(k)}[/math].Положить [math]k=k+1[/math] и перейти к пункту 2.


Замечания


1. Используя обозначение [math]\overline{P}_k= \frac{2a_{ij}^{(k)}}{a_{ii}^{(k)}-a_{jj}^{(k)}}[/math], можно в пункте 3 алгоритма вычислять элементы матрицы вращения по формулам


[math]\sin\varphi^{(k)}= \operatorname{sign}\overline{P}_k\cdot \sqrt{\frac{1}{2}\!\left(1-\frac{1}{\sqrt{1+\overline{P}_{k}^{\;2}}}\right)}\,,\qquad \cos\varphi^{(k)}= \sqrt{\frac{1}{2}\!\left(1+ \frac{1}{\sqrt{1+\overline{P}_{k}^{\;2}}}\right)}\,.[/math]

2. Контроль правильности выполнения действий по каждому повороту осуществляется путем проверки сохранения следа преобразуемой матрицы.


3. При [math]n=2[/math] для решения задачи требуется одна итерация.


Пример 2.5. Для матрицы [math]A=\begin{pmatrix} 2&1\\1&3 \end{pmatrix}[/math] методом вращений найти собственные значения и собственные векторы.


▼ Решение

1. Положим [math]k=0,~A^{(0)}=A= \begin{pmatrix} 2&1\\1&3 \end{pmatrix}\!,~ \varepsilon=10^{-10}[/math].


2°. Выше главной диагонали имеется только один элемент [math]a_{ij}=a_{12}=1[/math].


3°. Находим угол поворота матрицы по формуле (2.7), используя в расчетах 11 цифр после запятой в соответствии с заданной точностью:


[math]\operatorname{tg}2\varphi^{(0)}= \frac{2a_{ij}}{a_{ii}-a_{jj}}= \frac{2}{2-3}=-2\quad \Rightarrow\quad\! \left[\begin{aligned} \sin\varphi^{(0)}&\approx -0,\!52573111212\,;\\ \cos\varphi^{(0)}&\approx 0,\!85065080835\,.\end{aligned}\right.[/math]

4°. Сформируем матрицу вращения:


[math]H^{(0)}= \begin{pmatrix}\cos\varphi^{(0)}& -\sin\varphi^{(0)}\\[2pt] \sin\varphi^{(0)}& \cos\varphi^{(0)} \end{pmatrix}= \begin{pmatrix} 0,\!85065080835 & 0,\!52573111212 \\[2pt] -0,\!52573111212 & 0,\!85065080835 \end{pmatrix}\!.[/math]

5°. Выполним первую итерацию:


[math]\begin{aligned}A^{(1)}= \bigl(H^{(0)}\bigr)^T\cdot A^{(0)}\cdot H^{(0)}&= \begin{pmatrix}0,\!85065080835 & -0,\!52573111212 \\[2pt] 0,\!52573111212 & 0,\!85065080835 \end{pmatrix}\!\cdot A^{(0)}\cdot \begin{pmatrix}0,\!85065080835 & 0,\!52573111212 \\[2pt] -0,\!52573111212 & 0,\!85065080835\end{pmatrix}=\\ &= \begin{pmatrix} 1,\!38196601125& -4,\!04620781325\cdot10^{-12}\\ -4,\!04587474634\cdot10^{-12}& 3,\!61803398874 \end{pmatrix}\!. \end{aligned}[/math]

Очевидно, след матрицы с заданной точностью сохраняется, т.е. [math]\sum_{i=1}^{2}a_{ii}^{(0)}= \sum_{i=1}^{2}a_{ii}^{(0)}=5[/math]. Положим [math]k=1[/math] и перейдем к пункту 2.


2. Максимальный по модулю наддиагональный элемент [math]|a_{12}|= 4,\!04620781325\cdot10^{-12} < \varepsilon=10^{-10}[/math]. Для решения задачи (подчеркнем, что [math]n=2[/math]) с принятой точностью потребовалась одна итерация, полученную матрицу можно считать диагональной. Найдены следующие собственные значения и собственные векторы:


[math]\lambda_1= 1,\!38196601125\,.\quad \lambda_2=3,\!61803398874\,; \qquad X^1= \begin{pmatrix} 0,\!85065080835\\ -0,\!52573111212\end{pmatrix}\!;\quad X^2=\begin{pmatrix} 0,\!52573111212\\ 0,\!85065080835 \end{pmatrix}\!.[/math]

Пример 2.6. Найти собственные значения и собственные векторы матрицы [math]A=\begin{pmatrix}5&1&2\\ 1&4&1\\ 2&1&3 \end{pmatrix}[/math].


▼ Решение

1. Положим [math]k=0,~ A^{(0)}=A,~ \varepsilon= 0,\!001[/math].


2°. Выделим максимальный по модулю элемент в наддиагональнои части: [math]a_{13}^{(0)}=2[/math]. Так как [math]a_{13}=2> \varepsilon=0,\!001[/math], то процесс продолжается.


3°. Находим угол поворота:


[math]\varphi^{(0)}= \frac{1}{2}\operatorname{arctg} \frac{2a_{13}^{(0)}}{a_{11}^{(0)}-a_{33}^{(0)}}= \frac{1}{2}\operatorname{arctg}\frac{4}{5-3}= \frac{1}{2} \operatorname{arctg}2\approx 0,\!553574\quad \Rightarrow\quad\! \left[\begin{aligned} \sin\varphi^{(0)}&\approx 0,\!52573\,;\\ \cos\varphi^{(0)}&\approx 0,\!85065\,. \end{aligned}\right.[/math]

4°. Сформируем матрицу вращения: [math]H^{(0)}= \begin{pmatrix}0,\!85065&0&-0,\!52573\\ 0&1&0\\ 0,\!52573&0&0,\!85065 \end{pmatrix}[/math].


5°. Выполним первую итерацию: [math]A^{(1)}= \bigl(H^{(0)}\bigr)^T A^{(0)}H^{(0)}= \begin{pmatrix} 6,\!236&1,\!376&2,\!33\cdot10^{-6}\\ 1,\!376&4&0,\!325\\ 2,\!33\cdot10^{-6}&0,\!325&1,\!764 \end{pmatrix}[/math]. Положим [math]k=1[/math] и перейдем к пункту 2.


2(1). Максимальный по модулю наддиагональный элемент [math]a_{12}^{(1)}=1,\!376[/math]. Так как [math]a_{12}^{(1)}> \varepsilon=0,\!001[/math], процесс продолжается.


3(1). Найдем угол поворота:


[math]\varphi^{(1)}= \frac{1}{2} \operatorname{arctg} \frac{2a_{12}^{(1)}}{a_{11}^{(1)}-a_{22}^{(1)}}= \frac{1}{2} \operatorname{arctg}\frac{2\cdot1,\!376}{6,\!236-4}= \frac{1}{2} \operatorname{arctg}1,\!230769\approx 0,\!444239\quad \Rightarrow\quad\! \left[\begin{aligned} \sin\varphi^{(1)}\approx 0,\!429770\,;\\ \cos\varphi^{(1)}\approx 0,\!902937\,. \end{aligned}\right.[/math]

4(1). Сформируем матрицу вращения: [math]H^{(1)}= \begin{pmatrix} 0,\!902937&-0,\!429770&0\\ 0,\!429770&0,\!902937&0\\ 0&0&1 \end{pmatrix}[/math].


5(1). Выполним вторую итерацию: [math]A^{(2)}= \bigl(H^{(1)}\bigr)^T A^{(1)}H^{(1)}= \begin{pmatrix} 6,\!891& 2,\!238\cdot10^{-4}&0,\!14\\ 2,\!238\cdot10^{-4}& 3,\!345&0,\!293\\ 0,\!14&0,\!293&1,\!764 \end{pmatrix}[/math]. Положим [math]k=2[/math] и перейдем к пункту 2.


2(2). Максимальный по модулю наддиагональный элемент [math]a_{23}^{(2)}=0,\!293> \varepsilon=0,\!001[/math].


3(2). Найдем угол поворота:


[math]\varphi^{(2)}= \frac{1}{2} \operatorname{arctg} \frac{2a_{23}^{(2)}}{a_{22}^{(2)}-a_{33}^{(2)}}= \frac{1}{2} \operatorname{arctg}\frac{2\cdot0,\!293}{3,\!345-1,\!764}= \frac{1}{2} \operatorname{arctg}0,\!370651\approx 0,\!177476\quad \Rightarrow\quad\! \left[\begin{aligned} \sin\varphi^{(2)}\approx 0,\!1765460\,;\\ \cos\varphi^{(2)}\approx 0,\!9842924\,. \end{aligned}\right.[/math]

4(2). Сформируем матрицу вращения [math]H^{(2)}= \begin{pmatrix}1&0&0\\ 0&0,\!9842924& -0,\!1765460\\ 0& 0,\!1765460& 0,\!9842924\end{pmatrix}[/math].


5(2). Выполним третью итерацию и положим [math]k=3[/math] и перейдем к пункту 2:


[math]A^{(3)}=\bigl(H^{(2)}\bigr)^T\cdot A^{(2)}\cdot H^{(2)}= \ldots= \begin{pmatrix} 6,\!891&0,\!025&0,\!138\\ 0,\!025&3,\!398&3,\!375\cdot 10^{-7}\\ 0,\!138&3,\!375\cdot10^{-7}& 1,\!711 \end{pmatrix}\!.[/math]

2(3). Максимальный по модулю наддиагональный элемент [math]a_{13}^{(3)}= 0,\!138>\varepsilon[/math].


3(3). Найдем угол поворота:


[math]\varphi^{(3)}= \frac{1}{2}\operatorname{arctg} \frac{2\cdot a_{13}^{(3)}}{a_{11}^{(3)}-a_{33}^{(3)}}= \frac{1}{2}\operatorname{arctg} \frac{2\cdot0,\!138}{6,\!891-1,\!711}= \frac{1}{2}\operatorname{arctg}0,\!05328\approx 0,\!026615 \quad\Rightarrow\quad\! \left[\begin{aligned} \sin\varphi^{(3)}\approx 0,\!026611\,;\\ \cos\varphi^{(3)}\approx 0,\!999646\,. \end{aligned}\right.[/math]

4(3). Сформируем матрицу вращения: [math]H^{(3)}= \begin{pmatrix} 0,\!999646&0&-0,\!026611\\ 0&1&0\\ 0,\!026611&0&0,\!999646 \end{pmatrix}[/math].


5(3). Выполним четвертую итерацию и положим [math]k=4[/math] и перейдем к пункту 2:


[math]A^{(4)}=\bigl(H^{(3)}\bigr)^T\cdot A^{(3)}\cdot H^{(3)}= \ldots= \begin{pmatrix} 6,\!895&0,\!025&8,\!406\cdot10^{-6}\\ 0,\!025&3,\!398&-6,\!649\cdot10^{-6}\\ 8,\!406\cdot10^{-6}& -6,\!649\cdot10^{-6}&1,\!707 \end{pmatrix}\!.[/math]

2(4). Так как [math]a_{12}^{(4)}=0,\!025>\varepsilon[/math], процесс повторяется


3(4). Найдем угол поворота


[math]\varphi^{(4)}= \frac{1}{2}\operatorname{arctg} \frac{2\cdot a_{12}^{(4)}}{a_{11}^{(4)}-a_{22}^{(4)}}= \frac{1}{2}\operatorname{arctg} \frac{2\cdot0,\!025}{6,\!895-3,\!398}\approx 0,\!0071484 \quad\Rightarrow\quad\! \left[\begin{aligned} \sin\varphi^{(4)}\approx 0,\!0071483\,;\\ \cos\varphi^{(4)}\approx 0,\!9999744\,. \end{aligned}\right.[/math]

4(4). Сформируем матрицу вращения: [math]H^{(4)}= \begin{pmatrix} 0,\!9999744&-0,\!0071483&0\\ 0,\!0071483&0,\!9999744&0\\ 0&0&1 \end{pmatrix}[/math].


5(4). Выполним пятую итерацию и положим [math]k=5[/math] и перейдем к пункту 2:


[math]A^{(5)}=\bigl(H^{(4)}\bigr)^T\cdot A^{(4)}\cdot H^{(4)}= \ldots= \begin{pmatrix} 6,\!895&4,\!774\cdot10^{-7}&3,\!653\cdot10^{-6}\\ 4,\!774\cdot10^{-7}&3,\!398&-6,\!649\cdot10^{-4}\\ 3,\!653\cdot10^{-6}& -6,\!649\cdot10^{-4}&1,\!707 \end{pmatrix}\!.[/math]

2(5). Так как наибольший по модулю наддиагональный элемент удовлетворяет условию [math]\bigl|-6,\!649\cdot10^{-4}\bigr|<\varepsilon=0,\!001[/math], процесс завершается.


Собственные значения: [math]\lambda_1\cong a_{11}^{(5)}= 6,\!895\,,~ \lambda_2\cong a_{22}^{(5)}=3,\!398\,,~ \lambda_3\cong a_{33}^{(5)}=1,\!707\,,[/math]. Для нахождения собственных векторов вычислим


[math]\nu_5=H^{(0)}\cdot H^{(1)}\cdot H^{(2)}\cdot H^{(3)}\cdot H^{(4)}= \begin{pmatrix} 0,\!753&-0,\!458&-0,\!473\\ 0,\!432&0,\!886&-0,\!171\\ 0,\!497&-0,\!076&0,\!864 \end{pmatrix}\!,[/math]

отсюда [math]X^1=\begin{pmatrix} 0,\!753\\0,\!432\\0,\!497 \end{pmatrix}\!,~ X^2=\begin{pmatrix}-0,\!458\\ 0,\!886\\ -0,\!076 \end{pmatrix}\!,~ X^3=\begin{pmatrix}-0,\!473\\-0,\!171\\0,\!864 \end{pmatrix}[/math] или после нормировки


[math]X^1=\begin{pmatrix}1\\0,\!5737\\0,\!660\end{pmatrix}\!,\qquad X^2=\begin{pmatrix}-0,\!517\\ 1\\ -0,\!0858\end{pmatrix}\!,\qquad X^3=\begin{pmatrix}-0,\!5474\\-0,\!1979\\1 \end{pmatrix}\!.[/math]

Часовой пояс: UTC + 4 часа [ Летнее время ]


Яндекс.Метрика

Copyright © 2010-2016 MathHelpPlanet.com. All rights reserved