Подпрограмма:  EF02R (модуль EF02R_p)

Назначение

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

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

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

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

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

                   d      b   
     min        ∫   (  ∫  k(x, y)  u(y) dy - f(x) )2 dx
    v, u (a)    c      a 

на множестве кусочно - монотонных функций  v (y) = u' (y), определенных на отрезке [a, b]

                                 y
          u(y)  =  u(a) + ∫ v(t) dt ,   a ≤ y ≤ b ,
                               a
 и удовлетворяющих ограничениям
        v1(y) ≤ v(y) ≤ v2(y) ,   a ≤ y ≤ b , 

v₁ (y) и  v₂ (y) - заданные на отрезке [a, b] функции. Значение  u (a) может считаться как известным, так и подлежащим определению.

Численная дискретизация этой вариационной задачи после аппроксимации соответствующих интегралов с помощью квадратурных формул на сетках  ω_x = {x_i :  c = x₁ < ... < x_i < ... < x_M = d} и ω_y = {y_j :  a = y₁ < ... < y_j < ... < y_N = b} приводит к задаче квадратичного программирования

                        M            N
  (1)      min      ∑  pi  (   ∑  dj ki j uj  -  fi  )2
            v, u1    i=1         j=1 

при условии, что вектоp  v = (v₁, ..., v_N - 1) удовлетворяет ограничениям

                               j -1
  (2)     uj  =  u1  +   ∑   he ve ,     j = 2, ..., N
                               e=1
        
  (3)     v1 j ≤ vj ≤ v2 j ,     j = 1, ..., N-1

  (4)     mj ( vj + 1 - vj ) ≥ 0 ,     j = 1, ..., N-2 , 

где  k_i j = k (x_i, y_j),  f_i = f (x_i),  d_j > 0 - коэффициенты квадратурной формулы трапеции на сетке  ω_y ,  p_i > 0 - весовые коэффициенты,  h_j = y_j + 1 - y_j - шаг сетки  ω_y ,  v_1j = v₁ (y_j),  v_2j = v₂ (y_j), параметры  m_j определяются условиями кусочной монотонности функции  v = u' и принимают значения 1, 0 или -1.

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

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

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

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

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

procedure EF02R(var M :Integer; N :Integer; var A :Array of Real;
                var F :Array of Real; var NG :Integer;
                var V1 :Array of Real; var V2 :Array of Real;
                var MU :Array of Integer; var H :Array of Real;
                var P :Array of Real; var Q :Array of Real;
                var R :Array of Real; var E :Real; var ISP :Integer;
                LU :Integer; var U :Array of Real; var FU :Real;
                var KN :Integer; var B :Array of Real);

   Параметры 

   Версии:  нет 

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

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

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

Unit tef02r_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EF02R_p;

function tef02r: String;

implementation

function tef02r: String;
var
I,J,_i,KN :Integer;
E,SF,Y,S,FU :Real;
A :Array [0..120] of Real;
F :Array [0..10] of Real;
U :Array [0..10] of Real;
B :Array [0..89] of Real;
const
V1 :Array [0..10] of Real = ( 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 );
V2 :Array [0..10] of Real = ( 1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,0.0 );
MU :Array [0..10] of Integer = ( 1,1,1,1,1,-1,-1,-1,-1,0,0 );
H :Array [0..10] of Real = ( 0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.0 );
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 );
Q :Array [0..10] of Real = ( 0.05,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.05 );
R :Array [0..11] of Real = ( 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 
);
M :Integer = 11;
N :Integer = 11;
NG :Integer = 3;
LU :Integer = 1;
ISP :Integer = 20;
label
_1,_2,_3,_4,_5;
begin
Result := '';
E := 1.E-7;
for I:=1 to M do
 begin
  for J:=1 to N do
   begin
    if ( I-J ) < 0
     then goto _2
     else if ( I-J ) > 0
           then goto _2
           else goto _1;
_1:
    A[(I-1)+(J-1)*11] := 1.0/Q[I-1];
    goto _3;
_2:
    A[(I-1)+(J-1)*11] := 0.0;
_3:
   end;
 end;
U[0] := -11.0/24.0;
for J:=2 to N do
 begin
_4:
  U[J-1] := 0.0;
 end;
SF := 0.025;
Y := -0.5;
for I:=1 to M do
 begin
  S := SF*Sin(123.456*I+789.0);
  F[I-1] := Y-(IntPower(Y,3))/3.0+S;
_5:
  Y := Y+H[I-1];
 end;
EF02R(M,N,A,F,NG,V1,V2,MU,H,P,Q,R,E,ISP,LU,U,FU,KN,B);
Result := Result + Format('%s',[' FU=']);
Result := Result + Format('%20.16f',[FU]) + #$0D#$0A;
Result := Result + Format('%s',['  KN=']);
Result := Result + Format('%3d',[KN]) + #$0D#$0A;
Result := Result + Format('%s',[' U=']);
Result := Result + #$0D#$0A;
for _i:=0 to 10 do
 begin
  Result := Result + Format('%20.16f',[U[_i]]);
  if ( ((_i+1) mod 1)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('tef02r',Result);  { вывод результатов в файл tef02r.res }
exit;
end;

end.


Результат:

       U  =  ( -4.58E-1,   -4.11E-1,   -3.16E-1,   -2.21E-1,   -1.26E-1, 
                  -3.13E-2,    7.87E-2,    1.73E-1,    2.68E-1,    3.64E-1,   4.58E-1 )