Версия для печати
Нажмите сюда для просмотра этой темы в оригинальном формате |
Форум на Исходниках.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 |
Скорее всего, при нулевом условии угол должен быть нулевым и ничего не повернётся, всё останется прежним. А у вас как-то обратно выходит. Надо посмотреть и избавиться от неопределённости. Эх, вспоминать надо. Добавлено Хм, у вас прям как в вики. Ну тогда надо использовать atan2(y,x) и всё будет хорошо! Добавлено Ой, пардон - у вас же какой-то необычный язык... ну тогда лучше всего рассмотреть два случая и вызывать arctg и arcctg в каждом. |
Автор: Simon25 06.01.19, 16:15 |
это maple 17 |
Автор: amk 06.01.19, 19:59 |
Цитата Славян @ Не надо. Сам угол не нужен. Нужны тангенс и синус этого угла. Они вычисляются по значениям найденного внедиагонального элемента и соответствующих диагональных.Хм, у вас прям как в вики. Ну тогда надо использовать atan2(y,x) и всё будет хорошо! Но что помню, там на каждом цикле два раза приходится квадратный корень искать (или это когда синус и косинус используются?). Точнее, нужна функция hypot. Да, неплохо было бы язык, на котором это всё написано, указать. Чтобы отвечающие могли возникающие непонятки для себя прояснить. atan2 в нём скорее всего тоже есть. |