БЧА НИВЦ МГУ. ASP1R_P. Решение систем с прямоугольными матрицами полного ранга

Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/as_htm_p/asp1r_p.htm
Дата изменения: Fri Dec 4 10:45:24 2015
Дата индексирования: Sun Apr 10 02:40:53 2016
Кодировка: Windows-1251

Текст подпрограммы и версий
asp1r_p.zip , asp1e_p.zip , asp1c_p.zip

Тексты тестовых примеров
tasp1r_p.zip , tasp1e_p.zip , tasp1c_p.zip

Подпрограмма: ASP1R (модуль ASP1R_p)

Назначение

Нахождение нормального решения неопределенной системы линейных алгебраических уравнений полного ранга методом отражений.

Математическое описание

Находится нормальное решение системы АХ = В, где А - матрица полного ранга размера N на М (N ≤ М), В - заданный вектор длины N. Для решения используется нормализованное приведение матрицы системы к нижней треугольной форме L с помощью последовательности преобразований отражения

 (1)     SAQ₁Q₂...Q_N = L ,

где Q_i - соответствующие матрицы отражения, 1 ≤ i ≤ N , S - результирующая матрица перестановок, L - нижняя треугольная матрица размера N на М.

Из полученной треугольной системы

        LY = SB

находится нормальное решение Y, по которому затем определяется искомое решение Х = Q₁Q₂ ... Q_NY .

Д.К.Фадеев, В.Н.Кублановская, В.Н.Фадеева, О решении линейных алгебраических систем с прямоугольными матрицами. Тp. Мат. ин-та АН СССР, 1968, 96.

Использование

procedure ASP1R(var A :Array of Real; var B :Array of Real;
                var X :Array of Real; var T :Array of Real;
                var S :Array of Integer; N :Integer; M :Integer;
                L :Integer);

   Параметры

A -	вещественный двумерный массив размера N на М (N ≤ М), в котором задается исходная матрица А; в результате работы подпрограммы в массиве А на соответствуюших местах запоминаются поддиагональные элементы вычисленной нижней треугольной матрицы L; в остальной части массива в последовательных строках запоминаются векторы, порождающие соответствующие матрицы отражения;
B -	вещественный вектор длины N, в котором задается правая часть системы;
X -	вещественный вектор длины М , в котором запоминается найденное решение системы;
T -	вещественный вектор длины М, используемый подпрограммой как рабочий; в результате работы подпрограммы в первых N компонентах массива Т запоминаются диагональные элементы нижней треугольной матрицы L;
IS -	целый вектор длины N, в котором запоминаются перестановки строк при проведении нормализованного процесса (1); в IS (k) запоминается номер строки, переставленной с k - ой строкой на k - ом шаге преобразования;
N, M -	число строк и столбцов исходной матрицы А, N ≤ М (тип: целый);
L -	задает режим работы подпрограммы (тип: целый); при этом:

L = 1 -	если система с данной матрицей А решается в первый раз;
L ≠ 1 -	если система решается повторно с той же матрицей А, но с другой правой частью; в этом случае разложение (1) матрицы А не выполняется.

   Версии

ASP1E -	нахождение нормального решения недоопределенной системы линейных алгебраических уравнений полного ранга, заданной с расширенной (Extended) точностью, методом отражений.
ASP1C -	нахождение нормального решения комплексной недоопределенной системы линейных алгебраических уравнений полного ранга методом отражений.

   Вызываемые подпрограммы:  нет

   Замечания по использованию

	1.	В подпрограмме АSР1С параметры А, В, Х, Т имеют тип Complex.
	2.	В подпрограмме АSР1E параметры А, В, Х, Т имеют тип Extended.
	3.	При повторном решении системы с той же матрицей информация, полученная ранее в массивах А, Т, IS не должна портиться.

Пример использования

Unit TASP1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, ASP1R_p;

function TASP1R: String;

implementation

function TASP1R: String;
var
P,I,J,J1,I1,_i :Integer;
X :Array [0..4] of Real;
T :Array [0..4] of Real;
S :Array [0..3] of Integer;
A1 :Array [0..19] of Real;
B1 :Array [0..3] of Real;
const
A :Array [0..19] of Real = ( 1.0,-2.0,3.0,-4.0,-5.0,6.0,-7.0,8.0,9.0,10.0,
-11.0,12.0,0.0,1.0,-1.0,4.0,-5.0,-15.0,16.0,-20.0 );
B :Array [0..11] of Real = ( 1.0,-24.0,27.0,-36.0,6.0,13.0,-13.0,16.0,6.0,
-8.0,10.0,-12.0 );
label
_100,_110,_101,_102;
begin
Result := '';  { результат функции }
for I:=1 to 4 do
 begin
  for J:=1 to 5 do
   begin
_100:
    A1[(I-1)+(J-1)*4] := A[(I-1)+(J-1)*4];
   end;
 end;
J1 := 0;
for I1:=1 to 3 do
 begin
  for I:=1 to 4 do
   begin
_110:
    B1[I-1] := B[J1+I-1];
   end;
  Result := Result + Format('%s',
 ['          РЕШЕНИЕ НЕДООПРЕДЕЛЕННОЙ СИСТЕМЫ AX=B ПОЛНОГО PAHГA' + #$0D#$0A +
 '          МЕТОДОМ ОТРАЖЕНИЙ ДЛЯ ВЕЩЕСТВЕННЫХ A И B.' + #$0D#$0A]) + #$0D#$0A; 
  Result := Result + Format('%s',
 [' ИСХОДНЫЕ ДAHHЫE:' + #$0D#$0A]) + #$0D#$0A; 
  Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' A' + #$0D#$0A]);
  for I:=1 to 4 do
   begin
    for J:=1 to 5 do
     begin
      Result := Result + Format('%20.16f ',
 [A[(I-1)+(J-1)*4]]) + #$0D#$0A;
     end;
   end;
  Result := Result + #$0D#$0A;
  Result := Result + Format('%s',[' B' + #$0D#$0A]);
  Result := Result + #$0D#$0A;
  for _i:=0 to 3 do
   begin
    Result := Result + Format('%20.16f ',[B1[_i]]);
    if ( ((_i+1) mod 4)=0 )
     then Result := Result + #$0D#$0A;
   end;
  Result := Result + #$0D#$0A;
  Result := Result + Format('%s',[' PEЗYЛЬTAT:' + #$0D#$0A]) + #$0D#$0A; 
  P := 1;
  ASP1R(A,B1,X,T,S,4,5,P);
  Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' A' + #$0D#$0A]);
  for I:=1 to 4 do
   begin
    for J:=1 to 5 do
     begin
      Result := Result + Format('%20.16f ',
 [A[(I-1)+(J-1)*4]]) + #$0D#$0A;
     end;
   end;
  Result := Result + #$0D#$0A;
  Result := Result + Format('%s',[' T' + #$0D#$0A]);
  Result := Result + #$0D#$0A;
  for _i:=0 to 4 do
   begin
    Result := Result + Format('%20.16f ',[T[_i]]);
    if ( ((_i+1) mod 4)=0 )
     then Result := Result + #$0D#$0A;
   end;
  Result := Result + #$0D#$0A;
  Result := Result + Format('%s',[#$0D#$0A + ' S' + #$0D#$0A]);
  Result := Result + #$0D#$0A;
  for _i:=0 to 3 do
   begin
    Result := Result + Format('%10d ',[S[_i]]);
    if ( ((_i+1) mod 4)=0 )
     then Result := Result + #$0D#$0A;
   end;
  Result := Result + #$0D#$0A;
  Result := Result + Format('%s',[#$0D#$0A + ' B' + #$0D#$0A]);
  Result := Result + #$0D#$0A;
  for _i:=0 to 3 do
   begin
    Result := Result + Format('%20.16f ',[B1[_i]]);
    if ( ((_i+1) mod 4)=0 )
     then Result := Result + #$0D#$0A;
   end;
  Result := Result + #$0D#$0A;
  Result := Result + Format('%s',[#$0D#$0A + ' X' + #$0D#$0A]);
  Result := Result + #$0D#$0A;
  for _i:=0 to 4 do
   begin
    Result := Result + Format('%20.16f ',[X[_i]]);
    if ( ((_i+1) mod 4)=0 )
     then Result := Result + #$0D#$0A;
   end;
  Result := Result + #$0D#$0A;
  for I:=1 to 4 do
   begin
    for J:=1 to 5 do
     begin
_101:
      A[(I-1)+(J-1)*4] := A1[(I-1)+(J-1)*4];
     end;
   end;
  J1 := J1+4;
_102:
 end;
UtRes('TASP1R',Result);  { вывод результатов в файл TASP1R.res }
exit;
end;

end.


Результат:

      X  =   ( 0.00000000044,  -0.069999999928,  0.799999999989,  -2.7000000002,
                  3.4999999997,  -1.5400000002 )