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

Подпрограмма:  EC01R (модуль EC01R_p)

Назначение

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

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

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

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

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

(2)     Mα [x, y]  =  || Ax - y ||2 + α Ω[x]   =  min 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.
Значение параметра  α определяется из условия 

(3)     ρ(α) = ε || y || ,   ρ(α) ≡ || Axα - y || ,

где  xα - регуляризированное решение задачи (2),
        ε - заданный относительный уровень невязки.

Считается, что  k (t) ,  x (t),  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 - шаг сетки.

Численный алгоритм для решения уравнения типа свертки с дискретным ядром и правой частью методом регуляризации основан на применении дискретного преобразования Фурье временных рядов [2].

B связи с этим выполняются преобразования, заключающиеся в вычислении ДПФ функций  ks = k (ts) и  ys = y (ts), приводящие исходную точку к "каноническому" виду.

Приближенные значения  xsα регуляризованных решений в узлах  ts определяются путем минимизации в спектральной области дискретной аппроксимации сглаживающего функционала и последующего выполнения обратного ДПФ. Для приближенного решения уравнения (3) используется метод касательных Ньютона в соответствии с вычислительной схемой [3].

Качество регуляризованного решения  xα при значении параметра  α, удовлетворяющем критерию (3), оценивается следующими числовыми характеристиками : 

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

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

3) значением регуляризирующего функционала 
    φ(α)  =  [ ρ2(α) + αγ2(α) ]1/2  ,

4) функцией "чувствительности"    τ(α)  =  α Ω1/2 [ dxα/dα ] ,
    где  Ω [x] ≡ Ω [x]  при  p = 1/2
    (для вычисления "квазиоптимальных" значений
    параметра регуляризации),

5) найденным значением  α как решением уравнения (3),

6) числом итераций при решении уравнения (3),

7) практически достигнутым относительным уpовнем невязки  ε0.

Вычисление всех критериальных функций  ρ (α),  γ (α),  φ (α),  τ (α) происходит в спектральной области  λ и требует порядка  N операций.

Возможны следующие режимы вычислений (в зависимости от значения параметра  L).

I.   Решение интегрального уравнения с приведением задачи к "каноническому" виду (L = 1).
II.  Решение интегрального уравнения при условии, что задача приведена к "каноническому" виду (L = 2).
III. Решение интегрального уравнения при условии, что ДПФ ядра вычислено (L = 3).

Эти режимы целесообразно использовать, например, при повторном решении того же интегрального уравнения (L = 2) или при решении уравнения с тем же ядром, но с другой правой частью (L = 3), а также в случае, когда известно аналитическое преобразование Фурье ядра  K (λ)  (L = 3).

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

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

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

procedure EC01R(N :Integer; DT :Real; var A :Array of Real;
                var Y :Array of Real; P :Real; EPS :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;
P - заданное значение порядка регуляризации  p,  P ≥ 0 (тип: вещественный);
EPS - заданный относительный уровень невязки  ε ,  0 < EPS < 1 (тип: вещественный);
L - параметр, определяющий режим использования подпрограммы (тип: целый);
L = 1 - решение интегрального уравнения в режиме I,
L = 2 - решение интегрального уравнения в режиме II,
L = 3 - решение интегрального уравнения в режиме III;
H - вещественный вектоp длины 7, содержащий вычисленные характеристики решения при найденном значении  α: H (1) = ρ(α),  H (2) = γ(α),  H (3) = φ(α),  H (4) = τ(α); H (5) = α,  H (6) - число итераций,  H (7) - достигнутый относительный уровень невязки  ε0;
X - вещественный вектоp длины  N, содержащий вычисленные значения регуляризованного решения на сетке:  X (I) = xI- N/2- 1,  I = 1, 2, ..., N;
Z - вещественный вектоp длины  N, используемый в подпрограмме как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы
FTF1C - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины, равной степени двух, методом быстрого преобразования Фурье. 
   Замечания по использованию
  1. 

В результате работы подпрограммы при  L = 1 в массиве  A содержится ДПФ ядра  Km, а в массиве  Y - ДПФ правой части  Ym ("канонический" вид задачи).

C учетом симметрии вещественной части и антисимметрии мнимой части преобразований Фурье вещественных рядов относительно центра  m = N/2, хранятся только половины значений  K, Y,  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, EPS.
При повторном решении того же интегрального уpавнения (L = 2) можно менять только параметры  P, EPS, а при решении уравнения с тем же ядром, но с другой правой частью (L = 3) - только параметры  Y, P, EPS, сохраняя значения всех остальных параметров.

  3. 

Если относительный уровень невязки  ε таков, что  ε ≥ ρ (∞)/|| y || ,   где 

ρ(∞)  =  lim  ρ(α) ,
           α → ∞

и величины  ρ (∞),  || y || вычисляются приближенно, то в подпрограмме строится приближение к предельному решению 

x∞  =  lim  xα  ≡  0
        α → ∞

с характеристиками:  H (1) = ρ (∞),  H (2) = 0,  H (3) = ρ (∞),  H (4) = 0,  H (5) = 1018,  H (6) = 0,  H (7) = ε0.

  4.  Если относительный уровень невязки  ε слишком мал и не может быть достигнут, то в подпрограмме вычисляется (если это возможно) нерегуляризованное pешение  xsα при  α = 0.

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

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

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

Unit TEC01R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EC01R_p;

function TEC01R: String;

implementation

function TEC01R: String;
var
I,_i,L :Integer;
T,EPS :Real;
A :Array [0..7] of Real;
Y :Array [0..7] of Real;
H :Array [0..6] 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;
P :Real = 1.0;
PI :Real = 3.14159265;
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;
EPS := 0.08;
Result := Result + Format('%s',['          EPS= ']);
Result := Result + Format('%20.16f',[EPS]) + #$0D#$0A;
EC01R(N,DT,A,Y,P,EPS,L,H,X,Z);
Result := Result + Format('%s',['                    КРИТЕРИАЛЬНЫЕ ФYHKЦИИ']);
Result := Result + #$0D#$0A;
for _i:=0 to 6 do
 begin
  Result := Result + Format('%20.16f ',[H[_i]]);
  if ( ((_i+1) mod 2)=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 := 3;
EPS := 0.085;
Result := Result + Format('%s',['          EPS= ']);
Result := Result + Format('%20.16f',[EPS]) + #$0D#$0A;
EC01R(N,DT,A,Y1,P,EPS,L,H,X,Z);
Result := Result + Format('%s',['                    КРИТЕРИАЛЬНЫЕ ФYHKЦИИ']);
Result := Result + #$0D#$0A;
for _i:=0 to 6 do
 begin
  Result := Result + Format('%20.16f ',[H[_i]]);
  if ( ((_i+1) mod 2)=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('TEC01R',Result);  { вывод результатов в файл TEC01R.res }
exit;
end;

end.

Результаты:

       H  =  ( 0.122,  1.213,  0.164,  0.335,  0.008,  3,  0.080 )
       X  =  ( 0.339,  0.447,  0.712,  0.991,  1.113,  0.991,  0.712,  0.447 )

       L = 3
       EPS = 0.085
       CALL  EC01R (N, DT, A, Y1, P, EPS, L, H, X, Z)

Результаты:

       H  =  ( 0.102,  1.492,  0.149,  0.325,  0.005,  3,  0.085 )
       X  =  ( 0.094,  0.225,  0.549,  0.893,  1.046,  0.893,  0.549,  0.225 )