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

Назначение

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

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

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

                  +∞   +∞
(1)    Af  ≡   ∫     ∫   K(x - ξ, y - η) f(ξ, η) dξ dη = g(x, y)  ,
                 -∞    -∞
                                   -∞ < x, y < +∞ 

с гладким ядром K осуществляется методом однопараметрической регуляризации А.Н.Тихонова [1] путем сведения задачи к минимизации соответствующего сглаживающего функционала (см. описание подпрограммы ec02c_c в данной Библиотеке).

Предполагается, что известно аналитическое преобразование Фурье ядра  K (λ, ω) и правая часть уравнения (1)  g (x, y) задана на прямоугольнике:  x ∈ [- T₁/2, T₁/2),  y ∈ [- T₂/2, T₂/2) в узлах равномерной сетки

xs1 = s1*d1 ,     s1 = - N1/2, - N1/2 + 1, ...,-1, 0, 1, ..., N1/2 - 1
ys2 = s2*d2 ,     s2 = - N2/2, - N2/2 + 1, ...,-1, 0, 1, ..., N2/2 - 1 , 

где  N₁, N₂ - четные числа узлов,  d1 = T₁/N₁,  d2 = T₂/N₂ - шаги сетки по переменным  x, y.

Аналогично введем сетку по переменным  λ,  ω в частотной области

     λm1 = m1 * Δλ   ,   m1 = - N1/2, - N1/2 + 1, ...,-1, 0, 1, ..., N1/2 - 1 ,
     Δλ = 2π / N1 d1  ,
     ωm2 = m2 * Δω  ,     m2 = - N2/2, - N2/2 + 1, ...,-1, 0, 1, ..., N2/2 - 1 ,
     Δω = 2π / N2 d2 

и рассмотрим дискретные значения преобразования Фурье ядра  К_m1, m2 = K (λ_m1, ω_m2).

Далее, аппроксимируя все интегралы по формуле прямоугольников, получаем дискретные аналоги преобразования Фурье правой части ( n1 = N₁/2, n2 = N₂/2 )

                                 n1-1       n2-1
(2) Gm1,m2 = d1d2     ∑           ∑        gs1,s2 exp( -2π i (s1 m1 /N1 + s2 m2 /N2))
                               s1= -n1  s2= -n2

 и регуляризованного решения

                                                     n1-1          n2-1
(3)  f αs1, s2  =  1 / N1N2d1 d2      ∑             ∑
                                                     m1= -n1    m2= -n2
                  { K*m1,m2 Gm1,m2  exp( 2π i (s1 m1 /N1 + s2 m2 /N2) ) /
                    / [ | Km1,m2 |2 + α ( 1 + (λm12 + ωm22)p ) ]  } 

где  G_s1, s2 = G (x_s1, y_s2),  α > 0 - параметр регуляризации,  p ≥ 0 - порядок стабилизирующего функционала (p не предполагается целым).

B формуле (3) при  p = 0,  m1 = 0,  m2 = 0 полагается ( λ_m1² + ω_m2² )^p = 1.

Для эффективного вычисления сумм вида (2), (3) применяется алгоритм быстрого дискретного преобразования Фурье.

Подпрограмма реализует изложенный алгоритм решения уравнения (1) и вычисляет приближенные значения четырех критериальных функций - невязки  ρ (α), нормы решения  γ (α), регуляризирующего функционала  φ (α), функции "чувствительности"  τ (α) (см. описание подпрограммы ec02c_c), которые могут быть использованы для выбора параметра регуляризации  α.
Вычисление всех критериальных функций происходит через компоненты ДПФ, что требует порядка  N₁ N₂ операций и значительно облегчает задачу выбора  α.
Выполнение обратного ДПФ в формуле (3) целесообразно производить лишь для выбранных значений  α.
За счет применения БПФ полное время численного решения задачи пропорционально  N₁ N₂ (log₂ N₁ + log₂ N₂), а объем памяти ЭВМ пропорционален  N₁ N₂.

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

L = 1 - задача приводится к "каноническому" виду, т.е. производится вычисление ДПФ правой части уравнения (1),

L = 2 - решается интегральное уравнение и вычисляются значения критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) для заданного значения  α - при условии, что задача приведена к "каноническому" виду,

L = 3 - вычисляются только значения критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) для заданного значения  α без построения решения (3) - при условии, что задача приведена к "каноническому" виду.

Эти режимы целесообразно использовать, например, при повторном решении того же интегрального уравнения или при выборе значения параметра регуляризации  α.
Аналогичный алгоритм описан в [3].
Отметим, что по сравнению с подпрограммой ec02c_c в данной подпрограмме экономятся два двумерных массива размера  N₁ * N₂.

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

    int ec12c_c (C_fp f, real *gr, real *gi, integer *n1,
            integer *n2, real *d1, real *d2, real *alp, real *p, integer *l,
            real *h, real *fr, real *fi)

   Параметры 

      gr(i, j) = re g I-n1/2-1,  J-n2/2-1 ,

      gi(i, j) = im g I-n1/2-1,  J-n2/2-1 ,

      i = 1, 2, ..., n1,   j = 1, 2, ..., n2;  

      fr(i, j) = re f αI-n1/2-1, J-n2/2-1 ,

      fi(i, j) = im f αI-n1/2-1, J-n2/2-1 ,

      i = 1, 2, ..., n1 ,   j = 1, 2, ..., n2; 

   Версии:  нет 

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

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

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

Рассматривается задача решения интегрального уравнения ( с ядром  k (x, y), аналитическое преобразование Фурье которого имеет вид

              k(λ, Ω) = π exp[ -1/4 (λ2 + ω2) ] ,
   и правой частью
              g(x, y) = π/2 exp[ -1/2 (x2 + y2) ]

   ( точное решение   f(ξ, η) = exp [ - (ξ2 + η2) ] ) .

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

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

    static int n1 = 8;
    static int n2 = 8;
    static float d1 = .25f;
    static float d2 = .25f;
    static float alp = .03f;
    static float p = 1.f;

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

    /* Local variables */
    extern void f_c();
    static float h__[4];
    static int i__, i, j, l;
    static float t1, t2, fi[64]  /* was [8][8] */, gi[64] /* was [8][8] */,
                         fr[64]  /* was [8][8] */, gr[64] /* was [8][8] */;
    extern int ec12c_c(U_fp, float *, float *, int *, int *, float *,
            float *, float *, float *, int *, float *, float *, float *);
    int i__1, i__2;

#define gi_ref(a_1,a_2) gi[(a_2)*8 + a_1 - 9]
#define gr_ref(a_1,a_2) gr[(a_2)*8 + a_1 - 9]

    i__1 = n1;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n2;
        for (j = 1; j <= i__2; ++j) {
            t1 = d1 * (-n1 / 2 + i__ - 1);
            t2 = d2 * (-n2 / 2 + j - 1);
            gr_ref(i__, j) = (float)exp((float)(-(t1 * t1 + t2 * t2) / 2.f)) *
                    1.57079632679f;
/* l1: */
            gi_ref(i__, j) = 0.f;
        }
    }
    l = 1;
    ec12c_c((U_fp)f_c, gr, gi, &n1, &n2, &d1, &d2, &alp, &p, &l, h__, fr, fi);
    l = 2;
    ec12c_c((U_fp)f_c, gr, gi, &n1, &n2, &d1, &d2, &alp, &p, &l, h__, fr, fi);

    printf("\n %16.7e %16.7e %16.7e %16.7e \n",
           h__[0], h__[1], h__[2], h__[3]);
    for (i = 1; i <= 58; i += 8) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e ",
            fr[i-1], fr[i], fr[i+1], fr[i+2],
            fr[i+3], fr[i+4], fr[i+5], fr[i+6]);
    }
    printf("\n");
    for (i = 1; i <= 58; i += 8) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e ",
            fi[i-1], fi[i], fi[i+1], fi[i+2],
            fi[i+3], fi[i+4], fi[i+5], fi[i+6]);
    }
    printf("\n");
    return 0;
} /* main */


/* Complex */ void f_c(complex * ret_val, float *w1, float *w2)
{
    /* System generated locals */
    float r__1;
    complex q__1;

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

    r__1 = (float)exp((float)((*w1 * *w1 + *w2 * *w2) * -.25f)) *
            3.14159265359f;
    q__1.r = r__1, q__1.i = 0.f;
     ret_val->r = q__1.r,  ret_val->i = q__1.i;
    return ;
} /* f_c */


Результаты:

       h__  =  ( 0.406963,  1.260784,  0.461851,  0.847403 )

                | 0.051,  0.096,  0.206,  0.315,  0.360,  0.315,  0.206,  0.096 |
                | 0.096,  0.142,  0.252,  0.361,  0.407,  0.361,  0.252,  0.142 |
                | 0.206,  0.252,  0.362,  0.473,  0.519,  0.473,  0.362,  0.252 |
       fr   = | 0.315,  0.361,  0.473,  0.584,  0.631,  0.584,  0.473,  0.361 |
                | 0.360,  0.407,  0.519,  0.631,  0.677,  0.631,  0.519,  0.407 |
                | 0.315,  0.361,  0.473,  0.584,  0.631,  0.584,  0.473,  0.361 |
                | 0.206,  0.252,  0.362,  0.473,  0.519,  0.473,  0.362,  0.252 |
                | 0.096,  0.142,  0.252,  0.361,  0.407,  0.361,  0.252,  0.142 |
 
       fi - все значения имеют порядки  10-12 - 10-14.