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

Назначение

Решение одномерного интегрального уравнения первого рода с разностным ядром на прямой (- ∞, + ∞) методом регуляризации с применением алгоритма быстрого преобразования Фурье и с выбором регуляризации по невязке.

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

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

                   +∞
(1)    Ax  ≡   ∫  k(t - τ)  x(τ) dτ  =  y(t)  ,    -∞ < t < +∞ ,
                  -∞ 

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

(2)     Mα [x, y]  =  || Ax - y ||2 + α Ω[x]   =  min x ,
 где :
                                        +∞
       || Ax - y ||2  =  1/2π    ∫   | K(λ) X(λ) - Y(λ) |2 dλ  -  квадрат невязки,
                                       -∞
      α > 0 - параметр регуляризации,
                             +∞
       Ω[x]  =  1/2π  ∫  | λ |2p | X(λ) |2 dλ  -  стабилизирующий функционал
                            -∞
 порядка  p ≥ 0 (p не предполагается целым); 

    K (λ), X (λ), Y (λ) - преобразования Фурье ядра  k, решения  x и приближенной правой части  y.
Значение параметра  α определяется из условия

(3)     ρ(α) = ε || y || ,   ρ(α) ≡ || Axα - y || , 

где  x^α - регуляризированное решение задачи (2),
        ε - заданный относительный уровень невязки.

Считается, что  k (t) ,  x (t),  y (t) → 0 при  t → ±∞ и функции  k (t),  y (t) заданы на конечном отрезке [-Т/2, Т/2] в узлах равномерной сетки по  t :

     ts  =  s Δt ,   s =  -N/2, -N/2 + 1, ... , -1, 0, 1, ... , N/2 - 1 , 

где  N - четное число узлов,  Δ (t) = T/N - шаг сетки.

Численный алгоритм для решения уравнения типа свертки с дискретным ядром и правой частью методом регуляризации основан на применении дискретного преобразования Фурье временных рядов [2].

B связи с этим выполняются преобразования, заключающиеся в вычислении ДПФ функций  k_s = k (t_s) и  y_s = y (t_s), приводящие исходную точку к "каноническому" виду.

Приближенные значения  x_s^α регуляризованных решений в узлах  t_s определяются путем минимизации в спектральной области дискретной аппроксимации сглаживающего функционала и последующего выполнения обратного ДПФ. Для приближенного решения уравнения (3) используется метод касательных Ньютона в соответствии с вычислительной схемой [3].

Качество регуляризованного решения  x^α при значении параметра  α, удовлетворяющем критерию (3), оценивается следующими числовыми характеристиками :

1) невязкой    ρ(α)  =  || Axα - y ||  ,         

2) Ω - нормой решения     γ(α)  =  Ω1/2 [ xα ]  ,    

3) значением регуляризирующего функционала 
    φ(α)  =  [ ρ2(α) + αγ2(α) ]1/2  ,

4) функцией "чувствительности"    τ(α)  =  α Ω1/2 [ dxα/dα ] ,
    где  Ω [x] ≡ Ω [x]  при  p = 1/2
    (для вычисления "квазиоптимальных" значений
    параметра регуляризации),

5) найденным значением  α как решением уравнения (3),

6) числом итераций при решении уравнения (3),

7) практически достигнутым относительным уpовнем невязки  ε0. 

Вычисление всех критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) происходит в спектральной области  λ и требует порядка  N операций.

Возможны следующие режимы вычислений (в зависимости от значения параметра  L).

I.   Решение интегрального уравнения с приведением задачи к "каноническому" виду (L = 1).
II.  Решение интегрального уравнения при условии, что задача приведена к "каноническому" виду (L = 2).
III. Решение интегрального уравнения при условии, что ДПФ ядра вычислено (L = 3).

Эти режимы целесообразно использовать, например, при повторном решении того же интегрального уравнения (L = 2) или при решении уравнения с тем же ядром, но с другой правой частью (L = 3), а также в случае, когда известно аналитическое преобразование Фурье ядра  K (λ)  (L = 3).

Для оптимального вычисления прямого и обратного ДПФ применяется алгоритм быстрого преобразования Фурье, за счет которого реализованный алгоритм численного решения задачи является эффективным в смысле быстродействия, экономии памяти ЭВМ и влияния ошибок округления, при этом время вычислений пропорционально  N log₂ N, а объем памяти пропорционален  N.

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

    SUBROUTINE  EC01R (N, DT, A, Y, P, EPS, L, H, X, Z) 

   Параметры 

   Версии:  нет 

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

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

     A(I) = Re KI-1 ,         I = 1, 2, ..., N/2 + 1 ,
     A(I) = Im KI-N/2-1 ,   I = N/2 + 2, N/2 + 3, ..., N ; 

ρ(∞)  =  lim  ρ(α) ,
           α → ∞  

x∞  =  lim  xα  ≡  0
        α → ∞ 

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

Рассматривается решение интегрального уравнения типа свертки (1) с ядром  k (t) = exp (- t²) и правой частью  y (t) = (π/2) ^1/2 exp (- t²/2) ( точное решение  x (t) = exp (- t²) ) или y₁ (t) = (π/3) ^1/2 exp (- 2/3 t²) ( точное решение  x₁ (t) = exp (- 2t²) ).

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие значения критериальных функций при заданном относительном уpовне невязки  ε с использованием регуляризатора порядка  p = 1 и равномерной на отрезке [- 1, 1] сетки из  N = 8 точек.

       REAL  A(8), Y(8), H(7), X(8), Z(8), Y1(8)
       DATA  N /8/,  DT /0.25/,  P /1./,  PI /3.14159265359/
       DO 1  I = 1, N
       T = DT*(-N/2 + I - 1)
       A(I) = EXP(-T*T)
       Y(I) = SQRT(PI/2.) * EXP(-T*T/2.)
    1 Y1(I) = SQRT(PI/3.) * EXP(-T*T*2./3.)
       L = 1
       EPS = 0.08
       CALL  EC01R (N, DT, A, Y, P, EPS, L, H, X, Z)

Результаты:

       H  =  ( 0.122,  1.213,  0.164,  0.335,  0.008,  3,  0.080 )
       X  =  ( 0.339,  0.447,  0.712,  0.991,  1.113,  0.991,  0.712,  0.447 )

       L = 3
       EPS = 0.085
       CALL  EC01R (N, DT, A, Y1, P, EPS, L, H, X, Z)

Результаты:

       H  =  ( 0.102,  1.492,  0.149,  0.325,  0.005,  3,  0.085 )
       X  =  ( 0.094,  0.225,  0.549,  0.893,  1.046,  0.893,  0.549,  0.225 )