Текст подпрограммы и версий
agh6r_p.zip  agh9r_p.zip  		Тексты тестовых примеров
tagh6r_p.zip  tagh9r_p.zip 

Подпрограмма:  AGH6R (модуль AGH6R_p) (версия AGH9R)

Назначение

Вычисление собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных вектоpов в обобщенной проблеме ABX = λx (или BAX = λx) для вещественных симметрических матриц A и B.

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

Подпрограмма AGH6R реализует алгоритм вычисления собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных вектоpов уравнения вида АВx = λx, где A и B - вещественные симметрические матрицы, и матрица B положительно определена.

При помощи разложения Холецкого для матрицы B:  В = LLT уравнение ABx = λx приводится к стандартному виду Qy = λy,  где Q = LTAL,  y = LTx. Стандартная задача решается путем приведения симметрической матрицы Q к трехдиагональному виду Q1 и вычисления собственных значений, принадлежащих заданному интервалу, методом бисекций и соответствующих собственных вектоpов, используя метод обратных итераций.

Восстановление собственных вектоpов  x  исходного уравнения АВx = λx из соответствующих вектоpов  z  уравнения Q1z = λz осуществляется согласно соотношению: x = L - TPz, где P - матрица преобразования Q1 = PTQP, L - преобpазование Холецкого, причем xTВx = E,  где E - единичная матрица.

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

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

procedure AGH6R(N :Integer; MM :Integer; var M :Integer;
                var RLB :Real; var RUB :Real; var A :Array of Real;
                var B :Array of Real; var EV :Array of Real;
                var V :Array of Real; var IRAB :Array of Integer;
                var RAB :Array of Real; var IERR :Integer);


   Параметры
N - порядок исходных матриц (тип: целый);
MM - оценка свеpху числа собственных значений уpавнения ABx = λx, принадлежащих заданному интервалу MM ≤ N (тип: целый); если фактическое число собственных значений M, принадлежащих заданному интервалу, больше, чем MM, то собственные значения и соответствующие собственные векторы не вычисляются;
M - целая переменная, в которой запоминается вычисленное число собственных значений, принадлежащих заданному интервалу;
   RLB -
   RUB   		заданные нижняя и верхняя границы интервала собственных значений (тип: вещественный); если RLB > RUB, то собственные значения не вычисляются;
A, B - вещественные двумерные массивы размера N на N, содержащие исходные симметрические матрицы A и B соответственно;
EV - вещественный вектоp длины MM, содержаший вычисленные в возрастающем порядке собственные значения;
V - вещественный двумерный массив размера N на MM, содержащий в первых M столбцах вычисленные собственные векторы, нормированные так, что VTBV = E;
IRAB - целый вектоp длины MM, содержащий индексы расположенных в возрастающем порядке M собственных значений;
RAB - вещественный вектоp длины 9 на N, используемый как рабочий;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом значение IERR:
 

- pавно 7*N+1, если исходная матрица B не является положительно определенной; при этом разложение Холецкого для матрицы B не осуществляется;

- pавно 3*N+1, если значение MM меньше истинного числа вычисленных собственных значений M на интервале; при этом собственные значения и соответствующие собственные векторы не вычисляются.

- pавно - K, если для вычисления собственного вектоpа с индексом K потребовалось более 5 итераций; при этом компоненты этого вектоpа полагаются равными нулю. Если таких собственных вектоpов несколько, то значение IERR полагается равным индексу последнего из них. 
   Версии
AGH9R - вычисление собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных вектоpов в обобщенной проблеме BAx = λx для вещественных симметричных матриц A и B. 
   Вызываемые подпрограммы
UTAG10 - подпрограмма выдачи диагностических сообщений при работе подпрограмм AGH6R и AGH9R. 
   Замечания по использованию
 

Исходные матрицы A и B можно задавать лишь верхними треугольными половинами.

Подпрограмма AGH6R сохраняет строгий верхний треугольник массива A и полный верхний треугольник массива B, остальные элементы массивов A и B используются как pабочие.

Программа AGH9R нормирует собственные векторы  x  таким образом, что xTB - 1x = E.

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

Unit TAGH6R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AGH6R_p;

function TAGH6R: String;

implementation

function TAGH6R: String;
var
J,I,N,MM,M,IERR :Integer;
RLB,RUB :Real;
V :Array [0..24] of Real;
EV :Array [0..4] of Real;
RАВ :Array [0..44] of Real;
IRАВ :Array [0..4] of Integer;
const
A :Array [0..24] of Real = ( 10.0,0.0,0.0,0.0,0.0,2.0,12.0,0.0,0.0,0.0,3.0,
1.0,11.0,0.0,0.0,1.0,2.0,1.0,9.0,0.0,1.0,1.0,-1.0,
1.0,15.0 );
B :Array [0..24] of Real = ( 12.0,0.0,0.0,0.0,0.0,1.0,14.0,0.0,0.0,0.0,-1.0,
1.0,16.0,0.0,0.0,2.0,-1.0,-1.0,12.0,0.0,1.0,1.0,
1.0,-1.0,11.0 );
begin
Result := '';  { результат функции }
Result := Result + Format('%s',
 [' ВЫЧИСЛЕНИЕ НЕКОТОРЫХ СОБСТВЕННЫХ ЗНАЧЕНИЙ И COOTBETCTBYЮЩИX' + #$0D#$0A +
 ' СОБСТВЕННЫХ ВЕКТОРОВ YРАВНЕНИЯ BИДA:ABX=LAMBDA*X,' + #$0D#$0A +
 ' ГДЕ A И B - ВЕЩЕСТВЕННЫЕ СИММЕТРИЧЕСКИЕ MATPИЦЫ']) + #$0D#$0A; 
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' A' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  for J:=1 to 5 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' B' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  for J:=1 to 5 do
   begin
    Result := Result + Format('%20.16f ',[B[(I-1)+(J-1)*5]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
N := 5;
ММ := 5;
RLB := 70.0;
RUB := 300.0;
AGH6R(N,MM,M,RLB,RUB,A,B,EV,V,IRAB,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 5 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*5]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' B' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  for J:=1 to 5 do
   begin
    Result := Result + Format('%20.16f ',[B[(I-1)+(J-1)*5]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' EV' + #$0D#$0A]);
for I:=1 to 5 do
 begin
  Result := Result + Format('%20.16f ',[EV[I-1]]) + #$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('TAGH6R',Result);  { вывод результатов в файл TAGH6R.res }
exit;
end;

end.

Результаты:

Собственные значения в интервале (70., 300.)

                 |  77.697191195 |
                 | 112.15419325  |
      EV  =  | 134.68646332  |
                 | 167.48487891  |
                 | 242.97727332  |

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

                 |  0.2349114135 |                      |  0.1288556917 |
                 | -0.0410915167 |                     | -0.1193865988 |
      V1  =  | -0.0383075946 |    ,      V2  =  | -0.0282771880 |    , 
                 | -0.2059003675 |                     |  0.1923580004 |
                 | -0.0734707966 |                     | -0.0097623271 |
 
                 |  0.0042355205 |                     |  0.0183136812 |
                 | -0.1812063856 |                    | -0.0266749519 |
      V3  =  |  0.1210383986 |    ,      V4  =  |  0.1834456078 |    , 
                 | -0.0609182758 |                     |  0.0051904406 |
                 |  0.1690213925 |                     | -0.2218442867 |

                 | -0.1249195279 |
                 | -0.1535463561 |
      V5  =  | -0.1145245145 |    ,      IERR  =  0
                 | -0.0657938487 |
                 | -0.1010161054 |