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

Подпрограмма:  EC02C (модуль EC02C_p)

Назначение

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

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

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

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

procedure EC02C(var AR :Array of Real; var AI :Array of Real;
                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);


   Параметры
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 - вычисление двумерного дискретного или обратного дискретного преобразования Фурье комплексной матрицы размера  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.

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

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

   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. 

Unit TEC02C_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, EC02C_p;

function TEC02C: String;

implementation

function TEC02C: String;
var
I,J,L,_i :Integer;
T1,T2,T :Real;
AR :Array [0..63] of Real;
AI :Array [0..63] of 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;
XR :Array [0..63] of Real;
XI :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);
    T := T1*T1+T2*T2;
    AR[(I-1)+(J-1)*8] := Exp(-T);
    AI[(I-1)+(J-1)*8] := 0.0;
    GR[(I-1)+(J-1)*8] := 1.57079632679*Exp(-T/2.0);
    GI[(I-1)+(J-1)*8] := 0.0;
    XR[(I-1)+(J-1)*8] := Exp(-T);
_1:
    XI[(I-1)+(J-1)*8] := 0.0;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    ЯДPO' + #$0D#$0A]);
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[AR[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[AI[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    ПРАВАЯ ЧACTЬ' + #$0D#$0A]);
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[GR[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[GI[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    ТОЧНОЕ PEШEHИE' + #$0D#$0A]);
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[XR[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[XI[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
L := 1;
EC02C(AR,AI,GR,GI,N1,N2,D1,D2,ALP,P,L,H,FR,FI);
L := 2;
EC02C(AR,AI,GR,GI,N1,N2,D1,D2,ALP,P,L,H,FR,FI);
Result := Result + Format('%s',['                    КРИТЕРИАЛЬНЫЕ ФYHKЦИИ' + #$0D#$0A]);
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 + #$0D#$0A;
Result := Result + Format('%s',['                    PEГYЛЯРИЗОВАННОЕ PEШEHИE' + #$0D#$0A]);
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[FR[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for I:=1 to 8 do
 begin
  for J:=1 to 8 do
   begin
    Result := Result + Format('%16.7f ',[FI[(I-1)+(J-1)*8]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
UtRes('TEC02C',Result);  { вывод результатов в файл TEC02C.res }
exit;
end;

end.

Результаты:

       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.