Подпрограмма:  AES4R

Назначение

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

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

Метод Ланцоша для заданной симметричной матрицы А порядка 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 относительно некоторых векторов Ритца.

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

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

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

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

    SUBROUTINE  AES4R (N, KL, KR, FF, EPS, M, VJ, K, L, R,
                                             R1, R2, EV, ER, V, Y, IP) 

   Параметры 

   Версии 

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

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

       SUBROUTINE  F0 (X, Z, N, S)
       DIMENSION  X(N), Z(N), RAB(6)
       DATA  RAB /0., 1.E - 5, 2.E - 5, 3.E - 5, 0.1, 1./
       DO  10  I = 1, N
       Z(I) = (RAB(I) - S) * X(I) - Z(I)
  10 CONTINUE
       RETURN
       END

       DIMENSION  VJ(6), EV(2), ER(2), R(97), R1(6), R2(6), V(6,5), Y(6,2)
       EXTERNAL  F0
       DATA  VJ /6*1./
       M = 5
       IP = 0
       CALL  AES4R (6, 0, 2, F0, 1.E - 3, M, VJ, 4, 97, 
      *                         R, R1, R2, EV, ER, V, Y, IP)

Результаты:

         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