Главная Промышленная автоматика. Перед использованием процедуры 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.0017 |