Текст подпрограммы и версий ( Фортран )
ev01r.zip 		Тексты тестовых примеров ( Фортран )
tev01r.zip
Текст подпрограммы и версий ( Си )
ev01r_c.zip 		Тексты тестовых примеров ( Си )
tev01r_c.zip
Текст подпрограммы и версий ( Паскаль )
ev01r_p.zip 		Тексты тестовых примеров ( Паскаль )
tev01r_p.zip

Подпрограмма:  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 ]

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

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

где K - верхняя (нижняя) треугольная теплицева матрица размера  n на n, полученная при дискретизации (1):   Ki 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ое осуществляется методом обратной подстановки.

1.  L.Elden. An efficient algorithm for the regularization of ill - conditioned least squares problems with triangular toeplitz matrix. - REPORT  LITH - MAT - R - 1981 - 25, Sweden, 1981.
2.  А.Н.Тихонов, В.Я.Арсенин. Методы решения некорректных задач. - М., Наука, 1979.
3.  В.В.Воеводин. Численные методы алгебры. - M., Hаука, 1966.

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

    SUBROUTINE  EV01R (A, F, N, L, ALFA, H, B, C, U, IERR) 
   Параметры
A - вещественный вектоp длины N, содержащий первую стpоку (первый столбец) теплицевой матрицы K;
F - вещественный вектоp длины 2N, содержащий в первых N ячейках компоненты заданного вектоpа правой части; последующие N ячеек используются как рабочие;
N - заданная размерность вектоpов F и U и число узлов сетки разбиения (тип: целый);
L - заданный параметр - указатель типа матрицы K:
 

L = 0, если матрица верхняя треугольная,

L = 1, если матрица нижняя треугольная (тип: целый);

ALFA - заданный параметр регуляризации (тип: вещественный);
H - заданный шаг разбиения сетки (тип: вещественный);
B - вещественный вектоp длины N (N + 1)/2, используемый как рабочий;
C - вещественный вектоp длины N, используемый как рабочий;
U - вещественный вектоp длины N, содержащий на выходе из подпрограммы вычисленное решение интегрального уравнения в узлах заданной сетки;
IERR - целая переменная, сообщающая об ошибках, обнаруженных в ходе вычислений;
IERR= 0 - если ошибок не было обнаружено;
IERR=65 - в случае, если при заданном значении ALFA матрица системы вырождается. 
   Версии:  нет 
   Вызываемые подпрограммы
UTEV10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы EV01R. 
   Замечания по использованию:  нет

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

 Решается интегральное уравнение
                      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