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

Назначение

Приближенное решение интегрального уравнения Фредгольма I рода в классе ограниченных кусочно - монотонных функций.

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

Приближенное решение интегрального уравнения Фредгольма I рода

            b
           ∫ k(x, y) u(y) dy  =  f(x)  ,  c ≤ x ≤ d ,
          a 

при априорно известной информации о качественном поведении искомого решения  u = u (y) осуществляется методом квазирешений В.К.Иванова путем сведения к решению задачи квадратичного программирования с линейными ограничениями

(1)     || K D u - f ||p2  =   inf   || K D v - f ||p2 ,
                                      v∈G 

где  K - известная матрица порядка  M * N с элементами  k_i j , которая в зависимости от ядра  k (x, y) и способа его аппроксимации может быть матрицей общего вида, либо симметричной, циркулянтной, теплицевой и т.п.,

f = ( f₁, ..., f_M ) - известный  M -  мерный вектор ,

u = ( u₁, ..., u_N ) - искомый  N - мерный вектор приближенного решения,

D - диагональная матрица коэффициентов квадратурной формулы  d_j ,  j = 1, ..., N,

                        M 
       || Z ||p2  =  ∑  pi zi2 ,      pi > 0  -  весовые коэффициенты.
                       i=1 

Линейные ограничения  G представляют собой одно из множеств  G₁,  G₂ или  G₁ ∩ G₂, где

 G1 = { v = (v1, ..., vN) :  aj ≤ vj ≤ bj ,  j = 1, ..., N } ,
 G2 = { v = (v1, ..., vN) :  mj (vj + 1 - vj) ≥ 0 ,  j = 1, ..., N - 1 } , 

a = (a₁, ..., a_N),  b = (b₁, ..., b_N) - заданные  N - мерные векторы нижних и верхних ограничений в  G₁; параметры  m_j определяются условиями кусочной монотонности и принимают значения 1, 0 или -1.

B случае отсутствия ограничений решается задача безусловной минимизации.

Численное решение задачи (1) основано на методе проекции сопряженных градиентов, который позволяет проводить итерационный процесс в конечномерном пространстве  W_2, ρ¹ с нормой

                                    N
 || v ||   =  ( || v || 2Q  +  ∑   ρj (vj - vj - 1)2 )1/2 ,
                                   j=2
                     N
 || v || Q  =  (  ∑   q j vj2 )1/2 ,
                    j=1 

q _j > 0 и  ρ_j ≥ 0 - заданные числа. Если априори известно, что искомое решение достаточно гладкое, то задание параметров  ρ_j > 0, связанных с конкретной нормировкой  W_2, ρ¹ ( при  ρ_j≡0 эта ноpма  v  совпадает с  || v || _Q ), позволяет увеличить точность и улучшить качество приближенного решения.

Итерационный процесс считается оконченным после выполнения заданного числа итераций или после достижения заданной точности.

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

B самой подпрограмме EF01R матрица  K не используется, поэтому способ ее представления не регламентируется.

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

    SUBROUTINE  EF01R (M, N, A, F, NG, AU, BU, MU, P, Q,
                                            D, R, E, ISP, U, FU, KN, B, OP) 

   Параметры 

                N
      wi  =  ∑   ki j vj ,    i = 1, ..., M ,
               j=1 

                M
      wj  =  ∑   kj i vi ,    j = 1, ..., N .
               i=1 

   Версии:  нет 

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

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

                             xi+1    yj+1
      ki j  =  1/h1x       ∫       ∫     k (x, y)  dx dy ,   dj ≡ 1 .
                              xi      yj 

  1 CALL  EF01R (M, N, A, F, NG, AU, BU, MU, P, Q, D, R, E, ISP,
    *                        U, FU, KN, OP)
     IF (KN) 2, 1, 2
  2 CONTINUE 

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

Матрица  K имеет общий вид

       DIMENSION  A(3, 3), F(3), MU(3), P(3), Q(3), D(3), R(4), U(3), B(26), AU(3), BU(3)
       DATA  A /3., 1., 5., 1., -3., 2., 1., 2., 1./
       DATA  F(1), F(2), F(3) /4., 3., 4./,  MU(1), MU(2), MU(3) /2*1, 0/
       DATA  P(1), P(2), P(3) /3*1./,  Q(1), Q(2), Q(3) /3*1./
       DATA  D(1), D(2), D(3) /3*1./,  R(1), R(2), R(3), R(4) /4*0./
       DATA  U(1), U(2), U(3) /3*0./
       EXTERNAL  EF01R1
       M = 3
       N = 3
       NG = 2
       E = 1.E - 10
       ISP = 20
       CALL  EF01R (M, N, A, F, NG, AU, BU, MU, P, Q, D, R, E, ISP,
      *                        U, FU, KN, B, EF01R1)

Результат:    U  =  (-1., 2., 5.)