Главная Промышленная автоматика.

.или же для вычисления norml использовалась выборка из элементов ац, а не все недиагональные элементы, например

norml = -"yaii

4. Изменялся порядок присваивания значений элементам преобразованной матрицы, например, так, чтобы диагональные элементы после преобразования располагались примерно в порядке возрастания.

Заметим, что в процедуре }асоЫ2, приведенной в предыдущем «Замечании», все элементарные вращения являются собственными и определитель матрицы s всегда остается равным 1. Вариант, описанный в п. 4, этим свойством не обладает, и определитель матрицы s, полученной после такой процедуры, может оказаться равным -1 (т. е. ортогональная матрица может оказаться несобственной).

Замечание к алгоритму 85а В. Ф. Дейкин, Москва, май 1973 Алгоритм 85а (так же как и исходный алгоритм 85) дает правильное решение, строго говоря, не для всех симметричных матриц. Так, для тривиального случая, когда матрица а заведомо диагональна, алгоритм зацикливается.

Действительно, для указанного вида матриц .

Ш1=2 S 2Х#.П12=0,

i=2/=l

поэтому И norml = thr=sqrt(intl) =0, а талже

погт2=I {rho In) X norm 1=0. Вследствие этого условие If abs{d(p,q\)thr, находящееся в начале цикла после метки mainl, удовлетворяется всегда, и величине ind тоже всегда присваивается значение единица. Поэтому в результате проверки достижения конечной точности управление неизбежно передается на метку mainl, т. е. происходит зацикливание.

Поскольку при решении задач не всегда бывает известно, с каким видом матриц придется иметь дело, то, по-видимому, целесообразно внести в алгоритм следующие изменения.

1. После оператора

norml: = thr:=sqrt(intl); .

добавить оператор

if norml =0 then go to main2;

2. Конец процедуры

end jacobi;

заменить на

main2: end jacobi;

В этом случае на выходе из процедуры матрица а сохраняет свой первоначальный вид, а матрица s превращается в единичную.

Подтверждение к алгоритму 85 Дж. Хилмор (Н i 11 m о г е J. S. «САСМ», 1962, № 8) Оператор

omega:=(if гшОЛ then 1 else sign(mu))X(-v2)/sqrt(v2f2-l-muf2);



был заменен на оператор omega:= if mu=0.0 then-1.0 else -sign(mu) Xv2/sqrt(v2f2+mut2);

Бели mu=0, to первый оператор сводился к omega: = =-v2jsqrt{v2\2); и ошибка округления при вычислении квадратного корня могла сделать значение omega несколько большим единицы. В результате происходил аварийный останов при вычислении следующего оператора, когда делается попытка вычислить значение sqrt (1-omega\2).

В модифицированной форме алгоритм был успешно проверен с использованием транслятора Elliott ALGOL на машине National Elliott 803. Задавались матрицы до 15-го порядка, для которых получились собственные значения и собственные векторы с общей точностью равной семи десятичным разрядам.

Подтверждение к алгоритму 85

П. Наур (Н а U г Р. «САСМ», 1963, № 8)

Первоначально мы проверили алгоритм в системе GIER ALGOL с включением следующих поправок., .*

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

Эта .переделка особенно существенна в системах, имеющих сравнительно медленно работающий механизм вычисления индексов, таких как GIER ALGOL, так как она исключает более чем три из восьми ссылок на индексные переменные.

Процедура jacobi была опробована на двух различных последовательностях матриц с известными собственными значениями. В обоих случаях отладочная программа была составлена для отыскания интервала погрешностей в собственных значениях, вычисленных процедурой jacobi. Дополнительно для проверки были использованы отношения aXv-Kxv-0 (а - данная матрица, v-собственный вектор, К - соответствующее собственное значение) и а-{st)XtambdaXs=0 {s - матрица, столбцами которой служат собственные векторы, st - результат транспонирования матрицы s, а lambda-диагональная матрица собственных значений). Для отладки были использованы тест-матрицы пересмотренного алгоритма 52, приведенного в журнале «САСМ», 1963, № 1, р. 39, и матрица, предложенная Х.Б.Хансеном (Hansen Н. В.), с элементами

НВЩД=НВНМп+\-1 для ;>/,

имеющая собственное значение

0.5/(1-cos((2xi-l)Xpi/(2xn+l))).

Результаты приведены в табл. 19 и 20 (табл. 19 -для тест-матрицы Хансена; табл. 20-для тест-матрицы алгоритма 52). GIER ALGOL использует числа с плавающей запятой и 29 значащими цифрами.

* Далее указываются три поправки к алгоритму 85, учтенные при составлении-алгоритма 85а. (Прим. ред.)

** Далее приводится участок процедуры jacobi, использованный при составлении алгоритма 856. (Прим. ред.)



Абсолютные погрешности в собственных значениях

Отклонение от равенства axv - Xxa=0

Отклонение от равенства а- (st) X tambdaxs=0

Ошибка

Ошибка

Ошибка

Ошибка

Ошибка

Ошибка

rfto=l.0,0-3

-1.1,0-6

5.2,0-8

-1.7,0-4

-2.5,0-4

1.0,0-4

-7.9,0-Б

3.5,0-5

-3.3,0-3

3.0,0-3

-4.2,0-3

3.2,0-3

-9.2,0-5

3.7,0-5

-1.7,0-3

1.7.0-3

-1.5,0-3

1.8,0-3

rfto=l .0,0-5

-1.1,0-6

6.0,0-8

-1.3,0-7

4.1,0-8

-1.6,0-7

4.5,0-8

-1.2,0-5

2.2и-7

-2.7,0-5

2.2.0-5

-2.4,0-5

2.3,0-5

-3.5,0-5

3.9„-7

-6.4„-6

4.8,0-6

-5.3и,-6

4.7,0-6

rfto= 1.0,0-8

-1.1,0-6

6.0,0-8

-1.3,0-7

6.Б,о-9

-1.3,г-7

3.0,0-8

-1.2,0-5

2.2,0-7

-1.1,0-6

6.4,0-8

-5.7„-7

8.2,0-8

-3.5,0-5

3.9,0-7

-3.4,0-6

3.9,0-7

-1.3,0-6

8.9,0-8

Таблица 20

Абсолютные погрешности в собственных значениях

Отклонение от равенства axv-\xv=0

Отклонение от равенства a-(st)xlambduxs=4

Ошибка

Ошибка

Ошибка

Ошибка

Оши5ка

§•

Ошибка

UJ <§

0=1.0,0-5

-1.0,0-8

«

-З.Зю-8

4.3,0-8

-5.1,0-8

3.9,0-8

-1.1,0-8

-1.2,0-8

1.3,0-8

-5.1,0-9

2.0.О-8 1.3.0-8

-1.1,0-8

-9.3,0-9

9.4.0-9

-1.9,0-9

rfto=I.0.0-8

-7.5,0-9

3.7,0-9

-2.8,0-9

9.3,0-9

1.9,0-8

-5.6,0-9

-4.5,0-9

3.3,0-9

9.3,0-9

-1.0.0-8

-4.9,0-9

5.8,0-9

-7.5,0-9

7.5,0-9

-4.7,0-9

-2.8,0-9

3.6,0-9

-2.3,0-10

9.3,0-9

-5.1,0-9

-2.8,0-9

3.4,0-9

-1.2,0-10

7.5,0-9

-7.5,0-9

-6.0,0-9

3.2,0-9

-1.2,0-10

9.3,0-9

-4.4,0-9

- 0.0

-5.1,0-9

3.2,0-9

-7.5,ff-9

1.5,0-8

-1.5,0-8

-9.3,0-9

7.2,0-9

-2.3,0-9

2.0,0-8

-7.5,0-9

-6.5,0-9

3.0,0-9

-3.1,0-9

7.5,0-9

-5.0,0-9

-7.6,0-9 -6.9,0-9

2.4,0-9

-1.7,0-8

1.3,-8

-1.1,0-8

9.1,0-9

-3.0,0-8

3.2,0-8

-1.5,0-8

-1.1,0-8

6.7,0-9

-3.5,0-9

1.71-8

-1.1,0-8

-1.1,0-8

3.5,0-9

-3.0,0-9

7.5.0-9





0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 [22] 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

0.0018