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

Подпрограмма:  EF13R (модуль EF13R_p)

Назначение

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

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

Решение интегрального уравнения Фредгольма I рода: 

                     b
(1)    Au  ≡   ∫  K(y, x)  u(x) dx  =  f(y)  ,    y ∈ [c, d] ,
                    a

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

          d           b
(2)     ∫ P(y) [  ∫  K(y, x) u(x) dx - f(y) ]2 dy  + α Ω(u) ,
        c            a

где  P (y) > 0 - весовая функция,  α > 0 - параметр регуляризации,  Ω (u) - стабилизирующий функционал, определяемый параметрами  P1 > 0,  P2 ≥ 0,  P3 ≥ 0 : 

                     b
       Ω(u)  =  ∫  ( P1 u2 + P2 (du/dx)2 + P3(d2u/dx2)2 ) dx
                    a

Для дискретизации первого слагаемого в (2) используются квадратурные формулы трапеций по обеим переменным, а при аппроксимации второго слагаемого - формулы численного дифференцирования второго порядка для аппроксимации первой и второй производных функции  u (x). 

   Тогда дискретный аналог (2) имеет вид:
         M                N
(3)     ∑  Pi σi''  [  ∑  σj' Ki j uj - fi ]2  +  α u TC u .
         i=1               j=1

Здесь  Pi = P (yi),  uj = u (xj),  fi= f (yi),  Ki j = K (yi, xj);  a =x1 < x2 < ... < xN = b;  c = y1 < y2 < ... < yM = d - сетки узлов по  x и  y,  σj' и  σi'' - коэффициенты квадратурных формул по  x и  y соответственно; C - матрица квадратичной формы, аппроксимирующей стабилизирующий функционал  Ω (u);  u = (u1, u2, ..., uN) - искомое решение.

Пусть f = ( f1 √σ1'',  f2 √σ2'', ..., fM √σM'' ) и  A = (ai j),  ai j = Ki j Pi σj' √σi'',  1 ≤ i ≤ M,  1 ≤ j≤ N.

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

(4)     ( AT A + α C ) u  =  AT f ,
   где T означает операцию транспонирования.

При решении этой системы используется сведение задачи к более простой, канонической задаче с единичной матрицей  C и двухдиагональной матрицей  A согласно методу В.В.Воеводина [2]. Подпрограмма вычисляет значения пяти критериальных функций, которые могут быть использованы для выбора параметра регуляризации. A именно, в подпрограмме вычисляются приближенные: 

1)  Невязка
                      d           b
     ρ(α)  =  (  ∫ P(y) [  ∫  K(y, x) uα(x) dx - f(y) ]2 dy )1/2 ;
                     c           a
2)  Функция  "чувствительности"
     τ(α)  =  α Ω1/2 (duα / dα) ;

3)  Hоpма решения
      γ(α)  =  Ω1/2 (uα ) ;

4)  Значение регуляризирующего функционала
      φ(α)  =  ρ2(α) + αγ2(α) ;

5)  Функция отношения
        ψ(α)  =  ρ2(α) / ρ(α) ,
где
                        d           b
     ρ2(α)  =  {  ∫ P(y) [  ∫  K(y, x) ( α d uα(x)/dα - uα(x) ) dx - f(y) ]2 dy }1/2 .
                       c           a

Здесь  uα = uα (x) - функция, минимизирующая (2).

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

Описываемая подпрограмма реализует эти возможности и позволяет:

1)  Выполнить только приведение задачи к каноническому виду, не решая интегрального уравнения (это требует порядка  M * N2 операций).

При работе с задачей, приведенной к каноническому виду, подпрограмма позволяет:

2)  Находить приближенное решение интегрального уравнения и значения критериальных функций при заданном значении паpаметpа регуляризации (это требует порядка  N2 операций).
3)  При заданном значении  α находить только значения критериальных функций, не находя решения (это требует порядка N операций).
4)  Если необходимо решить уравнение (1) с другой правой частью, но с тем же ядром, то подпрограмма позволяет выполнить необходимое преобразование новой задачи к "канонической" за число операций порядка  M * N.

1. Тихонов A.H., O регуляризации некоppектно поставленных задач, ДАН CCCP, 1963, т.153, N 1.

2. Воеводин B.B., O методе регуляризации, ЖВМ и МФ, 1969, т.9, N 3.

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

procedure EF13R(var SX :Array of Real; var SY :Array of Real; N :Integer;
                M :Integer; var A :Array of Real; var U :Array of Real;
                var F :Array of Real; var P :Array of Real; P1 :Real;
                P2 :Real; P3 :Real; Q :Real; var H :Array of Real;
                var W :Array of Real; L :Integer);


   Параметры
SX - вещественный вектоp длины  N значений узлов сетки по переменной  X;  SX (I) = xI,  x1 = a,  xN = b;
SY - вещественный вектоp длины  M значений узлов сетки по переменной  Y;  SY (I) = YI,  Y1 = c,  YM = d;
N - число узлов сетки по переменной  X (тип: целый);
M - число узлов сетки по переменной  Y (тип: целый);
A - вещественный двумерный массив размера  M * N, содержащий значения ядра уравнения на заданных сетках:  A (I, J) = K (YI, XJ);
U - вещественный вектоp длины  N; в результате работы подпрограммы содержит значения приближенного решения в узлах сетки по  X:  U (I) = UI;
F - вещественный вектоp длины  M, содержащий значения на сетке по  Y правой части интегрального уравнения; в результате работы подпрограммы содержит значения правой части для "канонической" задачи;
P - вещественный вектоp длины  M, содержащий значения весовой функции в узлах сетки по  Y;
          P1 -
         P2  
         P3   		параметры, определяющие стабилизирующий функционал (тип: вещественный);
Q - заданный параметр регуляризации (тип: вещественный);
H - вещественный вектоp длины 50 ; в результате работы подпрограммы содержит вычисленные значения критериальных функций при заданном значении  Q:  H (1) = ρ (α),  H (2) = τ (α),  H (3) = γ (α),  H (4) = φ (α),  H (5) = ρ1 (α);
W - двумерный вещественный рабочий массив размеpа  N * 8;
L - параметр, определяющий режим использования подпрограммы (тип: целый);
L = 1 - выполняется сведение задачи к канонической;
L = 2 - вычисляется решение и значения критериальных функций;
L = 3 - вычисляются только значения критериальных функций;
L = 4 - выполняются дополнительные преобразования, приводящие задачу с новой правой частью к канонической: новая правая часть должна быть расположена на месте F. 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1) 

Подпрограмма EF13R обращается к вспомогательным подпрограммам с именами: EFFR, EFRS, EFUS, EFUF, EFSL, EFUD, EFRD.

  2)  При работе подпрограммы в режиме, выполняющем дополнительные преобразования, приводящие задачу с новой правой частью к "канонической" (p;4). Значения остальных перечисленных выше параметров (и рабочего массива  W) должны сохраняться, т.к. они используются в программе при всех значениях  L.

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

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

         1
         ∫ k(y, x) u(x) dx = F(y)
        0

с ядром  k(x, y) = y / (1 + y2 x2) и правой частью  f (y) = arctg y (точное решение  u ≡ 1) или  f (y) = ln (1 + y2) / 2y (точное решение  u ≡ x).

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие знначения критериальных функций при фиксированном параметре регуляризации, с весовой функцией  p ≡ 1 и использованием равномерных по  x и  y на отрезке [0,1] сеток из 11 точек. 

Unit TEF13R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EF13R_p;

function TEF13R: String;

implementation

function TEF13R: String;
var
I,J,_i :Integer;
U :Array [0..10] of Real;
F :Array [0..10] of Real;
A :Array [0..120] of Real;
H :Array [0..49] of Real;
W :Array [0..87] of Real;
const
P :Array [0..10] of Real = ( 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0 );
N :Integer = 11;
M :Integer = 11;
Q :Real = 1.0E-6;
X :Array [0..10] of Real = ( 0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0 );
Y :Array [0..10] of Real = ( 0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0 );
label
_1,_2;
begin
Result := '';
{ **** ТЕСТ ДЛЯ  S/R  EF13R В  БЧА НИВЦ МГУ ********* }
for I:=1 to M do
 begin
  F[I-1] := ArcTan(X[I-1]);
  for J:=1 to N do
   begin
_1:
    A[(I-1)+(J-1)*11] := Y[I-1]/(1.0+X[J-1]*X[J-1]*Y[I-1]*Y[I-1]);
   end;
 end;
EF13R(X,Y,N,M,A,U,F,P,1.0,0.0,0.0,Q,H,W,1);
EF13R(X,Y,N,M,A,U,F,P,1.0,0.0,0.0,Q,H,W,2);
Result := Result + Format('%s',[' PEШEHИE']);
Result := Result + #$0D#$0A;
for _i:=0 to 10 do
 begin
  Result := Result + Format('%20.16f ',[U[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 49 do
 begin
  Result := Result + Format('%20.16f ',[H[_i]]);
  if ( ((_i+1) mod 5)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
F[0] := 0.0;
for I:=2 to M do
 begin
_2:
  F[I-1] := Ln(1.0+Y[I-1]*Y[I-1])/(2.0*Y[I-1]);
 end;
Q := 1.E-6;
EF13R(X,Y,N,M,A,U,F,P,1.0,0.0,0.0,Q,H,W,4);
EF13R(X,Y,N,M,A,U,F,P,1.0,0.0,0.0,Q,H,W,2);
Result := Result + Format('%s',[' PEШEHИE']);
Result := Result + #$0D#$0A;
for _i:=0 to 10 do
 begin
  Result := Result + Format('%20.16f ',[U[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 49 do
 begin
  Result := Result + Format('%20.16f ',[H[_i]]);
  if ( ((_i+1) mod 5)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TEF13R',Result);  { вывод результатов в файл TEF13R.res }
exit;
end;

end.

Результаты:

       U  =  ( 0.99, 0.99, 0.99, 1.0, 1, 0, 1.0, 1.0, 1.0, 0.99, 0.98, 0.96 )
       H  =  ( 5.6*10-6,  5.7*10-4,  1.0,  1.0*10-6,  1.7*105 )