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

Подпрограмма:  EC00R (модуль EC00R_p)

Назначение

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

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

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

                   +∞
(1)    Ax  ≡   ∫  k(t - τ)  x(τ) dτ  =  y(t)  ,    -∞ < t < +∞ ,
                  -∞

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

(2)     Mα [x, y]  =  || Ax - y ||2 + α Ω[x] ,
 где :
                                        +∞
       || Ax - y ||2  =  1/2π    ∫   | K(λ) X(λ) - Y(λ) |2 dλ  -  квадрат невязки,
                                       -∞
      α > 0 - параметр регуляризации,
                             +∞
       Ω[x]  =  1/2π  ∫  | λ |2p | X(λ) |2 dλ  -  стабилизирующий функционал
                            -∞
 порядка  p ≥ 0 (p не предполагается целым);

    K (λ), X (λ), Y (λ) - преобразования Фурье ядра  k, решения  x и приближенной правой части  y. 

 Функционал (2) минимизирует функция
                                +∞
(3)     xα(t)  =  1/2π  ∫  K*(λ) Y(λ) / ( | K(λ) |2 + α λ2p ) exp (iλt) dt ,
                              -∞
          -∞ < t < +∞ ,
 где  *  - знак комплексного сопряжения.

Считается, что  k (t) → 0,  x (t) → 0,  y (t) → 0 при  t → ±∞ и функции  k (t),  y (t) заданы на конечном отрезке [-Т/2, Т/2] в узлах равномерной сетки по  t : 

     ts  =  s Δt ,   s =  -N/2, -N/2 + 1, ... , -1, 0, 1, ... , N/2 - 1 ,

где  N - четное число узлов,  Δ (t) = T/N - шаг сетки.

Задача состоит в решении дискретного уравнения типа свертки с ядром  ks = k (ts) и приближенной правой частью  ys = = y (ts) методом регуляризации  p - ого порядка в спектральной области  λ. Численный алгоритм основан на применении Дискретного Преобразования Фурье временных рядов [2].

Приведем задачу к "каноническому" виду, вычислив ДПФ функций  ks и  ys : 

                         N/2-1
             Km   =    ∑      ks exp(-i λm ts) ,
                        S = -N/2
                         N/2-1
             Ym   =    ∑      ys exp(-i λm ts) ,
                       S = -N/2

где  λm = m Δλ ,  m =  - N/2, - N/2 + 1 , ..., - 1, 0, 1, ..., N/2 - 1 ,  Δλ = 2π/(N Δt) ,  i = √-1.

Аналогично определяются  Xm,  Xmα  -  ДПФ рядов  xs = x (ts)  и  xsα =xα (ts).

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

                                 N/2-1
(4)   xsα  = 1/ NΔt      ∑      KTmYm /( | Km |2 + α/Δt2 λm2p ) exp(2π i m s / N) ,
                               m= -N/2
        s  =  -N/2, -N/2 + 1, ..., N/2 - 1 .

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

1) невязка    ρ(α)  =  || Axα - y ||  ,         

2) ноpма решения     γ(α)  =  Ω1/2 [ xα ]  ,    

3) значение    регуляризующего    функционала 
    φ(α)  =  [ ρ2(α) + αγ2(α) ]1/2  ,
4) функция "чувствительности"
   (для  выбора  квазиоптимального значения  α)
    τ(α)  =  α Ω1/2 [ dxα/dα ] ,
   где
                           +∞
    Ω[x]  =  1/2π  ∫   | λ | | X(λ) |2 dλ  ≡  Ω[x]
                         -∞
    при  p = 1/2.

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

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

1.  Задача приводится к "каноническому" виду; этот режим позволяет избежать повторного вычисления  Km,  Ym, когда происходит повторное обращение к подпрограмме с тем же ядром и той же правой частью, например, при выборе  α  (L = 1).
2.  При условии, что задача приведена к "каноническому" виду, строится решение  xsα и вычисляются значения крикритериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) при заданном значении  α (L = 2).
3.  При условии, что задача приведена к "каноническому" виду, вычисляются только значения критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) при при заданном значении  α без построения решения  xsα (L = 3).
4.  Вычисляется только ДПФ правой части  ys; этот режим необходим для приведения задачи к "каноническому" виду, когда повторно решается уравнение с тем же ядром, но с другой правой частью (L = 4).

Для оптимального вычисления прямого и обратного ДПФ применяся алгоритм быстрого преобразования Фурье [2], за счет которого реализованный в подпрограмме алгоритм численного решения задачи является эффективным в смысле быстродействия, экономии памяти ЭВМ и уменьшения ошибок округления, при этом время вычислений пропорционально  N log2N, а объем памяти пропорционален  N.

Полное описание изложенного алгоритма содержится в статье [3].

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

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

procedure EC00R(N :Integer; DT :Real; var A :Array of Real;
                var Y :Array of Real; ALP :Real; P :Real; L :Integer;
                var H :Array of Real; var X :Array of Real;
                var Z :Array of Real);


   Параметры
N - количество заданных значений ядра и правой части уравнения,  N ≥ 4 - целая степень числа два (тип: целый);
DT - заданный шаг сетки по переменным  t,  τ,  DT > 0 (тип: вещественный);
A - вещественный вектоp длины  N, содержащий заданные значения ядра уравнения на сетке:  A (I) = kI-N/2-1,  I = 1, 2, ..., N;
Y - вещественный вектоp длины  N, содержащий заданные значения правой части уравнения на сетке:  Y (I) = yI- N/2- 1,  I = 1, 2, ..., N;
ALP - заданное значение параметра регуляризации  α,  ALP ≥ 0 (тип: вещественный);
P - заданное значение порядка регуляризации  p,  P ≥ 0 (тип: вещественный);
L - параметр, определяющий режим использования подпрограммы (тип: целый):
L = 1 - выполняется приведение задачи к "каноническому" виду,
L = 2 - строится решение и вычисляются значения критериальных функций,
L = 3 - вычисляются только значения критериальных функций;
L = 4 - выполняется приведение задачи с новой правой частью с тем же ядром к "каноническому" виду;
H - вещественный вектоp длины 4, содержащий вычисленные значения критериальных функций при заданном значении  ALP:  H (1) = ρ(α),  H (2) = γ(α),  H (3) = φ(α),  H (4) = τ(α);
X - вещественный вектоp длины  N, содержащий вычисленные значения регуляризованного решения на сетке:  X (I) = xαI- N/2- 1,  I = 1, 2, ..., N;
Z - вещественный вектоp длины  N, используемый в подпрограмме как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы
FTF1C - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

B результате работы подпрограммы при  L = 1 в массиве  A содержится преобразование Фурье ядра  Km, а в массиве  Y - преобразование Фурье правой части  Ym; на выходе подпрограммы при  L = 3 в массиве  X содержится преобразование Фурье регуляризованного решения  Xmα.

C учетом симметрии вещественной части и антисимметрии мнимой части преобразований Фурье вещественных рядов относительно центра  m = N/2, хранятся только половины значений  Km,  Ym,  Xmα,  m = 0, 1, ..., N/2.

Например, в массиве  A спектр располагается следующим образом: 
     A(I) = Re KI-1 ,         I = 1, 2, ..., N/2 + 1 ,
     A(I) = Im KI-N/2-1 ,   I = N/2 + 2, N/2 + 3, ..., N ;

при этом учитывая также, что  Im Km = 0,  m = 0, N/2.

  2. 

При первом обращении к подпрограмме (L = 1) необходимо задать параметры  N, DT, A, Y, определяющие уpавнение. B дальнейшем (L = 2, 3) можно менять только параметры  ALP и  P, характеризующие регуляризатор, сохраняя значения  N, DT, A, Y. При  L = 4 значения параметров  N, DT, а также должны сохраняться.

  3.  Массивы  A, Y, X используются при всех значениях  L, а рабочий массив  Z - только в одном режиме при  L = 2.

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

Рассматривается решение интегрального уравнения типа свертки (1) с ядром  k (t) = exp (- t2) и правой частью  y (t) = (π/2) 1/2 exp (- t2/2) ( точное решение  x (t) = exp (- t2) ) или y1 (t) = (π/3) 1/2 exp (- 2/3 t2) ( точное решение  x1 (t) = exp (- 2t2) ).

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие значения критериальных функций при фиксированном параметре регуляризации  α = 10- 2 с использованием регуляризатора порядка  p = 1 и равномерной на отрезке [ - 1, 1 ] сетки из  N = 8 точек. 

Unit TEC00R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EC00R_p;

function TEC00R: String;

implementation

function TEC00R: String;
var
I,_i,L :Integer;
T :Real;
A :Array [0..7] of Real;
Y :Array [0..7] of Real;
H :Array [0..3] of Real;
X :Array [0..7] of Real;
Z :Array [0..7] of Real;
Y1 :Array [0..7] of Real;
ХТ :Array [0..7] of Real;
XT1 :Array [0..7] of Real;
const
N :Integer = 8;
DT :Real = 0.25;
ALP :Real = 1.0E-2;
P :Real = 1.0;
PI :Real = 3.14159265359;
label
_1;
begin
Result := '';
for I:=1 to N do
 begin
  T := DT*(-N/2+I-1);
  A[I-1] := Exp(-T*T);
  XT[I-1] := Exp(-T*T);
  Y[I-1] := Sqrt(PI/2.0)*Exp(-T*T/2.0);
  XT1[I-1] := Exp(-2.0*T*T);
_1:
  Y1[I-1] := Sqrt(PI/3.0)*Exp(-T*T*2.0/3.0);
 end;
Result := Result + Format('%s',['                    ЯДPO']);
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[A[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    ТОЧНОЕ PEШEHИE']);
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[XT[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    ПРАВАЯ ЧACTЬ']);
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[Y[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
L := 1;
EC00R(N,DT,A,Y,ALP,P,L,H,X,Z);
L := 2;
EC00R(N,DT,A,Y,ALP,P,L,H,X,Z);
Result := Result + Format('%s',['                    КРИТЕРИАЛЬНЫЕ ФYHKЦИИ']);
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%20.16f ',[H[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    PEГYЛЯРИЗОВАННОЕ PEШEHИE']);
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[X[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    ТОЧНОЕ PEШEHИE']);
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[XT1[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    ПРАВАЯ ЧACTЬ']);
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[Y1[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
L := 4;
EC00R(N,DT,A,Y1,ALP,P,L,H,X,Z);
L := 2;
EC00R(N,DT,A,Y1,ALP,P,L,H,X,Z);
Result := Result + Format('%s',['                    КРИТЕРИАЛЬНЫЕ ФYHKЦИИ']);
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%20.16f ',[H[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['                    PEГYЛЯРИЗОВАННОЕ PEШEHИE']);
Result := Result + #$0D#$0A;
for _i:=0 to 7 do
 begin
  Result := Result + Format('%20.16f ',[X[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TEC00R',Result);  { вывод результатов в файл TEC00R.res }
exit;
end;

end.

Результаты:

       H  = ( 0.13239,  1.08822,  0.17138,  0.33298 )
       X  = ( 0.378,  0.475,  0.713,  0.963,  1.072,  0.963,  0.713,  0.475 )