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

Перед использованием процедуры linearsystem нужно подставить в нее вместо eps наибольшее (для данной машины) число, такое что eps + \=\.

procedure linearsystem (a,n,b,m) result: (x,det,ex,cnr);

value n,m; real det,cnr; integer n,m,ex; array a,b,x; begin real max; integer i,j,k; array lu[l:n,l:n],y,res,mult[l:n];

integer array pivot [l:n];

comment Удаляются соответствующие множители из строк матрицы а;

equilibrate(a,n,mult); спг:=1;

comment Запоминаются результаты для возможного вычисления разностей во время итерации; for i:=l step 1 until n do

for j: = l step 1 until n do lu[i,j]: = a[i,j]; . . comment Матрица преобразуется к треугольной форме; crout(]u,n,pivot,det);

comment Если не было выхода из crout к метке signall35, то вычисляется определитель в форме detXlOf ех; for i:=l step ] until n do y[i]:=lu[i,i] Xmult[i]; -

det:=det X product (y, 1 ,n,ex);

comment Теперь начинается обработка правых частей; lor к:=1 step ] until m do

begin real normy.kr; integer count.limit; kr:=k;

comment Масштабируем правые части; for i:=l step 1 until n do

res [i] :=b fi.k] :=b [i.k] /mult [i]; comment Запоминается нервое приближение иегонормаЬ(1); normy:=0;

solve (lu,n,res,pivot,y); for i:=:l step 1 until n do

begin x[i,k]:=y[i]; normy:=normy-f abs (y [i])

end i;

comment Вход в цикл итерации. Итерация заканчивается по достижении счетчиком count целого значения limit, которое определяется по результатам первой итерации и по значению зависящего от машины вещественного числа eps; for count:=l,2 step 1 until limit do begin real t;

comment Вычисляются разности решения у;

residuals (a,n,b,k,x,res);

comment Находится очередное приращение в решении; solve (lu,n,res,pivot,y);

comment Устанавливаются условия окончания; if count=l then begin real normdy; normdy:=0;

for i:=l step 1 unti! n do normdy: = normdy-b abs (y [i]); - if normdy=0 then go to enditer;

t: = normy/normdy;



comment Вычисляется IIaXla-MI (спектральная норма) - мера обусловленности матрицы а. Она является мерой затруднений в нахождении решения исходного уравнения и входит, естественно, в границы погрешно сти решения. Параметр спг является прямой мерой погрешности и экспериментально аппроксимирует меру обусдовденности;

if t<2 then go to signal 135; . •

cnr:= ((kr-1) X СПГ+1 / (eps Xt))/кг; limit:.= ln(eps)/ln(l/t) end if;

comment Запоминается новое приближение; for i: = l step 1 until n do x[i,k] :=x[i,k] --y[i] end count; enditer: end к end linearsystem;

Процедура equilibrate {equilibrate - уравновешивать) выбирает максимальный по модулю элемент в каждой строке матрицы а[1:/г, 1:п], запоминает его в массиве mtiZ/[l:n] и делит на него все элементы соответствующей строки.

Если матрица плохо обусловлена, то решение очень чувствительно по отношению к возмущению входных данных, и тогда уравновешивающее деление нужно делать не на максимальный элемент, а на степень основания машинной системы счисления (2 или 10 для двоичных или десятичных машин соответственно), ближайшую к максимальному элементу, для того чтобы исключить ошибки округления.

Метод уравновешивания рассматривается в работе Дж. Вилкинсона [5i].

procedure equilibrate (a,n,mult);

value n; integer n; array a.mult; begin real max; integer i,j;

for i:=l step 1 until n do begin max:-0;

for j:=l step 1 until n do

if abs(a[i,j]) >max then max:=abs(a[i,j]).; , if max=0 then go to signall35;

mult[i]:=max;

for j:=l step 1 until n do a[i,j] :=a[i,j]/max end i end equilibrate;

Процедура crout выполняет метод Краута с вращением, как это сформулировано в алгоритмах Г. Форсайта (Forsythе G. Е. Algorithm 16. «САСМ», 1960, № 9 и алгоритм 436 [23]), для преобразования матрицы а[1:п, l:n] в треугольное представление LU со всеми L[,]=l. В массиве pivotyil запоминается индекс строки вращения на -м шаге исключения для использования в процедуре solve, sg. (сокращение от signum--знак) получает значение ±1 и меняет знак с каждой перестановкой строк матрицы а.

procedure crout (a,n,pivot,sg);

value n; integer n,sg; array a; integer array pivot;



begin real t.quot; integer i,j,k,imax,p; real procedure ipl (a,t,f,i,k);

value t,f,i,k; real t; integer f,i,k; array a;

comment Программист может получить значительный выигрыш, подставив вместо процедуры ipl более быструю и более точную процедуру скалярного произведения; begin real sum; integer p; sum:=t;

for p:=l step 1 until f do sum:=sum-f a[i,p] Xa[p,k]; ipl:=sum end ipl;

sg:=l; comment к - индекс шага исключения; for к:=1 step ] until n do begin t:=0;

for i:=k step 1 until n do

begin comment Вычисление L; . afi,k]:=-ipl(a,-a[i,k],k-l,i,k); if abs(a[i,k])>t then begin t: = abs(a[i,k]); imax: = i end end i;

if t=0 then go to signall35;

comment a[imax,k]-наибольший элемент из оставшихся в столбце к. Далее переставляем строки, если нужно, и фиксируем факт перестановки; pivot[k]: = imax; .if imaxk then begin sg:=-sg;

for j:=l step 1 until n do begin t:=a[k,j];

a[k,j]:=a[imax,j]; a[imax,j]:=t end j

end Затем вычисляем столбец множителей; quot:=l/a[k,k];

for i:=:k-bl step 1 until n do a[i,k]:=a[i,k] Xquot; comment Далее вычисляем строку U; for j:=k--l step 1 until n do

a [k,j] :,=-ipl (a,-a [k,j],k-l,k,j) end к end crout;

Процедура-функция product {product - произведение) принимает значение произведения компонент вектора lactors\s:f {factor - множитель, S и / - начальные буквы слов start и finisJi), предохраняя от пе-оеполнения и обращения произведения в машинный нуль тем, что результат каждого умножения нормализуется (т. е. 0.1результат<1), а каждая компонента вектора factors, меньшая чем 0.1, перед умножением увеличивается в 10 раз. Показатель порядка ех произведения соответственно корректируется.

real procedure product (factors,s,f) result: (ex);

value s,f; integer s,f,ex; array factors; begin real p,pl; integer i;

ex:=0; p:=l;





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.0033