4.6. Обращение матрицы и решение системы линейных алгебраических уравнений

Рассматриваемый здесь алгоритм ите­ра­ци­он­но­го ти­па объ­единяет задачи обращения мат­ри­цы Aи ре­ше­ния сис­те­мы линейных ал­геб­ра­и­чес­ких урав­нений, за­да­ваемых в мат­ричной форме как

A^Т x = b,(1.49)

где А^Т- транспонированная матрицаA; x- век­тор­-стро­ка не­известных членов; b- вектор-строка сво­бод­ных чле­нов. В случае если заданоb = (1, 0, ..., 0), то про­це­дураinve­qu будет вычислять об­ращенную мат­ри­цуA^-1, для другихbрешается урав­нение (1.49)

Формальные параметры про­цедуры.Входные:n(типin­­te­ger) - размер матрицыА;а[1:n, 1:n] (типre­al) - ре­зуль­тат транспонирования матрицыA;b[1:n] (типreal) - век­тор-строка b;count(типin­te­ger) - параметр, за­да­ющий ко­ли­­чество итераций (в ра­­боте Роека (1962), в ко­торой при­ве­­дена пер­во­на­­чаль­ная, не оп­ти­ми­зи­ро­ван­ная версия ал­го­рит­­ма, не рекомендуется задаватьcountбольше 6);eps(типre­­al) - число, характеризующее точ­ностьпо­­лу­чен­но­­го при­бли­жения к теоретическому ре­зуль­та­ту. Вы­­ход­ные:w[1:n, 1:n] (типreal) - мат­ри­цаA^-1в слу­чаеb= (1,0,...,0);x[1:n] (типreal) - ве­к­тор-строка ре­ше­ния урав­не­ния (1.23);c[1:n,1:n] (типreal) - матрица, слу­­жащая для конт­роля точ­ности. После выполнения про­­це­ду­ры в качестве пер­вой строки матрицыcбудет век­­торbвдоль главной ее ди­а­го­нали (1,0,...,0). В ну­ле­вом цик­ле первой итерации про­це­­ду­ра invequ оп­ре­де­ляет в ка­честве строк мат­ри­цыwстол­бцы матрицыa, для всех по­сле­ду­ю­­щих итераций стро­ки матрицыw, вы­­чис­лен­­ные вn-м цикле последней ите­рации, яв­ля­ют­ся стро­­ка­ми матрицыwнулевого цикла.

procedure invequ (a:mas2; b:mas1;

n,count: integer;eps:real;

var w: mas2; x: mas1; c:mas2);

var bh,z : real; i,j,p : integer;

label beg;

begin

for j:=1 to n do

for i:=1 to n do

w[j,i]:=a[i,j];

beg:

for j:=1 to n do

for i:=1 to n do

c[i,j]:=0;

bh:=b[1];

for j:=1 to n do

begin

for i:=1 to n do

c[j,j]:=c[j,j]+w[j,i]*a[i,j];

z:=bh/c[j,j];

for i:=1 to n do

w[j,i]:=z*w[j,i];

for i:=1 to n do

if i<>j then

begin

for p:=1 to n do

c[i,j]:=c[i,j]+a[p,j]*w[i,p];

if i=1 then z:=b[j] -c[i,j];

else z:=0-c[i,j];

for p:=1 to n do

w[i,p]:=bh*w[i,p]+z*w[j,p]

end;

bh:=1

end;

z:=abs(b[1]);

for j:=1 to n do

if abs(abs(c[j,j])-z) > eps then

begin

count:=count-1;

if count > 0 then go to beg

end;

for j:=1 to n do

x[j]:=w[1,j]

end.

Процедура invequполучена из процедурыsimult[Агеев, 1976]¹путем ее перевода с языкаALGOLна языкPASCAL.

Проверка процедуры invequ производилась в це­лях контроля для тех же входных параметров, что и в ра­­бо­те Агеева (1976). При этом для вход­ных парамет­ров

,

n = 3,count = 3, 6, 10 иeps= 1e - 3, 1e - 5, 1e - 8, при всех зна­­че­­нияхcountиepsбыли получены сле­ду­ю­щие ре­зуль­та­ты, совпадающие с полученными Агеевым (1976):

1) для b= (1, 0, 0)

;

2) для b= (9, -2, 7),x= (-1, 2, 3)

.

Для матрицы, транспонированной от мат­ри­цы a, ис­­пользованной в первом тесте, результаты так­­­же хо­ро­­шо совпали с результатами, при­ве­ден­­­ны­ми М.И. Аге­­евым (1976):

1) для b= (1, 0, 0):

;

2) для b= (9, -2, 7):x= (-0.42536367e-10, 3, 2),

,

точ­ное решение уравнения (1.23), как указано Аге­е­вым (1976), при этих входных данных x = (0, 3, 2).