Подпрограмма:  AES0R (модуль AES0R_p)

Назначение

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

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

Метод Ланцоша для заданной симметричной матрицы А порядка N и начального вектора v ∈ Е^N ,  v ≠ 0 ,  строит последовательность векторов Ланцоша  v₁, v₂, ..., v_j и последовательность трехдиагональных матриц  Т₁, Т₂, ..., Т_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-1v1 }  ,   j = 1, ..., N  , 

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

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

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

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

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

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

Выход из подпрограммы происходит, когда все требуемые собственные пары будут вычислены с точностью ЕРS по невязке, где ЕРS задаваемая величина, или по заданному максимально допустимому числу NI крупных шагов.

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

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

procedure AES0R(N :Integer; KL :Integer; KR :Integer;
                var KL0 :Integer; var KR0 :Integer; M :Integer;
                var NI :Integer; var VJ :Array of Real; FF :Proc_F4A;
                EPS :Real; L :Integer; var R :Array of Real;
                var EV :Array of Real; var Y :Array of Real;
                var ER :Array of Real; var IP :Integer);

   Параметры 

   Версии :   нет 

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

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

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

FAES0R - подпрограмма пользователя, вычисляющая W:  = (A - S*I)V - W
       для A = diag (0., 1.E-5, 2.E-5, 3.E-5, 0.1, 1.)

Unit TAES0R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FAES0R_p, AES0R_p;

function TAES0R: String;

implementation

function TAES0R: String;
var
KL0,KR0,IP,NI,_i,N :Integer;
Y :Array [0..11] of Real;
EV :Array [0..1] of Real;
ER :Array [0..1] of Real;
RАВ :Array [0..210] of Real;
const
VJ :Array [0..5] of Real = ( 1.0,1.0,1.0,1.0,1.0,1.0 );
begin
Result := '';  { результат функции }
KL0 := 0;
KR0 := 0;
IP := 0;
NI := 1;
AES0R(6,0,2,KL0,KR0,5,NI,VJ,FAES0R,1.E-6,211,RAB,EV,Y,ER,IP);
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;
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;
Result := Result + Format(' %5d %5d ',[IP,NI]) + #$0D#$0A;
N := 6;
Result := Result + #$0D#$0A;
for _i:=0 to 11 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TAES0R',Result);  { вывод результатов в файл TAES0R.res }
exit;
end;

end.

Результаты:

       EV = (1., 0.1) ,   ER = (2.2E-11, 2.2E-9) ,   IP = 0 ,   NI = 1 ,
  
       Y1  =  (-1.1E-11,  1.1E-11,  1.1E-11,  -1.2E-11,  -2.0E-13,  1.)T
       Y2  =  (-1.1E-8,  1.1E-8,  1.1E-8,  -1.1E-8,  -1.,  1.2E-10) .