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

Назначение

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

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

Приближенное решение одномepнoгo интегрального уравнения 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.

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

    int ec01r_c (integer *n, real *dt, real *a, real *y, real *p,
            real *eps, integer *l, real *h, real *x, real *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
        α → ∞ 

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

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

int main(void)
{
    /* Initialized data */

    static int n = 8;
    static float dt = .25f;
    static float p = 1.f;
    static float pi = 3.14159265f;

    /* Builtin functions */
    double exp(double), sqrt(double);

    /* Local variables */
    extern int ec01r_c(int *, float *, float *, float *,
                       float *, float *, int *, float *, float *, float *);
    static float a[8], h__[7];
    static int i__, i, l;
    static float t, x[8], y[8], z__[8], y1[8], xt[8], xt1[8], eps;
    int i__1;

    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        t = dt * (-n / 2 + i__ - 1);
        a[i__ - 1] = (float)exp((float)(-t * t));
        xt[i__ - 1] = (float)exp((float)(-t * t));
        y[i__ - 1] = (float)sqrt((float)(pi / 2.f)) *
                (float)exp((float)(-t * t / 2.f));
        xt1[i__ - 1] = (float)exp((float)(t * -2.f * t));
/* l1: */
        y1[i__ - 1] = (float)sqrt((float)(pi / 3.f)) *
                 (float)exp((float)(-t * t * 2.f / 3.f));
    }
    for (i = 0; i <= 6; i+= 2) {
        printf("\n %16.7e %16.7e \n",a[i], a[i+1]);
    }
    for (i = 0; i <= 6; i+= 2) {
        printf("\n %16.7e %16.7e \n",xt[i], xt[i+1]);
    }
    for (i = 0; i <= 6; i+= 2) {
        printf("\n %16.7e %16.7e \n",y[i], y[i+1]);
    }
    l = 1;
    eps = .08f;
    printf("\n\n %16.7e \n",eps);
    ec01r_c(&n, &dt, a, y, &p, &eps, &l, h__, x, z__);

    for (i = 0; i <= 4; i+= 2) {
        printf("\n %16.7e %16.7e ",h__[i], h__[i+1]);
    }
        printf("\n %16.7e \n",h__[6]);
    for (i = 0; i <= 6; i+= 2) {
        printf("\n %16.7e %16.7e \n",x[i], x[i+1]);
    }
    for (i = 0; i <= 6; i+= 2) {
        printf("\n %16.7e %16.7e \n",xt1[i], xt1[i+1]);
    }
    for (i = 0; i <= 6; i+= 2) {
        printf("\n %16.7e %16.7e \n",y1[i], y1[i+1]);
    }
    l = 3;
    eps = .085f;
    printf("\n\n %16.7e \n",eps);
    ec01r_c(&n, &dt, a, y1, &p, &eps, &l, h__, x, z__);

    for (i = 0; i <= 4; i+= 2) {
        printf("\n %16.7e %16.7e ",h__[i], h__[i+1]);
    }
        printf("\n %16.7e \n",h__[6]);
    for (i = 0; i <= 6; i+= 2) {
        printf("\n %16.7e %16.7e ",x[i], x[i+1]);
    }
    return 0;
} /* main */


Результаты:

       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 = .085f;
       ec01r_c(&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 )