Подпрограмма:  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 ∈ [- 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), которые могут быть использованы для выбора параметра регуляризации  α.
Вычисление всех критериальных функций происходит через компоненты ДПФ, что требует порядка  N₁ N₂ операций и значительно облегчает задачу выбора  α.
Выполнение обратного ДПФ в формуле (3) целесообразно производить лишь для выбранных значений  α.
За счет применения БПФ полное время численного решения задачи пропорционально  N₁ N₂ (log₂ N₁ + log₂ N₂), а объем памяти ЭВМ пропорционален  N₁ N₂.

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

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

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

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

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

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

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);

   Параметры 

      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; 

   Версии:  нет 

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

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

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

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