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

Назначение

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

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

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

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

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

(2)     Mα [x, y]  =  || Ax - y ||2 + α Ω[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.

 Функционал (2) минимизирует функция
                                +∞
(3)     xα(t)  =  1/2π  ∫  K*(λ) Y(λ) / ( | K(λ) |2 + α λ2p ) exp (iλt) dt ,
                              -∞
          -∞ < t < +∞ ,
 где  *  - знак комплексного сопряжения. 

Считается, что  k (t) → 0,  x (t) → 0,  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 - шаг сетки.

Задача состоит в решении дискретного уравнения типа свертки с ядром  k_s = k (t_s) и приближенной правой частью  y_s = = y (t_s) методом регуляризации  p - ого порядка в спектральной области  λ. Численный алгоритм основан на применении Дискретного Преобразования Фурье временных рядов [2].

Приведем задачу к "каноническому" виду, вычислив ДПФ функций  k_s и  y_s :

                         N/2-1
             Km   =    ∑      ks exp(-i λm ts) ,
                        S = -N/2
                         N/2-1
             Ym   =    ∑      ys exp(-i λm ts) ,
                       S = -N/2 

где  λ_m = m Δλ ,  m =  - N/2, - N/2 + 1 , ..., - 1, 0, 1, ..., N/2 - 1 ,  Δλ = 2π/(N Δt) ,  i = √-1.

Аналогично определяются  X_m,  X_m^α  -  ДПФ рядов  x_s = x (t_s)  и  x_s^α =x^α (t_s).

Далее, аппроксимируя интегралы в pавенствах (2), (3) и в преобразованиях Фурье по формуле прямоугольников, получаем дискретные аналоги сглаживающего функционала  M^α [x_s, y_s] и минимизирующего его регуляризованного решения

                                 N/2-1
(4)   xsα  = 1/ NΔt      ∑      KTmYm /( | Km |2 + α/Δt2 λm2p ) exp(2π i m s / N) ,
                               m= -N/2
        s  =  -N/2, -N/2 + 1, ..., N/2 - 1 . 

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

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

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

3) значение    регуляризующего    функционала 
    φ(α)  =  [ ρ2(α) + αγ2(α) ]1/2  ,
4) функция "чувствительности"
   (для  выбора  квазиоптимального значения  α)
    τ(α)  =  α Ω1/2 [ dxα/dα ] ,
   где
                           +∞
    Ω[x]  =  1/2π  ∫   | λ | | X(λ) |2 dλ  ≡  Ω[x]
                         -∞
    при  p = 1/2. 

Вычисление всех критериальных функций происходит в спектральной области  λ (через компоненты ДПФ). Это требует порядка  N операций и значительно облегчает задачу выбора параметра регуляризации. Восстановление искомого решения по формуле (4) целесообразно производить лишь для выбранных значений  λ.

Подпрограмма предусматривает работу в следующих режимах (в зависимости от значения параметра  L).

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

Полное описание изложенного алгоритма содержится в статье [3].

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

    SUBROUTINE  EC00R (N, DT, A, Y, ALP, P, 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 ; 

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

Рассматривается решение интегрального уравнения типа свертки (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²) ).

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

       REAL  A(8), Y(8), H(4), X(8), Z(8), Y1(8)
       DATA  N /8/,  DT /0.25/,  ALP /1.E - 2/,  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
       CALL  EC00R (N, DT, A, Y, ALP, P, L, H, X, Z)
       L = 2
       CALL  EC00R (N, DT, A, Y, ALP, P, L, H, X, Z)

Результаты:

       H  = ( 0.13239,  1.08822,  0.17138,  0.33298 )
       X  = ( 0.378,  0.475,  0.713,  0.963,  1.072,  0.963,  0.713,  0.475 )


       L = 4
       CALL  EC00R (N, DT, A, Y1, ALP, P, L, H, X, Z)
       L = 2
       CALL  EC00R (N, DT, A, Y1, ALP, P, L, H, X, Z)

Результаты:

       H  = ( 0.13270,  1.11089,  0.17306,  0.33989 )
       X  = ( 0.211,  0.310,  0.554,  0.809,  0.919,  0.809,  0.554,  0.310 )