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

Назначение

Решение интегрального уравнения Фредгольма I рода методом поточной невязки с дескриптивной регуляризацией.

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

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

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

с гладким ядром и правой частью  f (x), заданной с погрешностью  δ (x), при априорно известной информации об ограниченности первой призводной исходного решения  u (y) и об участках знакопостоянства его кривизны сводится методом поточной невязки с дескриптивной регуляризацией [1,2] к отысканию общей точки  u (y) множеств  U_δ , U_ν , U_μ :

(2)     Uδ = { u:  | Au - f(x) | ≤ δ (x) }  ,
(3)     Uν = { u:  v1(y) ≤ u'(y) ≤ v2(y)   ∀y ∈ [a, b] }  ,
(4)     Uμ = { u:  μ(y) u''(y) ≥ 0   ∀y ∈ [a, b] ,  μ(y) = sign u''(y) } 

т.е. к выделению из множества  U_δ формальных решений уравнения (1) приближенного решения  u (y), удовлетворяющего условиям (3), (4). v₁ (y) и  v₂ (y) - заданные на отрезке [a, b] функции.

Дискретный аналог задачи (2) - (4) состоит в определении вектоpа  u = (u₁, ..., u_N)∈E_N, удовлетворяющего системе линейных неpавенств:

(5)     | (K Du)i - fi | ≤ δi ,      i = 1, 2, ..., M  ,
(6)     v1 j ≤ (uj+1 - uj) / hj ≤ v2 j ,      j = 1, 2, ..., N-1  ,
(7)     μj ( (uj+1 - uj) / hj - (uj - uj-1) / hj-1 ) ≥ 0 ,      j = 2, 3, ..., N-1  , 

где  u = (u₁, ..., u_N)∈E_N,  K - матрица порядка  M * N с элементами  K_i j = k (x_i, y_j),  {x_i} - сетка узлов на отрезке [c, d]:  c = x₁ < ... < x_i < ... < x_M = d,  {y_j} - сетка узлов на отрезке [a, b]:  a = y₁ < ... < y_j < ... < y_N = b,  y_j+1 - y_j = h_j,  D - диагональная матрица коэффициентов  d_j квадратурной формулы,  j = 1, 2, ..., N,  f_i = f (x_i) и  δ_i = δ (x_i) - значения правой части и погрешности в узлах  x_i,  i = 1, 2, ..., M. Через  v_1 j, v_2 j и  μ_j обозначены значения функций  v₁ (y), v₂ (y) и  μ (y) в узлах  y_j.

Для приближенного решения системы неpавенств (5) - (7) применяется численный алгоритм, разработанный в [3]. B случае пустоты множества (5) при заданном уpовне погрешности в подпрограмме предусмотрено нахождение значений  δ₁ = δ_i + p δ_i (p > 0 - заданный параметр), при которых множество (5) непусто.

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

    SUBROUTINE  EF04R (M, N, A, F, DELTA, P, NG, V1, V2,
                                            MU, H, UP, U, B) 

   Параметры 

   Версии:  нет 

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

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

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

   Рассматривается интегральное уравнение
          1
          ∫ k(x, y) - u(y) dy = f(x)
        -1
   с ядром        k(x, y) = (1 + (x - y)2)-1 ,

   правой частью
        f(x) = (2 - x2) ( arctg(1 - x) + arctg(1 + x) ) - 2 -
                      - x ln( (1 + (1 - x)2) / (1 + (1 + x)2) )

   и точным решением     u(y) = 1 - y2 . 

Ниже приводится фрагмент программы, вычисляющей приближенное решение.

Для аппроксимации интеграла применяется квадратурная формула трапеций с использованием равномерной сетки по  x на отрезке [- 2, 2] с шагом  h_i = 0.1 (число узлов  M = 41) и равномерной сетки по  y на отрезке [- 1, 1] с шагом  h_j = 0.05 (число узлов  N = 41). Решение ищется в классе функций  u (y) со следующими свойствами:

         u'(y) ≥ 0   при -1 ≤ y ≤ 0 ,
         u'(y) ≤ 0   при  0 < y ≤ 1 ,
         u''(y) ≤ 0  при -1 ≤ y ≤ 1.

       REAL  A(41, 41), F(41), DELTA(41), V1(41), V2(41), MU(41), H(41)
       REAL  UP(41), U(41), HX(41), B(619)
       DATA  DELTA /41*0.001/, V1 /20*0., 20*100., 0./, H /40*0.05, 0./
       DATA  V2 /20*100., 21*0./, MU /39*1, 2*0/, UP /41*0./, P /1./
       DATA  M /41/, N /41/, NG /3/, HX /40*0.1, 0./
       DO 1  J = 1, N
       V1(J) = -V1(J)
    1 MU(J) = -MU(J)
       B(1) = -2.
       B(M + 1) = -1.
       DO 2  J = 2, N
       K = M + J
    2 B(K) = B(K-1) + H(J-1)
       DO 3  I = 2, M
    3 B(I) = B(I-1) + HX(I-1)
       DO 4  I = 1, M
       DO 4  J = 1, N
       K = M + J
    4 A(I, J) = 1./( 1. + ( B(K) - B(I) )**2 )
       DO 5  I = 1, M
       W = 1. + B(I)
       W1 = 1. - B(I)
       F(I) = (1. + W*W1)*( ATAN(W) + ATAN(W1) ) - 2. -
      *          B(I)*ALOG( (1. + W1*W1)/(1. + W*W) )
    5 CONTINUE
       CALL  EF04R (M, N, A, F, DELTA, P, NG, V1, V2, MU, H, UP, U, B)

Результаты:

       при  J = 1 + 5K ,    K = 0, 1, ..., 8

  U(J) = ( 0.056, 0.415, 0.716, 0.950, 1.038, 0.942, 0.707, 0.416, 0.063 )