Версия для печати
Нажмите сюда для просмотра этой темы в оригинальном формате
Форум на Исходниках.RU > Прочие языки программирования > Собственные значения и векторы


Автор: Simon25 06.01.19, 11:59
Здравствуйте!!! Написал код, для поиска собственных значений и векторов методом вращений Якоби
но препод сказал, что он коряв с точки зрения программирования.
Сможете ответить, в чём проявляется корявость?
Скрытый текст
<{CODE_COLLAPSE_OFF}><{CODE_WRAP_OFF}>
    restart;
    with(linalg):
    n:=10;
    N:=5;
    j:=1.5+0.1*n;
    k:=n;
    l:=n;  
    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*j,j]]);
    eigenvectors(A);
    B:=A: k:=0: flag:=0: eps:=0.001:
    while (flag = 0) do
    m:=abs(B[1,2]):
    iStr:=1: iCol:=2:
    for i from 1 to N do
    for j from 1 to N do
    if (i < j) and (j > 2) then
    if (m < abs(B[i,j])) then
    m:=abs(B[i,j]):
    iStr:=i: iCol:=j:
    end if: end if: end do: end do:
    if abs(m) <= eps then
    flag:=1:
    else
    phi:=evalf((1/2)*arctan(2*B[iStr,iCol]/(B[iStr,iStr]-B[iCol,iCol]))):
    H[k]:=diag(1,1,1,1,1):
    H[k][iStr,iStr]:=cos(phi):
    H[k][iStr,iCol]:=-sin(phi):
    H[k][iCol,iStr]:=sin(phi):
    H[k][iCol,iCol]:=cos(phi):
    if k = 0 then
    G:=H[0]:
    end if:
    #print(B);
    #print(H[k]);
    B:=multiply(transpose(H[k]), multiply(B,H[k])):
    if k <> 0 then
    G:=multiply(G,H[k]):
    end if:
    k:=k+1:
    end if: end do:
    #print(k);
    G:=multiply(seq(H[i],i=0..k-1)):
    for i from 1 to N do
    lambda[i]:=B[i,i];
    e[i]:=col(-G,i);
    end do;

Автор: Славян 06.01.19, 13:12
Из недостатков я бы отметил такие:
1. Всё без разметки - тяжело читать (у меня так, по крайней мере, отображается).
2. Я так и не понял где гарантия, что при арктангенсе на нуль не поделим?
3. Синус и косинус дважды считается (мелочь, впрочем).

Автор: Simon25 06.01.19, 13:18
Цитата Славян @
тяжело читать

а как подкорректировать код, чтобы был легко читаем

Автор: Славян 06.01.19, 13:44
Ну пробелы там порасставлять, и в таком стиле:
<{CODE_COLLAPSE_OFF}><{CODE_WRAP_OFF}>
    if A
         if B
         end if
    end if

Автор: Simon25 06.01.19, 13:50
Цитата Славян @
где гарантия, что при арктангенсе на нуль не поделим?

там условия if нужно?

Добавлено
так правильно?
<{CODE_COLLAPSE_OFF}><{CODE_WRAP_OFF}>
    r:=B[iStr,iStr]-B[iCol,iCol]:
    if (r<>0) then
    phi:=evalf((1/2)*arctan(2*B[iStr,iCol]/r)):
    end if:

Автор: amk 06.01.19, 15:00
Насколько помню, в методе Якоби ищется максимальный по модулю внедиагональный элемент. Пока он ненулевой, то тангенс угла поворота тоже ненулевой. А когда он становится нулевым повороты прекращаются.

Автор: Simon25 06.01.19, 15:12
совершенно правильно
значит мое условие правильно?

Автор: Славян 06.01.19, 16:00
Скорее всего, при нулевом условии угол должен быть нулевым и ничего не повернётся, всё останется прежним. А у вас как-то обратно выходит. Надо посмотреть и избавиться от неопределённости.
Эх, вспоминать надо. :blush:

Добавлено
Хм, у вас прям как в вики. Ну тогда надо использовать atan2(y,x) и всё будет хорошо! ;)

Добавлено
Ой, пардон - у вас же какой-то необычный язык... ну тогда лучше всего рассмотреть два случая и вызывать arctg и arcctg в каждом.

Автор: Simon25 06.01.19, 16:15
это maple 17

Автор: amk 06.01.19, 19:59
Цитата Славян @
Хм, у вас прям как в вики. Ну тогда надо использовать atan2(y,x) и всё будет хорошо!
Не надо. Сам угол не нужен. Нужны тангенс и синус этого угла. Они вычисляются по значениям найденного внедиагонального элемента и соответствующих диагональных.
Но что помню, там на каждом цикле два раза приходится квадратный корень искать (или это когда синус и косинус используются?). Точнее, нужна функция hypot.

Да, неплохо было бы язык, на котором это всё написано, указать. Чтобы отвечающие могли возникающие непонятки для себя прояснить.
atan2 в нём скорее всего тоже есть.

Powered by Invision Power Board (https://www.invisionboard.com)
© Invision Power Services (https://www.invisionpower.com)