Текст подпрограммы и версий ( Фортран )
ec12c.zip 		Тексты тестовых примеров ( Фортран )
tec12c.zip
Текст подпрограммы и версий ( Си )
ec12c_c.zip 		Тексты тестовых примеров ( Си )
tec12c_c.zip
Текст подпрограммы и версий ( Паскаль )
ec12c_p.zip 		Тексты тестовых примеров ( Паскаль )
tec12c_p.zip

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

Назначение

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

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

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

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

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

Предполагается, что известно аналитическое преобразование Фурье ядра  K (λ, ω) и правая часть уравнения (1)  g (x, y) задана на прямоугольнике:  x ∈ [- T1/2, T1/2),  y ∈ [- T2/2, T2/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 ,

где  N1, N2 - четные числа узлов,  d1 = T1/N1,  d2 = T2/N2 - шаги сетки по переменным  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 = N1/2, n2 = N2/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 ) ]  }

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

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

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

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

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

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

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

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

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

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

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

    SUBROUTINE  EC12C (F, GR, GI, N1, N2, D1, D2, ALP, P, L, H,
                                            FR, FI) 
   Параметры
F(W1,W2)- комплексная подпрограмма - функция вычисления преобразования Фурье ядра  K (λ, ω) в точке  λ = W1,  ω = W2;
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ок матрицы правой части  Gs1, s2,  N1 ≥ 4 - целая степень числа два (тип: целый);
N2 - заданное число столбцов матрицы правой части  Gs1, s2,  N2 ≥ 4 - целая степень числа два (тип: целый);
D1 - заданный шаг сетки по переменным  x, ξ, D1 > 0 (тип: вещественный);
D2 - заданный шаг сетки по переменным  y, η, D2 > 0 (тип: вещественный);
ALP - заданное значение параметра регуляризации  α,  ALP ≥ 0 (тип: вещественный);
P - заданное значение порядка регуляризатора  p,  P ≥ 0 (тип: вещественный);
L - параметр, определяющий режим использования подпрограммы (тип: целый);
L = 1 - вычисление ДПФ правой части уравнения (1),
L = 2 - построение решения и вычисление значений критериальных функций в предположении, что вещественные и мнимые части элементов ДПФ правой части содержатся соответственно в массивах GR, GI,
L = 3 - вычисление только значений критериальных функций в том же предположении, что при  L = 2;
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 - подпрограмма вычисления двумерного дискретного или обратного дискретного преобразования Фурье комплексной матрицы размера  N1 * N2,  N1 = 2 J1,  N2 = 2 J2 методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

B результате работы подпрограммы при  L = 1 в массиве  GR содержится вещественная часть, а в массиве  GI мнимая часть элементов ДПФ правой части.

  2. 

Комплексная подпрограмма - функция  F (W1, W2), вычисляющая значения преобразования Фурье ядра  K (λ, ω) в точке  λ = W1,  ω = W2, должна быть написана пользователем.

  3. 

Массивы  GR, GI используются при всех значениях  L, а массивы  FR, FI - только при  L = 2 или 3.

  4.  Если каждый раз при обращении к подпрограмме EC12C ДПФ правой части вычислять заново, то общее число массивов, используемых для решения задачи, можно сократить, совмещая параметры  GR и  FR,  GI и  FI.

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

Рассматривается задача решения интегрального уравнения (1) с ядром  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) сетки из N1*N2 точек, где  N1 = 8 , N2 = 8. 

       REAL  GR(8, 8), GI(8, 8), H(4), FR(8, 8), FI(8, 8)
       DATA  N1 /8/, N2 /8/, D1 /0.25/, D2 /0.25/, ALP /3.E - 02/, P /1./
       EXTERNAL  F
       DO 1  I = 1, N1
       DO 1  J = 1, N2
       T1 = D1 * (-N1/2 + I - 1)
       T2 = D2 * (-N2/2 + J - 1)
       GR(I, J) = 1.57079632679 * EXP( -(T1*T1 + T2*T2)/2. )
    1 GI(I, J) = 0.
       L = 1
       CALL  EC12C (F, GR, GI, N1, N2, D1, D2, ALP, P, L, H, FR, FI)
       L = 2
       CALL  EC12C (F, GR, GI, N1, N2, D1, D2, ALP, P, L, H, FR, FI)
       END

       COMPLEX  FUNCTION  F(W1, W2)
       F = CMPLX( EXP( -0.25 * (W1*W1 + W2*W2) ) * 3.14159265359, 0. )
       RETURN
       END

Результаты:

       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.