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

Назначение

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

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

Решение интегрального уравнения Фредгольма I рода:

                     b
(1)    Au  ≡   ∫  K(y, x)  u(x) dx  =  f(y)  ,    y ∈ [c, d] ,
                    a 

методом регуляризации А.Н.Тихонова [1] сводится к минимизации параметрического функционала:

          d           b
(2)     ∫ P(y) [  ∫  K(y, x) u(x) dx - f(y) ]2 dy  + α Ω(u) ,
        c            a 

где  P (y) > 0 - весовая функция,  α > 0 - параметр регуляризации,  Ω (u) - стабилизирующий функционал, определяемый параметрами  P1 > 0,  P2 ≥ 0,  P3 ≥ 0 :

                     b
       Ω(u)  =  ∫  ( P1 u2 + P2 (du/dx)2 + P3(d2u/dx2)2 ) dx
                    a 

Для дискретизации первого слагаемого в (2) используются квадратурные формулы трапеций по обеим переменным, а при аппроксимации второго слагаемого - формулы численного дифференцирования второго порядка для аппроксимации первой и второй производных функции  u (x).

   Тогда дискретный аналог (2) имеет вид:
         M                N
(3)     ∑  Pi σi''  [  ∑  σj' Ki j uj - fi ]2  +  α u TC u .
         i=1               j=1 

Здесь  P_i = P (y_i),  u_j = u (x_j),  f_i= f (y_i),  K_i j = K (y_i, x_j);  a =x₁ < x₂ < ... < x_N = b;  c = y₁ < y₂ < ... < y_M = d - сетки узлов по  x и  y,  σ_j' и  σ_i'' - коэффициенты квадратурных формул по  x и  y соответственно; C - матрица квадратичной формы, аппроксимирующей стабилизирующий функционал  Ω (u);  u = (u₁, u₂, ..., u_N) - искомое решение.

Пусть f = ( f₁ √σ₁'',  f₂ √σ₂'', ..., f_M √σ_M'' ) и  A = (a_i j),  a_i j = K_i j P_i σ_j' √σ_i'',  1 ≤ i ≤ M,  1 ≤ j≤ N.

Тогда задача минимизации (3) по  u эквивалентна решению системы линейных алгебраических уравнений:

(4)     ( AT A + α C ) u  =  AT f ,
   где T означает операцию транспонирования. 

При решении этой системы используется сведение задачи к более простой, канонической задаче с единичной матрицей  C и двухдиагональной матрицей  A согласно методу В.В.Воеводина [2]. Подпрограмма вычисляет значения пяти критериальных функций, которые могут быть использованы для выбора параметра регуляризации. A именно, в подпрограмме вычисляются приближенные:

1)  Невязка
                      d           b
     ρ(α)  =  (  ∫ P(y) [  ∫  K(y, x) uα(x) dx - f(y) ]2 dy )1/2 ;
                     c           a
2)  Функция  "чувствительности"
     τ(α)  =  α Ω1/2 (duα / dα) ;

3)  Hоpма решения
      γ(α)  =  Ω1/2 (uα ) ;

4)  Значение регуляризирующего функционала
      φ(α)  =  ρ2(α) + αγ2(α) ;

5)  Функция отношения
        y(α)  =  ρ2(α) / ρ(α) ,
где
                        d           b
     ρ2(α)  =  {  ∫ P(y) [  ∫  K(y, x) ( α d uα(x)/dα - uα(x) ) dx - f(y) ]2 dy }1/2 .
                       c           a 

Здесь  u_α = u_α (x) - функция, минимизирующая (2).

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

Описываемая подпрограмма реализует эти возможности и позволяет:

При работе с задачей, приведенной к каноническому виду, подпрограмма позволяет:

1. Тихонов A.H., O регуляризации некоppектно поставленных задач, ДАН CCCP, 1963, т.153, N 1.

2. Воеводин B.B., O методе регуляризации, ЖВМ и МФ, 1969, т.9, N 3.

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

    SUBROUTINE  EF13R (SX, SY, N, M, A, U, F, P, P1, P2, P3,
                                           Q, H, W, L) 

   Параметры 

   Версии:  нет 

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

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

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

Рассматривается решение интегрального уравнения

         1
         ∫ k(y, x) u(x) dx = f(y)
        0 

с ядром  k(x, y) = y / (1 + y² x²) и правой частью  f (y) = arctg y (точное решение  u ≡ 1) или  f (y) = ln (1 + y²) / 2y (точное решение  u ≡ x).

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие знначения критериальных функций при фиксированном параметре регуляризации, с весовой функцией  p ≡ 1 и использованием равномерных по  x и  y на отрезке [0,1] сеток из 11 точек.

       REAL  SX(11), SY(11), P(11), U(11), F(11), A(11, 11), H(50), W(11, 8)
       DATA  P /11*1./,  N /11/,  M /11/,  Q /1.E - 6/
       DATA  SX /0., 0.1, .2, .3, .4, .5, .6, .7, .8, .9, 1./
       DATA  SY /0., .1, .2, .3, .4, .5, .6, .7, .8, .9, 1./
       DO 1  I = 1, M
       F(I) = ATAN(SX(I))
       DO 1  J = 1, N
    1 A(I, J) = SY(I) / ( 1. + SX(J) * SX(J) * SY(I) * SY(I) )
       CALL  EF13R (SX, SY, N, M, A, U, F, P, 1., 0., 0., Q, H, W, 1)
       CALL  EF13R (SX, SY, N, M, A, U, F, P, 1., 0., 0., Q, H, W, 2)

Результаты:

       U  =  ( 0.99, 0.99, 0.99, 1.0, 1, 0, 1.0, 1.0, 1.0, 0.99, 0.98, 0.96 )
       H  =  ( 5.6*10-6,  5.7*10-4,  1.0,  1.0*10-6,  1.7*105 )

       F(1) = 0.
       DO 2  I = 2, M
    2 F(I) = ALOG( 1. + SY(I) * SY(I) / (2. * SY(I) )
       Q = 1.E - 11
       CALL  EF13R (SX, SY, N, M, A, U, F, P, 1., 0., 0., Q, H, W, 4)
       CALL  EF13R (SX, SY, N, M, A, U, F, P, 1., 0., 0., Q, H, W, 2)

Результаты:

       U  = ( 0.04, 0.10, 0.19, 0.296, 0.39, 0.48, 0.59, 0.72, 0.85, 0.92, 0.84 )
       H  = ( 9.8*10-8, 1.7*10-2, 0.577, 3.*10-12, 4.6*10-6 )