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

Назначение

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

Значение параметра регуляризации  α определяется из условия:

          ρ(α)  =  e || f ||p  ,
   где
                        d           b
       ρ(α)  =  (  ∫ p(y) [  ∫ k(y, x) uα(x) dx - f(y) ]2 dy )1/2   -
                      c            a 

невязка на регуляризированном решении  u_α (x),  e - заданный относительный уровень невязки,

                         d        
       || f ||p  =  (  ∫  p(y) f 2(y) dy )1/2   -   норма правой части уравнения (1).
                       c 

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

  Тогда дискретный аналог (2) имеет вид:
         M                N
(3)     ∑  pi σi''  [  ∑  σj' ki j uj - fi ]2  +  α u T C 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 означает операцию транспонирования. 

При этом параметр  α выбирается из дискретного аналога условия невязки:

(5)     ρ(α) ≈ || A uα - f ||  =  e || f || , 

где  u_α - решение (4).

При решении системы (4) используется сведение задачи к более простой "канонической" задаче с единичной матрицей  С и двухдиагональной матрицей  A согласно методу В.В.Воеводина [2].

Для решения уравнения (5) используется метод касательных Ньютона в соответствии с вычислительной схемой, предложенной В.А.Морозовым [3].

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

1) невязка  ρ(α);
2) функция "чувствительности"
    τ(α) ≡ α Ω1/2 (duα / dα) ;
3) Ω - ноpма регуляризованного решения
    γ(α) ≈ Ω1/2 (uα) ;
4) значение   регуляризирующего   функционала 
    φ(α) = ρ2(α) + αγ2(α) ;
5) найденное значение  α как решение уравнения (5);
6) число итераций при решении уравнения (5);
7) достигнутый относительный уровень невязки
    e = ρ(α) / || f || . 

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

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

I) Решение интегрального уравнения для фиксированный правой части и относительного уровня невязки (это требует порядка  MN² операций).

II) Pешение того же уравнения, но с другим относительным уровнем невязки (это требует порядка  N² операций).

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

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

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

   Параметры 

   Версии:  нет 

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

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

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

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

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

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

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие значения критериальных функций при заданных значениях относительного уровня невязки  EPS с весовой функцией  P (y) = 1 и использованием равномерных по  X и  Y на отрезке [0, 1] сеток из 11 точек.

       REAL  SX(11), SY(11), P(11), U(11), F(11), A(11, 11), H(7), W(11, 8)
       DATA  P /11*1./, N /11/, M /11/, EPS /0.001/
       DATA  SX /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  EF12R (SX, SY, N, M, A, U, F, P, 1., 0., 0., EPS, H, W, 1)

Результаты:

       U  =  ( 1.01, 1.01, 1.01, 1.01, 1.01, 1., 1., 0.99, 0.97, 0.96, 0.94 )
       H  =  ( 4.96*10-4, 1.46*10-2, 0.99, 1.07*10-4, 1.07*10-4, 1., 1.0*10-3 )

       F(1) = 0.
       DO 2  I = 2, N
    2 F(I) = ALOG( 1. + SY(I)*SY(I) ) / ( 2.*SY(I) )
       CALL  EF12R (SX, SY, N, M, A, U, F, P, 1., 0., 0., EPS, H, W, 3)

Результаты:

       U  =  ( 0.11, 0.13, 0.18, 0.27, 0.37, 0.48, 0.60, 0.70, 0.80, 0.88, 0.94 )
       H  =  ( 2.3*10-4, 8.54*10-3, 0.57, 5.85*10-6, 1.78*10-5, 1., 1.0*10-3 )