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

Назначение

Решение интегрального уравнения Вольтерра со стационарным ядром методом рягуляризации первого порядка без выбора параметра регуляризации.

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

Данная подпрограмма реализует алгоритм, изложенный в [1].

Численное решение интегрального уравнения Вольтерра 1 рода вида

                T
(1)           ∫   k (t - x)  u (x)  dx  =  f (t)
              t
                t
           (   ∫   k (t - x)  u (x)  dx  =  f (t) )
              0

         t ∈ [ 0, T ] 

на равномерной сетке  t_i = i H,  i = 1, 2, ..., n,  H = Т/n, сводится к задаче минимизации

(2)           || K u - f || E - min u , 

где K - верхняя (нижняя) треугольная теплицева матрица размера  n на n, полученная при дискретизации (1):   K_i j  =  H * k (H ( i - j )).

Используя метод регуляризации А.Н.Тихонова [2], приходим к следующей задаче

(3)           || K u - f || E2  +  α || D u || E2  -  min u , 

где α > 0 - заданный параметр регуляризации, D - дискретный аналог оператора дифференцирования   d /dx.

Задача (3) эквивалентна задаче

                ||     |    K      |           |  f   |     ||  
(4)           ||     |             |  u   -   |      |     ||      -   min u ,                 
                ||     | α1/2 D |           |  0  |     || E 

для решения которой при помощи метода вращений [3] определим ортогональную матрицу Q, такую, что

              |     K    |            |  Kα  |                             |   f    |            |   f1    |
    QT     | α1/2 D |     =     |   0    |  ,               QT     |   0   |     =     |   f2    | 

где K_α - верхняя (нижняя) треугольная матрица. Тогда решение задачи (4) сводится к решению системы линейных алгебраических уравнений с треугольной матрицей:

(5)           Kα u  =  f1 , 

котоpое осуществляется методом обратной подстановки.

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

    SUBROUTINE  EV01R (A, F, N, L, ALFA, H, B, C, U, IERR) 

   Параметры 

   Версии:  нет 

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

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

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

 Решается интегральное уравнение
                      t
                     ∫   K (t - s)  f (s)  ds  =  g (t)   ,      0 ≤ t ≤ 1
                   0
с ядром
                              ∞
      K(t)  =  0.5 π   ∑    (- 1)n (2nH )  exp (-(2n + 1)2 π2 t /8) .
                             n=0 

Точное значение функции  f (s) задано таблицей. Используя формулу трапеции, вычисляем приближенное значение правой части  g (t), затем по  g (t) восстанавливаем значение функции  f (s) = u с помощью программы EV01R.

Параметр регуляризации  α = 10 ^- 8

       REAL  R(50), F(50), G(100), U(50), WORK(1325), K(50)
       DATA  F /0.025, 0.215, 0.405, 0.595, 0.785, 0.975, 1., 0.92,
     *                0.71, 0.5, 0.29, 0.08, 0., 0.08, 0.29, 0.5, 0.71, 0.92,
     *                1., 0.975, 0.785, 0.595, 0.405, 0.215, 0.025, 25*0./
       N = 50
       ALFA = 1.E - 8
       H = 2.E - 2
       PI = 3.14159265
       N1 = N - 1
       DO 11  I = 1, N
       G(N + I) = 0
       S = EXP( - (PI*PI*I*H/8) )
       DO 12  J = 1, 15
       RAB = EXP( - (2*J + 1)**2*PI*PI*I*H/8 )
  12 S = S + (- 1)**J*(2*J + 1)*RAB
  11 R(I) = 0.5*PI*S
       K(1) = R(1)*H/2
       K(N) = R(N)*H/2
       DO 13  I = 2, N1
  13 K(I) = R(I)*H
       DO 14  I = 1, N
       S = 0
       DO 15  J = 1, I
  15 S = S + F(J)*K(I - J + 1)
  14 G(I) = S
       L = 1
       CALL  EV01R (R, G, N, L, ALFA, H, WORK(N + 1), WORK, U, IERR)

Результаты:

          ALFA  =  1.00000000-08

               |  4.18939426-02   9.47262731-01   2.69209987-01    4.94652786-01
      U  =  |  8.07735236-01   9.34217208-04   8.94488047-03   -2.22705150-03
               | -3.81489539-04   5.87342988-03