Подпрограмма:  EF01R (модуль EF01R_p)

Назначение

Приближенное решение интегрального уравнения Фредгольма I рода в классе ограниченных кусочно - монотонных функций.

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

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

            b
           ∫ k(x, y) u(y) dy  =  f(x)  ,  c ≤ x ≤ d ,
          a 

при априорно известной информации о качественном поведении искомого решения  u = u (y) осуществляется методом квазирешений В.К.Иванова путем сведения к решению задачи квадратичного программирования с линейными ограничениями

(1)     || K D u - f ||p2  =   inf   || K D v - f ||p2 ,
                                      v∈G 

где  K - известная матрица порядка  M * N с элементами  k_i j , которая в зависимости от ядра  k (x, y) и способа его аппроксимации может быть матрицей общего вида, либо симметричной, циркулянтной, теплицевой и т.п.,

f = ( f₁, ..., f_M ) - известный  M -  мерный вектор ,

u = ( u₁, ..., u_N ) - искомый  N - мерный вектор приближенного решения,

D - диагональная матрица коэффициентов квадратурной формулы  d_j ,  j = 1, ..., N,

                        M 
       || Z ||p2  =  ∑  pi zi2 ,      pi > 0  -  весовые коэффициенты.
                       i=1 

Линейные ограничения  G представляют собой одно из множеств  G₁,  G₂ или  G₁ ∩ G₂, где

 G1 = { v = (v1, ..., vN) :  aj ≤ vj ≤ bj ,  j = 1, ..., N } ,
 G2 = { v = (v1, ..., vN) :  mj (vj + 1 - vj) ≥ 0 ,  j = 1, ..., N - 1 } , 

a = (a₁, ..., a_N),  b = (b₁, ..., b_N) - заданные  N - мерные векторы нижних и верхних ограничений в  G₁; параметры  m_j определяются условиями кусочной монотонности и принимают значения 1, 0 или -1.

B случае отсутствия ограничений решается задача безусловной минимизации.

Численное решение задачи (1) основано на методе проекции сопряженных градиентов, который позволяет проводить итерационный процесс в конечномерном пространстве  W_2, ρ¹ с нормой

                                    N
 || v ||   =  ( || v || 2Q  +  ∑   ρj (vj - vj - 1)2 )1/2 ,
                                   j=2
                     N
 || v || Q  =  (  ∑   q j vj2 )1/2 ,
                    j=1 

q _j > 0 и  ρ_j ≥ 0 - заданные числа. Если априори известно, что искомое решение достаточно гладкое, то задание параметров  ρ_j > 0, связанных с конкретной нормировкой  W_2, ρ¹ ( при  ρ_j≡0 эта ноpма  v  совпадает с  || v || _Q ), позволяет увеличить точность и улучшить качество приближенного решения.

Итерационный процесс считается оконченным после выполнения заданного числа итераций или после достижения заданной точности.

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

B самой подпрограмме EF01R матрица  K не используется, поэтому способ ее представления не регламентируется.

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

procedure EF01R(var M :Integer; var N :Integer; var A :Array of Real;
                var F :Array of Real; NG :Integer; var AU :Array of Real;
                var BU :Array of Real; var MU :Array of Integer;
                var P :Array of Real; var Q :Array of Real;
                var D :Array of Real; var R :Array of Real; E :Real;
                ISP :Integer; var U :Array of Real; var FU :Real;
                var KN :Integer; var B :Array of Real; EF01R0 :Proc_EF01R);

   Параметры 

                N
      wi  =  ∑   ki j vj ,    i = 1, ..., M ,
               j=1 

                M
      wj  =  ∑   kj i vi ,    j = 1, ..., N .
               i=1 

   Версии:  нет 

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

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

                             xi+1    yj+1
      ki j  =  1/h1x       ∫       ∫     k (x, y)  dx dy ,   dj ≡ 1 .
                              xi      yj 

  
    _1:
    EF01R (M, N, A, F, NG, AU, BU, MU, P, Q, D, R, E, ISP,
           U, FU, KN, OP);
    if (KN<0)
     then goto _2
     else if (KN=0)
           then goto _1
           else goto _2;
    _2:
     ...  

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

Матрица  K имеет общий вид

Unit tef01r_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EF01R_p, EF01R1_p;

function tef01r: String;

implementation

function tef01r: String;
var
M,N,NG,ISP,_i,KN :Integer;
E,FU :Real;
B :Array [0..25] of Real;
const
A :Array [0..8] of Real = ( 3.0,1.0,5.0,1.0,-3.0,2.0,1.0,2.0,1.0 );
F :Array [0..2] of Real = ( 4.0,3.0,4.0 );
MU :Array [0..2] of Integer = ( 1,1,0 );
P :Array [0..2] of Real = ( 1.0,1.0,1.0 );
Q :Array [0..2] of Real = ( 1.0,1.0,1.0 );
D :Array [0..2] of Real = ( 1.0,1.0,1.0 );
R :Array [0..3] of Real = ( 0.0,0.0,0.0,0.0 );
U :Array [0..2] of Real = ( 0.0,0.0,0.0 );
AU :Array [0..2] of Real = ( -10.0,-10.0,-10.0 );
BU :Array [0..2] of Real = ( 10.0,10.0,10.0 );
label
_1,_2;
begin
Result := '';
M := 3;
N := 3;
NG := 2;
E := 1.E-7;
ISP := 20;
_1:
EF01R(M,N,A,F,NG,AU,BU,MU,P,Q,D,R,E,ISP,U,FU,KN,B,EF01R1);
Result := Result + Format('%s',['  U= ']);
Result := Result + #$0D#$0A;
for _i:=0 to 2 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 + Format('   %20.16f',[FU]) + #$0D#$0A;
if ( KN ) < 0
 then goto _2
 else if ( KN ) > 0
       then goto _2
       else goto _1;
_2:
UtRes('tef01r',Result);  { вывод результатов в файл tef01r.res }
exit;
end;

end.


Unit ef01r1_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc;

procedure ef01r1(M :Integer; N :Integer; var A :Array of Real;
                 var V :Array of Real; var W :Array of Real; NOP :Integer);

implementation

procedure ef01r1(M :Integer; N :Integer; var A :Array of Real;
                 var V :Array of Real; var W :Array of Real; NOP :Integer);
var
I,J :Integer;
label
_1,_2,_3,_4,_5;
begin
case NOP of
 1: goto _1;
 2: goto _3;
end;
_1:
for I:=1 to M do
 begin
  W[I-1] := 0.0;
  for J:=1 to N do
   begin
    W[I-1] := W[I-1]+A[(I-1)+(J-1)*M]*V[J-1];
_2:
   end;
 end;
goto _5;
_3:
for J:=1 to N do
 begin
  W[J-1] := 0.0;
  for I:=1 to M do
   begin
    W[J-1] := W[J-1]+A[(I-1)+(J-1)*M]*V[I-1];
_4:
   end;
 end;
_5:
end;

end.

Результат:    U  =  (-1., 2., 5.)