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

Подпрограмма:  EC12C (модуль EC12C_p)

Назначение

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

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

procedure EC12C(F :Func_F_E; var GR :Array of Real; var GI :Array of Real;
                N1 :Integer; N2 :Integer; D1 :Real; D2 :Real;
                ALP :Real; P :Real; L :Integer; var H :Array of Real;
                var FR :Array of Real; var FI :Array of Real);

Параметры

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.

Unit tec12c_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FEC12C_p, EC12C_p;

function tec12c: String;

implementation

function tec12c: String;
var
I,J,L,_i :Integer;
T1,T2 :Real;
GR :Array [0..63] of Real;
GI :Array [0..63] of Real;
H :Array [0..3] of Real;
FR :Array [0..63] of Real;
FI :Array [0..63] of Real;
const
N1 :Integer = 8;
N2 :Integer = 8;
D1 :Real = 0.25;
D2 :Real = 0.25;
ALP :Real = 3.0E-02;
P :Real = 1.0;
label
_1;
begin
Result := '';  { результат функции }
for I:=1 to N1 do
 begin
  for J:=1 to N2 do
   begin
    T1 := D1*(-N1/2+I-1);
    T2 := D2*(-N2/2+J-1);
    GR[(I-1)+(J-1)*8] := 1.57079632679*Exp(-(T1*T1+T2*T2)/2.0);
_1:
    GI[(I-1)+(J-1)*8] := 0.0;
   end;
 end;
L := 1;
EC12C(F,GR,GI,N1,N2,D1,D2,ALP,P,L,H,FR,FI);
L := 2;
EC12C(FEC12C,GR,GI,N1,N2,D1,D2,ALP,P,L,H,FR,FI);
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%16.7f ',[H[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[]);
Result := Result + #$0D#$0A;
for _i:=0 to 63 do
 begin
  Result := Result + Format('%16.7f ',[FR[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[]);
Result := Result + #$0D#$0A;
for _i:=0 to 63 do
 begin
  Result := Result + Format('%16.7f ',[FI[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('tec12c',Result);  { вывод результатов в файл tec12c.res }
exit;
end;

end.

Unit fec12c_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

function fec12c(W1 :Real; W2 :Real): Complex;

implementation

function fec12c(W1 :Real; W2 :Real): Complex;
begin
{ Result - прототип имени функции F на FORTRANe }
Result := Cmplx(Exp(-0.25*(W1*W1+W2*W2))*3.14159265359,0.0);
exit;
end;

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.