Текст подпрограммы и версий
ec02c_c.zip 		Тексты тестовых примеров
tec02c_c.zip

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

Назначение

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

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

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

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

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

(2)     M α [f, g]  =  || Af - g ||2 + α Ω [f] ,
   где
                                          +∞   +∞
     || Af - g ||2  =  1/(2π)2     ∫     ∫   | K(λ, ω) F(λ, ω) - G(λ, ω) |2 dλ dω -
                                         -∞   -∞
     квадрат невязки,
     α > 0 - параметр регуляризации,
                                 +∞   +∞
     Ω[f]  =  1/(2π)2     ∫     ∫   [ 1 + (λ2 + ω2)P ] | F(λ, ω) |2 dλ dω -
                               -∞   -∞
     стабилизирующий функционал порядка  P ≥ 0
     (P не предполагается целым);

    K(λ, ω), F(λ, ω), G(λ, ω) - преобразование Фурье ядра  k, решения  f и приближенной правой части  g. 

 Функционал (2) минимизирует функция
                                        +∞   +∞
(3)  f α(ξ, η)  =  1/(2π)2     ∫     ∫    [ K*(λ, ω) G(λ, ω) /
                                       -∞   -∞
            / | K(λ, ω) |2 + α[1 + (λ2 + ω2)P ] ]  e i (λξ + ωη) dλ dω

                            -∞ < ξ, η < +∞ ,

где * - знак комплексного сопряжения,  i = √-1.

Считается, что  k (x, y)→0,  f (x, y)→0,  g (x, y)→0 при  x→±∞,  y→±∞ и функции  k (x, y), g (x, y) заданы на конечном прямоугольнике:

x ∈ [- T1/2, T1/2 ),  y ∈ [- T2/2, T2/2 ) в узлах одной и той же равномерной сетки по  x, y: 

     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,
где  N1 - четное число узлов,  d1 = T1/N1 - шаг сетки по  x,
        N2 - четное число узлов,  d2 = T2/N2 - шаг сетки по  y.

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

B связи с этим выполняются преобразования, заключающиеся в вычислении ДПФ функций  ks1, s2 = k (xs1, ys2) и  gs1, s2 = g (xs1, ys2), приводящие исходную задачу к "каноническому" виду ( n1 = N1/2,  n2 = N2/2 ): 

                       n1-1       n2-1
    Km1,m2 =     ∑           ∑        ks1,s2 exp( -i (λ m1 xs1 + ωm2 ys2) )  =
                     s1= -n1  s2= -n2

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

                       n1-1       n2-1
    Gm1,m2 =     ∑           ∑        gs1,s2 exp( -i (λ m1 xs1 + ωm2 ys2) )  =
                     s1= -n1  s2= -n2

                       n1-1       n2-1
                =     ∑           ∑        gs1,s2 exp( -2π i (s1 m1 / N1 + s2 m2 / N2) ) ,
                     s1= -n1  s2= -n2
 где
        λm1 = m1 Δλ ,     m1 = -N1/2, -N1/2 + 1,...,-1,0,1,..., N1/2 - 1,
        Δλ = 2π / N1d1 ,
        ωm2 = m2 Δω ,    m2 = -N2/2, -N2/2 + 1,...,-1,0,1,..., N2/2 - 1,
        Δω = 2π / N2d2 .

Аналогично определяются  Fm1, m2,  F αm1, m2 - ДПФ матриц с элементами  fs1, s2 = f (xs1, ys2)  и  f αs1, s2 = f α (xs1, ys2).

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

                                                     n1-1          n2-1
(4)  f αs1, s2  =  1 / N1N2d1d2        ∑             ∑
                                                   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 ] / d12d22 ]  }

B формуле (4) предполагается, что при  P = 0,  m1 = 0,  m2 = 0    (λ2m1 + ω2m2) P = 1.

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

  1) невязка  ρ(α)  =  || Af α - g ||  ,
  2) Ω - ноpма решения  γ(α)  =  Ω1/2 [f α] ,

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

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

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

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

1.  L = 1. Задача полностью приводится к "каноническому" виду, т.е. производится вычисление ДПФ каждой из двух исходных матриц с элементами  ks1, s2,  gs1, s2.
2.  L = 2. Решается интегральное уравнение и вычисляются значения критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) для заданного значения  α - при условии, что задача приведена к "каноническому" виду.
3.  L = 3. Вычисляются только значения критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) для заданного значения  без построения решения (4) - при условии, что задача приведена к "каноническому" виду.
4.  L = 4. Задача приводится к "каноническому" виду в предположении, что ДПФ ядра  Km1, m2 уже вычислено.

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

Для эффективного вычисления прямого и обратного ДПФ применяется двумерное быстрое преобразование Фурье [2]. При этом время вычислений пропорционально  N1 N2 (log2 N1 + log2 N2), а объем памяти пропорционален  N1 N2.

Аналогичный алгоритм для решения одномерного интегрального уравнения типа свертки описан в [3].

1.  А.Н.Тихонов. O регуляризации некоppектно поставленных задач, ДАН CCCP, 1963, т.153, N 1, 49 - 52.
2.  В.А.Морозов, Н.Н.Кирсанова, А.Ф.Сысоев. Комплекс алгоритмов быстрого преобразования Фурье дискретных рядов, сб. "Численный анализ на ФОРТРАНе", вып.15, Изд - во МГУ, M., 1976.
3.  M.B. Арефьева. Быстрый регуляризирующий алгоритм решения интегральных уравнений Фредгольма I рода типа свертки, сб. "Численный анализ на ФОРТРАНе. Методы и алгоритмы", Изд - во МГУ, M., 1978.

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

    int ec02c_c(real *ar, real *ai, real *gr, real *gi,
        integer *n1, integer *n2, real *d1, real *d2, real *alp, real *p,
        integer *l, real *h, real *fr, real *fi)

   Параметры
ar, ai - двумерные вещественные массивы размера  n1 * n2, содержащие соответственно вещественные и мнимые части элементов заданной матрицы ядра на сетке: 
    ar(i, j)  =  re k I-n1/2-1, J-n2/2-1,
    ai(i, j)  =  im k I-n1/2-1, J-n2/2-1,
     i = 1, 2, ..., n1,  j = 1, 2, ..., n2;
gr, gi - двумерные вещественные массивы размера  n1 * n2, содержащие соответственно вещественные и мнимые части элементов заданной матрицы правой части на сетке: 
    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;
n1 - заданное число стpок матриц ядра  ks1, s2 и правой части  gs1, s2,  n1 ≥ 4 - целая степень числа два (тип: целый);
n2 - заданное число столбцов матриц ядра  ks1, s2 и правой части  gs1, s2,  n1 ≥ 4 - целая степень числа два (тип: целый);
d1 - заданный шаг сетки по переменным  x, ξ, d1 > 0 (тип: вещественный);
d2 - заданный шаг сетки по переменным  y, η, d2 > 0 (тип: вещественный);
alp - заданное значение параметра регуляризации  α,  alp ≥ 0 (тип: вещественный);
p - заданное значение порядка регуляризации  p,  p ≥ 0 (тип: вещественный);
l - параметр, определяющий режим использования подпрограммы (тип: целый):
l = 1 - полное приведение задачи к "каноническому" виду;
l = 2 - построение решения и вычисление значений критериальных функций в предположении, что вещественные и мнимые части элементов ДПФ ядра и правой части содержатся соответственно в массивах  ar, ai, gr, gi;
l = 3 - вычисление только значений критериальных функций в том же предложении, что при  l = 2;
l = 4 - приведение задачи к "каноническому" виду в предположении, что вещественные и мнимые части элементов ДПФ ядра содержатся соответственно в массивах  ar, ai;
h - вещественный вектоp длины 4, содержащий вычисленные значения критериальных функций при заданном значении  alp:  h (1) = ρ (α),  h (2) = γ (α),  h (3) = φ (α),  h (4) = τ (α);
fr, fi - двумерные вещественные массивы размера  n1 * 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. 
   Версии:  нет 
   Вызываемые подпрограммы
ftftc_c - вычисление двумерного дискретного или обратного дискретного преобразования Фурье комплексной матрицы размера  n1 * n2,  n1 =2 j1,  n2 = 2 j2 методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

B результате работы подпрограммы при  l = 1 в массивах  ar, gr содержатся вещественные части, а в массивах  ai, gi мнимые части соответственно элементов ДПФ ядpа  km1, m2 и правой части  gm1, m2.

  2. 

При первом обращении к подпрограмме (l = 1) необходимо задать значения параметров:  ar, ai, gr, gi, n1, n2, d1, d2, определяющих уравнение.

B дальнейшем, например, при повторном решении того же интегрального уравнения можно менять только параметры  alp, p, а при решении уравнения с тем же ядpом, но с другой правой частью - только параметры  gr, gi, alp, p, сохраняя значения всех остальных паpаметpов.

  3. 

Массивы  ar, ai, gr, gi используются при всех значениях  l, а массивы  fr, fi - только при  l = 2, 3.

  4.  Если каждый раз при обращении к подпрограмме ДПФ правой части  gm1, m2 вычислять заново, то общее число массивов, используемых для решения задачи, можно сократить, допуская совпадение параметров  gr, fr и  gi, fi.

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

Рассматривается решение интегрального уравнения типа свертки ( с ядром, определяемым функцией 

   k(x, y)  =  exp [-(x2 + y2)] ,
   и  правой  частью   g(x, y)  =  π/2 exp [-1/2 (x2 + y2)].
   Точное решение имеет вид      f(ξ, η)  =  exp [- (ξ2 + η2)].

Ниже приводится фрагмент программы, вычисляющей решение и соответствующие значения критериальных функций при фиксированном параметре регуляризации  α = 3 * 10- 2 с использованием регуляризатора порядка  p = 1 и равномерной на квадрате [- 1, 1) * [- 1, 1) сетки из  n1 * n2 точек, где  n1 = 8,  n2 = 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 */
    static float h__[4];
    static int i__, i, j, l;
    static float t, t1, t2,
                 ai[64]  /* was [8][8] */, fi[64]  /* was [8][8] */,
                 gi[64]  /* was [8][8] */, ar[64]  /* was [8][8] */,
                 fr[64]  /* was [8][8] */, gr[64]  /* was [8][8] */,
                 xi[64]  /* was [8][8] */, xr[64]  /* was [8][8] */;
    extern int ec02c_c(float *, float *, float *, float *, int *, int *,
            float *, float *, float *, float *, int *, float *, float *,
            float *);
    int i__1, i__2;

#define ai_ref(a_1,a_2) ai[(a_2)*8 + a_1 - 9]
#define fi_ref(a_1,a_2) fi[(a_2)*8 + a_1 - 9]
#define gi_ref(a_1,a_2) gi[(a_2)*8 + a_1 - 9]
#define ar_ref(a_1,a_2) ar[(a_2)*8 + a_1 - 9]
#define fr_ref(a_1,a_2) fr[(a_2)*8 + a_1 - 9]
#define gr_ref(a_1,a_2) gr[(a_2)*8 + a_1 - 9]
#define xi_ref(a_1,a_2) xi[(a_2)*8 + a_1 - 9]
#define xr_ref(a_1,a_2) xr[(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);
            t = t1 * t1 + t2 * t2;
            ar_ref(i__, j) = (float)exp((float)(-t));
            ai_ref(i__, j) = 0.f;
            gr_ref(i__, j) = (float)exp((float)(-t / 2.f)) * 1.57079632679f;
            gi_ref(i__, j) = 0.f;
            xr_ref(i__, j) = (float)exp((float)(-t));
/* l1: */
            xi_ref(i__, j) = 0.f;
        }
    }
    for (i = 1; i <= 8; ++i) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e \n",
            ar_ref(i, 1), ar_ref(i, 2), ar_ref(i, 3), ar_ref(i, 4),
            ar_ref(i, 5), ar_ref(i, 6), ar_ref(i, 7), ar_ref(i, 8));
    }
    printf("\n");
    for (i = 1; i <= 8; ++i) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e \n",
            ai_ref(i, 1), ai_ref(i, 2), ai_ref(i, 3), ai_ref(i, 4),
            ai_ref(i, 5), ai_ref(i, 6), ai_ref(i, 7), ai_ref(i, 8));
    }
    printf("\n");
    for (i = 1; i <= 8; ++i) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e \n",
            gr_ref(i, 1), gr_ref(i, 2), gr_ref(i, 3), gr_ref(i, 4),
            gr_ref(i, 5), gr_ref(i, 6), gr_ref(i, 7), gr_ref(i, 8));
    }
    printf("\n");
    for (i = 1; i <= 8; ++i) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e \n",
            gi_ref(i, 1), gi_ref(i, 2), gi_ref(i, 3), gi_ref(i, 4),
            gi_ref(i, 5), gi_ref(i, 6), gi_ref(i, 7), gi_ref(i, 8));
    }
    printf("\n");
    for (i = 1; i <= 8; ++i) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e \n",
            xr_ref(i, 1), xr_ref(i, 2), xr_ref(i, 3), xr_ref(i, 4),
            xr_ref(i, 5), xr_ref(i, 6), xr_ref(i, 7), xr_ref(i, 8));
    }
    printf("\n");
    for (i = 1; i <= 8; ++i) {
    printf("\n %16.7e %16.7e %16.7e %16.7e \n %16.7e %16.7e %16.7e %16.7e \n",
            xi_ref(i, 1), xi_ref(i, 2), xi_ref(i, 3), xi_ref(i, 4),
            xi_ref(i, 5), xi_ref(i, 6), xi_ref(i, 7), xi_ref(i, 8));
    }
    l = 1;
    ec02c_c(ar, ai, gr, gi, &n1, &n2, &d1, &d2, &alp, &p, &l, h__, fr, fi);
    l = 2;
    ec02c_c(ar, ai, gr, gi, &n1, &n2, &d1, &d2, &alp, &p, &l, h__, fr, fi);

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



Результаты:

       h__  =  ( 0.328307,  1.652517,  0.435557,  0.828122 ) ,

                | 0.133,  0.186,  0.317,  0.454,  0.514,  0.454,  0.317,  0.186 |
                | 0.186,  0.240,  0.372,  0.510,  0.571,  0.510,  0.372,  0.240 |
                | 0.317,  0.372,  0.508,  0.649,  0.710,  0.649,  0.508,  0.372 |
       fr   = | 0.454,  0.510,  0.649,  0.793,  0.856,  0.793,  0.649,  0.510 |
                | 0.514,  0.571,  0.710,  0.856,  0.920,  0.856,  0.710,  0.571 |
                | 0.454,  0.510,  0.649,  0.793,  0.856,  0.793,  0.649,  0.510 |
                | 0.317,  0.372,  0.508,  0.649,  0.710,  0.649,  0.508,  0.372 |
                | 0.186,  0.240,  0.372,  0.510,  0.571,  0.510,  0.372,  0.240 |

       fi - все значения имеют порядки 10-12 - 10-17.