Подпрограмма:  AES4R (модуль AES4R_p)

Назначение

Вычисление нескольких крайних собственных значений и соответствующих собственных векторов вещественной симметричной разреженной или имеющей регулярную структуру матрицы методом Ланцоша с выборочной ортогонализацией.

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

Метод Ланцоша для заданной симметричной матрицы А порядка N и начального вектора  v ∈ E^N,  v ≠ 0, строит последовательность векторов  v₁, v₂,..., v_j и последовательность трехдиагональных матриц T₁, T₂,..., T_j по следующим формулам

(1)       v1 = v / (vTv)1/2  ,     β1 = 0  , 
            ui = Avi - βivi -1  ,     αi = uiTvi  , 

  T1 =  α1 ,                            i = 1 ,

           |               |   0  |
           |               |   .   |
  Ti =   |    Ti-1     |  0   |  ,          i > 1 ,
           | ______ |  βi  |
           | 0 ... 0  βi  αi  |

   wi +1 = ui - αivi ,    βi +1 = (wTi +1wi +1)1/2 ,

   vi +1 = wi +1 / βi +1  ,    i = 1,..., j 

Метод Ланцоша может быть рассмотрен как реализация для спектральной задачи Az = λz метода проекций Ритца - Галеркина на последовательность подпространств Крылова

     K( j )  =  span { v1,  Av1,... , Aj -1 v1 }   ,     j = 1,... , N   , 

при этом числа Ритца (собственные значения сужения оператора P_jA на K ^( j ), где P_j - ортогональный проект на K ^( j ) ) совпадают с собственными значениями θ_i ^( j ), i = 1,..., j, матрицы T_j, получаемой на  j - ом шаге метода Ланцоша, а векторы Ритца равны V_js_i ^( j ), где  s_i ^( j ) - соответствующие θ_i ^( j ) собственные векторы T_j, а V_j - матрица, столбцами которой являются векторы Ланцоша  v₁, v₂,..., v_j.

Полезной особенностью метода Ланцоша является то, что часто уже при  j ~ 2√N крайние пары Ритца ( θ_i ^( j ), V_j s_j ^( j ) ) оказываются хорошими приближениями к искомым крайним собственным парам матрицы А.

Метод Ланцоша использует исходную матрицу только в операциях умножения матрицы на вектор, поэтому информацию о матрице удобно задавать в виде подпрограммы умножения матрицы на вектор, в которой пользователь может максимально учесть специфику матрицы.

В подпрограмме AES4R используется метод Ланцоша с выборочной ортогонализацией - устойчивая к влиянию ошибок округления модификация простого метода Ланцоша, отличающаяся от него тем, что на некоторых шагах процесса (1) проводится ортогонализация  w_j + 1 относительно некоторых векторов Ритца.

В подпрограмме AES4R предполагается, что || A || ₂ ~ 1, где || A || ₂ - спектральная норма матрицы. Гарантированной точностью по невязке для вычисляемых собственных векторов является величина  k || A ||₂√ε, где  ε - относительная машинная точность,  k ~ 1, но часто возможно получение лучшей точности.

Выход из подпрограммы происходит, когда все требуемые собственные векторы будут вычислены с заданной точностью по невязке или если продолжение процесса (1) окажется невозможным или бесполезным.

Б.Парлетт. Симметричная проблема собственных значений. Численные методы. М.: Мир, 1983.

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

procedure AES4R(N :Integer; KL :Integer; KR :Integer;
                FF :Proc_F4A; EPS :Real; var M :Integer;
                var VJ :Array of Real; K :Integer; L :Integer;
                var R :Array of Real; var R1 :Array of Real;
                var R2 :Array of Real; var EV :Array of Real;
                var ER :Array of Real; var V :Array of Real;
                var Y :Array of Real; var IP :Integer);

   Параметры 

   Версии 

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

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

Unit TAES4R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FAES4R_p, AES4R_p;

function TAES4R: String;

implementation

function TAES4R: String;
var
M,IP,_i,J,I :Integer;
EV :Array [0..1] of Real;
ER :Array [0..1] of Real;
R :Array [0..96] of Real;
R1 :Array [0..5] of Real;
R2 :Array [0..5] of Real;
V :Array [0..29] of Real;
Y :Array [0..11] of Real;
const
VJ :Array [0..5] of Real = ( 1.0,1.0,1.0,1.0,1.0,1.0 );
label
_55;
begin
for _i:=0 to 96 do
 R[_I] := 0.0;  { РАБОЧИЙ МАССИВ }
ResULT := '';  { РЕЗУЛЬТАт функции }
M := 5;
IP := 0;
AES4R(6,0,2,FAES4R,1.E-3,M,VJ,4,97,R,R1,R2,EV,ER,V,Y,IP);
Result := Result + Format('%s',['  IP=']);
Result := Result + Format('%5d ',[IP]);
Result := Result + Format('%s',['   M=']);
Result := Result + Format('%5d ',[M]) + #$0D#$0A;
if ( IP < 0 )
 then goto _55;
Result := Result + Format('%s',['  EV=' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[EV[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  Y=C.B.' + #$0D#$0A]);
for I:=1 to 6 do
 BЕgin
  for J:=1 to 2 do
   begin
    Result := Result + Format('  %20.16f ',
 [Y[(I-1)+(J-1)*6]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  ER=' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[ER[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
_55:
UtRes('TAES4R',Result);  { вывод результатов в файл TAES4R.res }
exit;
end;

end.

Результаты:

         IP  =  0   ,     M  =  3
        
         Собственные значения:

         EV(1)  =  1.   ,     EV(2)  =  0.1
         
         Собственные векторы:  

                   |   5.E - 7           2.E - 4  |
                   | - 6.E - 7           5.E - 5  |
                   | - 2.E - 6         - 6.E - 5  |
         Y  =   | - 3.E - 6         - 2.E - 4  |
                   |   1.E - 6         - 1.          |
                   |   1.                   3.E - 6  |
                     
         Невязки:     ER(1)  =  4.E - 6   ,     ER(2)  =  2.E - 5