Собственные значения и векторы

Здравствуйте!!! Написал код, для поиска собственных значений и векторов методом вращений Якоби
но препод сказал, что он коряв с точки зрения программирования.
Сможете ответить, в чём проявляется корявость?

Скрытый текст

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;

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

Цитата Славян @ 06.01.19, 13:12

тяжело читать

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

Ну пробелы там порасставлять, и в таком стиле:

if A
     if B
     end if
end if

Сообщение отредактировано: Славян - 06.01.19, 13:45

Цитата Славян @ 06.01.19, 13:12

где гарантия, что при арктангенсе на нуль не поделим?

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

Добавлено 06.01.19, 14:14
так правильно?

r:=B[iStr,iStr]-B[iCol,iCol]:
if (r<>0) then
phi:=evalf((1/2)*arctan(2*B[iStr,iCol]/r)):
end if:

Сообщение отредактировано: Simon25 - 06.01.19, 14:19

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

совершенно правильно
значит мое условие правильно?

Сообщение отредактировано: Simon25 - 06.01.19, 15:13

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

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

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

это maple 17

Цитата Славян @ 06.01.19, 16:00

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

Не надо. Сам угол не нужен. Нужны тангенс и синус этого угла. Они вычисляются по значениям найденного внедиагонального элемента и соответствующих диагональных.
Но что помню, там на каждом цикле два раза приходится квадратный корень искать (или это когда синус и косинус используются?). Точнее, нужна функция hypot.

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