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

for i:=:s step 1 until i do begin pl:=factors [i]; if abs(pl)<0.1 then

begin pl: = 10Xpl; ex:=ex-l end; p:=pXpl; if p=0 then

begin ex:=0; go to fin end; numl: if abs(p) <0.1 then

begin p:=pX10; ex:=;ex-1; go to numl end; num2: if abs (p) 1 then

begin p:=0.1 Xp; ex:=ex+1; go to num 2 end end i; fin: product:,=p end product;

Процедура residuals (residuals - разности) вычисляет столбец разностей с-аХу, где с - к-й столбец матрицы правых частей Ь, а у- k-и столбец матрицы х.

procedure residuals (a,n,b,k,x) result: (res);

value n,k; integer n,k; array a,b,x,res; begin integer i;

real procedure ip2(a,i,n,x,k,t);

value i,n,k,t; real t; integer i,n,k; array a,x;

comment Процедура ip2 получает значение скалярного произ ведения i-й строки матрицы а на k-й столбец матрицы решений X, добавляя к нему слагаемое t. Важно, чтобы ip2 была «накапливающей» или пользовалась арифметикой с двойной точностью, как об этом говорится в работе Дж. Вилкинсона [5i];

for i:=:=i step I until n do res[i]:=-ip2(a,i,n,x,k,-b[i,k]) end residuals;

Процедура solve обрабатывает правые части b и затем находит обратное решение у, используя представление LU, предусмотренное методом Краута.

procedure solve (a,n,b,pivot) result: (у);

value n; integer n; array a,b,y; integer array pivot; begin real t; integer k,p;

for k:=l step 1 until n do begin t:z=:b[pivot[k]]; b [pivot[k]]:=b[k];

for p:=l step 1 until k-1 do t:=t-a[k,p] Xb[p]; b[k]:=t

end модификации b; .

for k:=n step -1 until 1 do begin t:=b[k];

for p:=k-l-1 step 1 until n do t:=t-a[k,p] Xy[p]; y[k]:=t/a[k,k] end к end solve;



Левые части

Правые части

2x+2y+2z=

2x+2y+3z=

При этом в процедуру linearsystem задавались параметры m=n=3,

В результате было получено

, detOA, ел=1, СП/-=0.64028427.

4 2 2 2 2 2 2 2 3

2-13

3 1 2

4 2 3

- 0.5 -1 0.5 1 0.5-0.5

1 1 1

При этом процедура ip2 (в теле процедуры residuals) имела вид

real procedure ip2(a,i,n,x,k,t);

value i,n,k,t; real t; integer i,n,k; array a,x; begin real sum; integer p; sum:=t;

for p:=l step 1 until n do sum:=sum-l-a[i,p] Xx[p,k]; ip2:=sum end ip2;

В тело процедуры вместо eps и ln{eps) были подставлены значения 10-10 и -23.02585 соответственно.

2. Для вычисления определителя с помощью алгоритма 1356 в нем задавались параметры т=1 и Ь=(1,0,0). Остальные параметры были такими же, как и в предыдущей задаче. Был получен правильный результат det=OA и ех=1.

3. Для обращения матрицы задавались параметры т=п=3;

* с такой же системой уравнений отлаживался алгоритм 197а. См. свидетельство к алгоритму 197а. (Прим. ред.)

Алгоритм 1356 ло существу не отличается от алгоритма 135а, если не считать того, что в процедуре residuals параметр t был внесен в список значений.

Алгоритм 1356 был транслирован в системе ТА-1М на мащине М-220, и с его помощью были рещены следующие задачи.

1. Система трех уравнений с тремя вариантами правых частей*



Был получен результат

0.5 -0.5 - 0.5 2

О -1

- 1 1

В правильности которого легко убедиться умножением х на с, дающим единичную матрицу.

Свидетельство к алгоритму 135а

Алгоритм 135а получен в результате ординарной переработки и некоторых очевидных сокращений алгоритма 135 (МсКеетап W. М. «САСМ», 1962, № 11), а также внесения в него поправок, указанных в «Замечании к алгоритму 135» В. М. Мак-Кимена (МсКеетап W. М. «САСМ», 1964, № 7).

Ошибка в процедуре crout, указанная Л. Мейснером в его «Заме-, чании к алгоритму 135» (см. перевод ниже), была исправлена путем внесения индексов i и k в список формальных параметров процедуры crout.

Кроме того, в процедуре residuals была обнаружена и исправлена ошибка, заключающаяся в том, что описание процедуры ip2 было помещено не в начале тела процедуры residuals, а в совокупности ее спецификаций.

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

В. М. .Мак-Кимен (МсКеетап W. М. «САСМ», 1962, № II)

Алгоритм был проверен на точность и на время вьшолнения на машине Burroughs 220 с транслятором с языка BALGOL. Точный обращенный сегмент матрицы Гильберта (порядка 6) был записан в ячейках машины В 220 с восемью десятичными цифрами и с плавающей запятой и был использован в тестах. Сегмент матрицы Гильберта Не очень плохо обусловлен (спектральная норма ЯеХЯе--1.3Х 10). Следовательно, число требующихся итераций не может считаться типичным.

Таблица 19

Вид решения

„Точное" уравновешивание (по степени основания 10)

Уравновешивание по максимальному элементу в строке

Начальное решение Итерация

первая

вторая

третья

четвертая

пятая

0.092587535

0.094091506

0.090877240 0.090909695 0.090909080 0.090909091 Окончание

0.091498265 0.091570311 0.091568310 0.091568365 0.091568364 Окончание

В табл. 19 приводятся значения элемента с индексами [п,п] (математически равного 0.090909...). Эти значения могут характеризовать поведение остальных элементов.





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

0.0017