|
DirectSolverN Библиотека "JINRLIB" F099 
(N=2,3,4,5,6)
Автор: А.П.Сапожников
Язык: Фортран
РЕШЕНИЕ СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ ОБЩЕГО ВИДА
ПРЯМЫМ МЕТОДОМ (ПО ПРАВИЛУ КРАМЕРА)
Подпрограмма-функция DirectSolverN (N=2,3,4,5,6) выполняет
решение системы линейных уравнений A*X = B общего вида
прямым методом, т.е. по правилу Крамера: X(i)=DX(i)/D
Общеизвестно, что для решения системы порядка N методами типа
Гаусса-Жордана требуется O(N**3) операций.
В то же время прямое использование правила Крамера требует
вычисления N+1 определителя N порядка, т.е. (N+1)! операций.
Однако при малых N обе оценки вполне соизмеримы, а при N<4
просто-напросто (N+1)! < N**3.
Более того, экспериментально установлено, что вплоть до 6 порядка
система решается прямым методом гораздо быстрее (при N=6 - в 2 раза).
Поэтому мы предоставляем всем желающим семейство программ
DirectSolverN(A,B,X) (N=2,3,4,5,6)
для решения систем линейных уравнений A*X=B, где
A - матрица системы NxN, B - массив длиной N,
X - массив длиной N, куда будет записано решение.
Все программы оформлены как функции, возвращающие детерминант
матрицы системы или 0, если система несовместна.
Особенностью этого семейства программ является то, что они
генерировались автоматически с помощью специально построенной
для этой цели рекурсивной процедуры DeterGen, написанной на
языке Паскаль.
В самом деле - кому захочется вручную "расписывать" определитель
6 порядка! Да и ошибиться тут нетрудно. Текст упомянутой
процедуры публикуется ниже, мы полагаем, что он представляет
самостоятельную ценность.
Замечание:
----------
В связи с ограничениями фортранного компилятора максимальное значение
размерности системы N для Microsoft Fortran 5.00 равно 5.
Пример:
-------
Implicit Real*8 (a-h,o-z)
dimension a(4,4),b(4),x(4)
do i=1,4
b(i)=0.0d0
do j=1,4
a(i,j)=1.0d0/(i+j) ! Hilbert matrix
b(i)=b(i)+a(i,j) ! solution = all 1
enddo
enddo
d=DirectSolver4(a,b,x)
write(*,*) ' Determinant for Hilbert matrix = ',d
Результат:
----------
Determinant for Hilbert matrix = 2.362055933475061E-009
Текст процедуры DeterGen:
-------------------------
Procedure DeterGen(N:integer; Lines:TStrings);
Type Matrix=array[1..9,1..9] of string[6];
var i,j,k,N0:integer;
a:Matrix; Indent:string;
Function Dig(d:integer):char;
begin Result:=Chr(d+Ord('0')); end;
Procedure Determ(n:integer; a:Matrix; Lines:TStrings);
var s,s1,sg:string; i,j,k:integer;
b:Matrix;
begin sg:=''; s:='D'+Dig(n);
if (n=2) and (N0=2) then begin Lines.Add(Indent+'D2 = '+a[1][1]+'*'+a[2][2]+'-'+a[1][2]+'*'+a[2][1]); Exit; end;
if n=3 then
begin s1:=Indent+' '; s1[6]:='&';
Lines.Add(Indent+'D3 = '+a[1][1]+'*('+a[2][2]+'*'+a[3][3]+'-'+a[3][2]+'*'+a[2][3]+')-'+a[1][2]+'*('+a[2][1]+'*'+a[3][3]+'-');
Lines.Add(s1+a[2][3]+'*'+a[3][1]+')+'+a[1][3]+'*('+a[2][1]+'*'+a[3][2]+'-'+a[2][2]+'*'+a[3][1]+')');
Exit;
end;
for j:=1 to n do
begin
for i:=1 to n-1 do for k:=1 to n do b[i][k]:=a[i+1][k];
for i:=1 to n-1 do for k:=j to n-1 do b[i][k]:=b[i][k+1];
Determ(n-1,b,Lines);
s1:=Indent+s+' = '; if j<>1 then s1:=s1+s+sg;
Lines.Add(s1+a[1][j]+' * D'+Dig(n-1));
if sg=' - ' then sg:=' + ' else sg:=' - ';
end;
end;
begin N0:=N; Indent:=' ';
Lines.Clear; if n>9 then n:=9;
Lines.Add(' Function DirectSolver'+Dig(n)+'(A,B,X) ! for Linear System A*X=B');
Lines.Add('C *** Generated by "DeterGen" procedure ***');
Lines.Add(' Implicit Real*8 (A-H,O-Z)');
Lines.Add(' Dimension A('+Dig(n)+','+Dig(n)+') ! system Matrix');
Lines.Add(' Dimension B('+Dig(n)+') ! system Right Part');
Lines.Add(' Dimension X('+Dig(n)+') ! system Solution');
for i:=1 to n do for j:=1 to n do a[i][j]:='A('+Dig(i)+','+Dig(j)+')';
Determ(n,a,Lines);
Lines.Add(' D = D'+Dig(n));
Lines.Add(' if(D.ne.0) then');
for k:=1 to n do
begin Indent:=' ';
for i:=1 to n do a[i][k]:='B('+Dig(i)+')';
Determ(n,a,Lines);
Lines.Add(Indent+'X('+Dig(k)+') = D'+Dig(n)+' / D');
for i:=1 to n do a[i][k]:='A('+Dig(i)+','+Dig(k)+')';
end;
Lines.Add(' endif');
Lines.Add(' DirectSolver'+Dig(n)+'=D');
Lines.Add(' return');
Lines.Add(' End');
end;
|
