Текст подпрограммы и версий
aej1r_p.zip , aej1e_p.zip
Тексты тестовых примеров
taej1r_p.zip , taej1e_p.zip

Подпрограмма:  AEJ1R (модуль AEJ1R_p)

Назначение

Вычисление всех собственных значений и соответствующих собственных векторов матрицы Якоби неявным Q*L - алгоpитмом.

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

В подпрограмме используется следующее свойство матрицы Якоби: произведения пар соответствующих элементов побочных диагоналей все неотрицательны и равны нулю только, если оба сомножителя равны нулю - для ее симметризации диагональными подобными преобразованиями, после чего вычисляются собственные значения и собственные векторы исходной матрицы с помощью неявного Q*L - алгоритма.

Дж.Х.Уилкинсон, Алгебраическая проблема собственных значений, "Наука", М., 1970.

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

procedure AEJ1R(N :Integer; var A :Array of Real;
                var EV :Array of Real; var V :Array of Real;
                var RAB1 :Array of Real; var IERR :Integer);

Параметры

N - порядок исходной матрицы (тип: целый);
A - вещественный двумерный массив размерности N на 3, содержащий в последних N - 1 компонентах первого столбца элементы нижней диагонали, во втором столбце - элементы главной диагонали, в первых N - 1 компонентах третьего столбца - элементы верхней диагонали;
EV - вещественный одномерный массив размерности N, содержащий вычисленные собственные значения в возрастающем порядке;
V - вещественный двумерный массив размерности N на N, содержащий вычисленные ортонормированные собственные векторы исходной матрицы;
RAB - вещественный одномерный массив размерности N, используемый как рабочий;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом значение IЕRR:
  а) больше N, если не все попарные произведения соответствующих элементов побочных диагоналей неотрицательны;
  б) больше 2*N, если имеется равное нулю произведение соответствующих элементов побочных диагоналей, причем сомножители равны нулю неодновременно. В этом случае нет симметризации, необходимой для правильного вычисления собственных векторов;
  в) полагается равным номеру собственного значения для вычисления которого потребовалось более 30 итераций, при этом вычислены правильно, но неупорядочены, собственные значения с индексами 1, 2, ..., IЕRR - 1 и соответствующие собственные векторы.

Версии

AEJ1E - вычисление всех собственных значений и соответствующих собственных векторов матрицы Якоби с расширенной (Extended) точностью. Массивы А, ЕV, V, RАВ имеют тип Extended.

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

UTAE10 - подпрограмма выдачи диагностических сообщений при работе подпрограмм АЕJ1R и АЕJ1E.

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

  Подпрограма АЕJ1R сохраняет исходную матрицу.

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

Unit TAEJ1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AEJ1R_p;

function TAEJ1R: String;

implementation

function TAEJ1R: String;
var
J,I,_i,IERR :Integer;
EV :Array [0..4] of Real;
RАВ :Array [0..4] of Real;
V :Array [0..24] of Real;
const
A :Array [0..14] of Real = ( 0.0,1.0,1.0,1.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,
1.0,1.0,0.0 );
begin
Result := '';  { результат функции }
Result := Result + #$0D#$0A;
Result := Result + Format('%s',
 [' ВСЕ СОБСТВЕННЫЕ ЗНАЧЕНИЯ И СОБСТВЕННЫЕ ВЕКТОРЫ ВЕЩЕСТВЕННОЙ  МАТРИЦЫ ЯKOБИ' + #$0D#$0A +
 ' A' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
AEJ1R(5,A,EV,V,RAB,IERR);
Result := Result + Format('%s',[' PEЗYЛЬTAT']) + #$0D#$0A; 
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' A' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' EV' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 4 do
 begin
  Result := Result + Format('%20.16f ',[EV[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' V' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  for J:=1 to 5 do
   begin
    Result := Result + Format('%20.16f ',[V[(I-1)+(J-1)*5]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' IERR' + #$0D#$0A]);
Result := Result + Format('%8d ',[IERR]) + #$0D#$0A;
UtRes('TAEJ1R',Result);  { вывод результатов в файл TAEJ1R.res }
exit;
end;

end.

Результаты:

      EV(1)  =  -1.7320508 , 
      EV(2)  =  -1.0 , 
      EV(3)  =  6.28*10-12 , 
      EV(4)  =  1.0 , 
      EV(5)  =  1.7320508 ;

      V1    =  (0.2887,  -0.5,  0.5774,  -0.5,  0.2887) , 
      V2    =  (-0.5,  0.5,  -1.7*10-12,  -0.5,  0.5) , 
      V3    =  (-0.5774,  2.625*10-12,  0.5774,  4.*10-12, -0.5774) , 
      V4    =  (-0.5,  -0.5,  -2.035*10-12,  0.5,  0.5) , 
      V5    =  (-0.2887,  -0.5,  -0.5774,  -0.5,  -0.2887) ;

      IERR  =  0