БЧА НИВЦ МГУ. AEJ3R_P. Линейная проблема собственных значений для матриц специального вида

Текст подпрограммы и версий
aej3r_p.zip , aej3e_p.zip

Тексты тестовых примеров
taej3r_p.zip , taej3e_p.zip

Подпрограмма: AEJ3R (модуль AEJ3R_p)

Назначение

Вычисление собственных значений матрицы Якоби, принадлежащих данному интервалу, их номеров и соответствующих собственных векторов.

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

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

                   F = V ^T*A*V

(V - диагональная матрица), для которой методом бисекций вычисляются собственные значения, принадлежащие данному интервалу, и их номера. Соответственные собственные векторы матрицы F вычисляются методом обратных итераций. Связь между векторами матриц А и F осуществляется согласно соотношению: Х = V*Y.

Уилкинсон,Райнш, Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра, М., "Машиностроение", 1976.

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

procedure AEJ3R(N :Integer; var A :Array of Real;
                var V :Array of Real; var EV :Array of Real;
                var RLB :Real; var RUB :Real; MM :Integer;
                var M :Integer; var IRAB1 :Array of Integer;
                var RAB1 :Array of Real; var IERR :Integer);

   Параметры

N -	порядок исходной матрицы (тип: целый);
A -	вещественный двумерный массив размера N на 3, содержащий заданную в компактной форме матрицу Якоби;
V -	вещественный двумерный массив размерности N на ММ, содержащий в первых М столбцах вычисленные ортонормированные собственные векторы, соответствующие М собственным значениям;
EV -	вещественный вектор длины ММ, содержащий вычисленные в возрастающем порядке собственные значения, принадлежащие данному интервалу;
RLB - RUB	заданные нижняя и верхняя границы интервала собственных значений (тип: вещественный); если RLВ > RUВ, то собственные значения не вычисляются;
MM -	оценка сверху числа собственных значений исходной матрицы, принадлежащих данному интервалу, ММ ≤ N (тип: целый); если вычисленное число М собственных значений, принадлежащих данному интервалу, больше, чем MM, то собственные значения и соответствующие собственные векторы не вычисляются;
M -	целая переменная, в которой запоминается вычисленное число собственных значений, принадлежащих данному интервалу;
IRAB -	целый вектор длины MM, содержащий индексы расположенных в возрастающем порядке М собственных значений;
RAB -	вещественный вектор длины 8 на N, используемый как рабочий;
IERR -	целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом значение IЕRR

	- больше N, если не все попарные произведения соответствующих элементов побочных диагоналей неотрицательны; при этом исходная матрица не симметризируется;
	- больше 2*N, если имеется равное нулю произведение соответствующих элементов побочных диагоналей, причем сомножители равны нулю неодновременно. В этом случае нет симметризации, необходимой для правильного вычисления собственных векторов;
	- равно 3*N+1, если значение ММ меньше числа М вычисленных на данном интервале собственных значений, при этом собственные значения и соответствующие собственные векторы не вычисляются;
	- равно - k, если для вычисления собственного вектора с индексом k потребовалось более 5 итераций; при этом компоненты этого вектора полагаются равными нулю. Если таких собственных векторов несколько, то значение IЕRR полагается равным индексу последнего из них.

   Версии

AEJ3E -

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

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

UTAE10 -

подпрограмма выдачи диагностических сообщений при работе подпрограмм АЕJ3R и АЕJ3E.

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

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

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

Unit TAEJ3R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AEJ3R_p;

function TAEJ3R: String;

implementation

function TAEJ3R: String;
var
J,I,N,MM,_i,M,IERR :Integer;
RLB,RUB :Real;
V :Array [0..14] of Real;
EV :Array [0..2] of Real;
RAB1 :Array [0..39] of Real;
IRAB1 :Array [0..2] of Integer;
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
for _i:=0 to 14 do
 V[_i] := 0.0; { рабочий массив }
Result := '';  { результат функции }
Result := Result + #$0D#$0A;
Result := Result + Format('%s',
 [#$0D#$0A + ' ВЫЧИСЛЕНИЕ НЕКОТОРЫХ СОБСТВЕННЫХ ЗНАЧЕНИЙ И СОБСТВЕННЫХ BEKTOPOB' + #$0D#$0A +
 ' МАТРИЦЫ Я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;
N := 5;
ММ := 3;
RLB := -2.0;
RUB := 0.0;
AEJ3R(N,A,V,EV,RLB,RUB,MM,M,IRAB1,RAB1,IERR);
Result := Result + Format('%s',[#$0D#$0A + ' PEЗУЛЬTAT']) + #$0D#$0A; 
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[#$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;
Result := Result + Format('%s',[#$0D#$0A + ' EV' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[EV[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[#$0D#$0A + ' V' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  for J:=1 to 3 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',[#$0D#$0A + ' IERR' + #$0D#$0A]);
Result := Result + Format('%8d ',[IERR]) + #$0D#$0A;
UtRes('TAEJ3R',Result);  { вывод результатов в файл TAEJ3R.res }
exit;
end;

end.

Результаты:

Собственные значения в интервале (-2., 0.) :

      EV  =  ( -1.732050807,  -0.999999999 ) , 

Собственные векторы, соответствующие вычисленным в интервале
(-2., 0.) собственным значениям:

                   |  0.288675135 |                  | -0.499999999  |
                   | -0.499999999 |                  |  0.500000000  |
      V₁    =  |  0.577350269 |      V₂    =  | -3.18*10¹²       | , 
                   | -0.500000000 |                  | -0.499999999  |
                   |  0.288675135 |                  | -0.500000000  |

      IERR  =  0