Подпрограмма:  EC21R (модуль EC21R_p)

Назначение

Решение одномерного интегрального уравнения Винера - Хопфа первого рода с симметричным неотрицательно определенным ядром методом регуляризации первого порядка без выбора параметра регуляризации.

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

Приближенное решение интегрального уравнения Винера - Хопфа первого рода

                    ∞   
(1)    Au  ≡   ∫    K(x - z) u(t) dt  =  f(x)  ,
                    0
          x ∈ [0, +∞) ,    K(x - t)  =  K(t - x)  , 

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

             ∞   ∞
(2)        ∫     ∫   K(x - t) u(x) u(t) dx dt  -
           0     0
                           ∞                             ∞
                    -  2  ∫   f(x) u(x) dx  +  α ∫   u ' 2(x) dx  ,
                          0                              0
   где  α > 0  -  параметр регуляризации ,
      ∞
      ∫   u ' 2(x) dx  -  стабилизирующий функционал. 
     0 

Считается, что  K (t)→0,  u (t)→0,  f (t)→0 (достаточно быстро) при  t→∞, так что функции  K (t),  f (t) можно считать заданными на конечном отрезке [0, T].

Тогда вместо (2) будем иметь:

             T    T
(3)        ∫     ∫   K(x - t) u(x) u(t) dx dt  -
           0     0
                            T                             T
                    -  2  ∫   f(x) u(x) dx  +  α ∫   u ' 2(x) dx  ,
                          0                              0 

Для дискретизации первого и второго слагаемого (3) используется квадратурная формула трапеций на равномерной сетке  0 = x₁ < x₂ < ... < x_N = T ,  x_i + 1 - x_i = Δx , а при аппроксимации третьего слагаемого - разностные отношения:

     du/dx (xi)  =  ( u(xi+1) - u(xi) ) / Δx ,        i = 1, 2, ..., N-1 . 

B предположении  du/dx (x_N) = 0 дискретный аналог (3) имеет вид:

(4)    (D K D u, u - 2(D f, u) + α (C u, u).

   Здесь: 

D = diag {d₁₁, d₂₂, ..., d_NN},  причем  d₁₁ = d_NN = Δx /2 ,  d_i i = Δx ,  i = 2, 3, ..., N-1 ;

K - симметричная теплицева матрица, определенная элементами первой строки  K_1 i = K (x_i) ;

u = (u₁, u₂, ..., u_N) ,  u_i = u (x_i) ,

f = (f₁, f₂, ..., f_N) ,  f_i = f (x_i) ,  i = 1, 2, ..., N ,

C - трехдиагональная матрица, аппроксимирующая стабилизирующий оператор.

Задача минимизации (4) по  u эквивалентна решению системы линейных алгебраических уравнений:

(5)     D K D u + α C u  =  D f . 

При ее решении используется схема, предложенная в [2], в основе которой лежит существенное использование теплицевости матрицы  K и "квазитеплицевости" матрицы  C.

Алгоритм, описанный в [3], был адаптирован для решения систем уравнений с симметричной теплицевой матрицей. Эта адаптация дает почти двойной выигрыш по времени.

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

procedure EC21R(var A :Array of Real; var X :Array of Real;
                var BA :Array of Real; var B :Array of Real;
                var C :Array of Real; N :Integer; T :Real; J1 :Integer;
                ALFA :Real);

   Параметры 

   Версии:  нет 

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

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

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

Рассматривается решение интегрального уравнения (1) с ядром

    K(t) = exp( -| t | )
    и правой частью  F(t) = 2/5 exp(-t) ( cos(t) + 2 sin(t) )
    ( точное решение  u(t) = exp(-t) cos(t) ). 

Ниже приводится фрагмент программы, вычисляющей регуляризованные решения:

Unit TEC21R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EC21R_p;

function TEC21R: String;

implementation

function TEC21R: String;
var
N,J1,I :Integer;
ALFA,T,H :Real;
A :Array [0..199] of Real;
X :Array [0..199] of Real;
ВА :Array [0..199] of Real;
B :Array [0..199] of Real;
C :Array [0..199] of Real;
UTON :Array [0..199] of Real;
label
_4,_1;
begin
Result := '';
T := 10.0;
N := 200;
H := T/(N-1);
J1 := 0;
ALFA := 1.0/(IntPower(10.0,9));
for I:=1 to N do
 begin
  UTON[I-1] := Exp((1.0-I)*H)*Cos((I-1.0)*H);
_4:
 end;
Result := Result + Format('%20.16f ',[UTON[0]]) + #$0D#$0A;
Result := Result + #$0D#$0A;
I := 10;
WHILE ( I<=N ) do
 begin
  Result := Result + Format('%20.16f ',[UTON[I-1]]) + #$0D#$0A;
  inc(I,30);
 end;
Result := Result + #$0D#$0A;
for I:=1 to N do
 begin
  A[I-1] := Exp((1.0-I)*H);
  BA[I-1] := 2.0/5.0*Exp((1.0-I)*H)*(Cos((I-1.0)*H)+2.0*Sin((I-1.0)*H));
_1:
 end;
EC21R(A,X,BA,B,C,N,T,J1,ALFA);
Result := Result + Format('%20.16f ',[X[0]]) + #$0D#$0A;
Result := Result + #$0D#$0A;
I := 10;
WHILE ( I<=N ) do
 begin
  Result := Result + Format('%20.16f ',[X[I-1]]) + #$0D#$0A;
  inc(I,30);
 end;
Result := Result + #$0D#$0A;
UtRes('TEC21R',Result);  { вывод результатов в файл TEC21R.res }
exit;
end;

end.

Результаты:

       X(1)      =   0.9830 ,
       X(10)    =   0.5722 ,
       X(40)    =  -0.0534 ,
       X(70)    =  -0.0296 ,
       X(100)  =   0.0018 ,
       X(130)  =   0.0015 ,
       X(160)  =  -0.0000 ,
       X(190)  =  -0.0001 .