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

Форум пользователей пакета Maple

Модератор: Admin

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

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

Сообщение Кли » Чт окт 04, 2018 9:51 pm

Правильно ли я понимаю, для этого метода, такие формул нужны?
https://ru.wikipedia.org/wiki/Метод_Яко ... х_значений
и если да, то эти формулы из в цикл надо вбивать?

Kitonum
Сообщения: 2050
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Re: Метод вращения Якоби

Сообщение Kitonum » Пт окт 05, 2018 11:06 am

Вот простая процедура, которая делает одну итерацию в методе вращения Якоби
для поиска собственных значений симметрической вещественной матрицы.
Сначала в двойном цикле находится максимальный внедиагональный элемент матрицы. Затем использованы формулы из вики для нахождения матрицы G .
Применяя эту процедуру последовательно несколько раз (в примере ниже 20 раз)
, получаем почти диагональную матрицу, по диагонали которой
идут приближённые собственные числа.

Jacobi:=proc(A::Matrix)
local n, a, i, j, p, c, s, G;
uses LinearAlgebra;
n:=RowDimension(A);
a:=[[0,0], 0];
for i to n do
for j from i+1 to n do
if abs(A[i,j])>=abs(a[2]) then a:=[[i,j], A[i,j]] fi;
od; od;
p:=2*a[2]/(A[a[1,1],a[1,1]]-A[a[1,2],a[1,2]]);
c:=evalf(cos(arctan(p)/2));
s:=evalf(sin(arctan(p)/2));
G:=Matrix(n,{op(subsop(a[1,1]=((a[1,1],a[1,1])=c), a[1,2]=((a[1,2],a[1,2])=c),[seq((i,i)=1, i=1..n)])), (a[1,1],a[1,2])=-s, (a[1,2],a[1,1])=s});
Transpose(G).A.G;
end proc:


Пример использования (для сравнения найдены точные собственные числа):
A:=Matrix(4, [1,2,3,4,2,5,5,6,3,5,7,7,4,6,7,8]);
(Jacobi@@20)(A);
LinearAlgebra[Eigenvalues](A);
evalf(%);

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Пт окт 05, 2018 4:04 pm

благодарю, а как увидеть матрицу G?

Kitonum
Сообщения: 2050
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Re: Метод вращения Якоби

Сообщение Kitonum » Пт окт 05, 2018 4:28 pm

Вставьте в процедуру строку print(G) после строки G:=...

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Пт окт 05, 2018 4:33 pm

увидел, спасибо большое

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Пт окт 05, 2018 4:55 pm

Код: Выделить всё

n:=10;
j:=1.5+0.1*n;
k:=n;
l:=n;
Jacobi:=proc(A::Matrix)
local n, a, i, j, p, c, s, G;
uses LinearAlgebra;
n:=RowDimension(A);
a:=[[0,0], 0];
for i to n do
for j from i+1 to n do
if abs(A[i,j])>=abs(a[2]) then a:=[[i,j], A[i,j]] fi;
od; od;
p:=2*a[2]/(A[a[1,1],a[1,1]]-A[a[1,2],a[1,2]]);
c:=evalf(cos(arctan(p)/2));
s:=evalf(sin(arctan(p)/2));
G:=Matrix(n,{op(subsop(a[1,1]=((a[1,1],a[1,1])=c), a[1,2]=((a[1,2],a[1,2])=c),[seq((i,i)=1, i=1..n)])), (a[1,1],a[1,2])=-s, (a[1,2],a[1,1])=s});
Transpose(G).A.G;
end proc:
A:=Matrix([[j,0.5*j,0,0.2*l,0],[0.5*j,j,0.3*j,0,0.1*l],[0,0.3*j,10,-0.3*j,0.5*l],[0.2*k,0,-0.3*j,j,-0.1*j],[0,0.1*k,0.5*k,-0.1*k,k]]);(Jacobi@@20)(A);
LinearAlgebra[Eigenvalues](A);
evalf(%);

у меня выводит ошибку Error, (in Jacobi) numeric exception: division by zero
у меня 10 вариант

Kitonum
Сообщения: 2050
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Re: Метод вращения Якоби

Сообщение Kitonum » Пт окт 05, 2018 9:45 pm

Проверьте исправленный вариант:
Jacobi:=proc(A::Matrix)
local n, a, i, j, p, c, s, G, r;
uses LinearAlgebra;
n:=RowDimension(A);
a:=[[0,0], 0];
for i to n do
for j from i+1 to n do
if abs(A[i,j])>=abs(a[2]) then a:=[[i,j], A[i,j]] fi;
od; od;
r:=A[a[1,1],a[1,1]]-A[a[1,2],a[1,2]];
p:=`if`(r<>0,2*a[2]/r, infinity);
c:=evalf(cos(arctan(p)/2));
s:=evalf(sin(arctan(p)/2));
G:=Matrix(n,{op(subsop(a[1,1]=((a[1,1],a[1,1])=c), a[1,2]=((a[1,2],a[1,2])=c),[seq((i,i)=1, i=1..n)])), (a[1,1],a[1,2])=-s, (a[1,2],a[1,1])=s});
Transpose(G).A.G;
end proc:

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Пт окт 05, 2018 10:10 pm

все работает, спасибо, что вы сделали?

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Пт окт 05, 2018 11:31 pm

можете объяснить что выводит a[[0,0],0] и r:=A[a[1,1],a[1,1]]-A[a[1,2],a[1,2]];

Kitonum
Сообщения: 2050
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Re: Метод вращения Якоби

Сообщение Kitonum » Сб окт 06, 2018 8:09 am

1.a это имя списка [[0,0],0] , который мы вводим для дальнейшего поиска позиции (номера строки и столбца) максимального по модулю элемента матрицы А .

2. r это выражение, стоящее в знаменателе, когда ищем нужный угол для поворота (см. формулу в вики). Например, если оказалось, что a=[[2,3],7] , то r=A[2,2]-A[3,3] , а A[2,3]=7

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Вс окт 07, 2018 10:57 pm

то есть a=[[0,0],0] третья координата это угол?

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Вт окт 16, 2018 1:13 pm

Кли писал(а):Источник цитаты то есть a=[[0,0],0] третья координата это угол?

можете ответить пожалуйста

Kitonum
Сообщения: 2050
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Re: Метод вращения Якоби

Сообщение Kitonum » Вт окт 16, 2018 4:32 pm

Нет.

Кли
Сообщения: 116
Зарегистрирован: Пт ноя 04, 2016 4:54 pm

Re: Метод вращения Якоби

Сообщение Кли » Вс окт 28, 2018 4:32 pm

Извините, но я не вижу собственные векторы