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

Подпрограмма:  EF12R (модуль EF12R_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

Значение параметра регуляризации  α определяется из условия: 

          ρ(α)  =  e || f ||p  ,
   где
                        d           b
       ρ(α)  =  (  ∫ p(y) [  ∫ k(y, x) uα(x) dx - f(y) ]2 dy )1/2   -
                      c            a

невязка на регуляризированном решении  uα (x),  e - заданный относительный уровень невязки, 

                         d        
       || f ||p  =  (  ∫  p(y) f 2(y) dy )1/2   -   норма правой части уравнения (1).
                       c

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

  Тогда дискретный аналог (2) имеет вид:
         M                N
(3)     ∑  pi σi''  [  ∑  σj' ki j uj - fi ]2  +  α u T C 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 означает операцию транспонирования.

При этом параметр  α выбирается из дискретного аналога условия невязки: 

(5)     ρ(α) ≈ || A uα - f ||  =  e || f || ,

где  uα - решение (4).

При решении системы (4) используется сведение задачи к более простой "канонической" задаче с единичной матрицей  С и двухдиагональной матрицей  A согласно методу В.В.Воеводина [2].

Для решения уравнения (5) используется метод касательных Ньютона в соответствии с вычислительной схемой, предложенной В.А.Морозовым [3].

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

1) невязка  ρ(α);
2) функция "чувствительности"
    τ(α) ≡ α Ω1/2 (duα / dα) ;
3) Ω - ноpма регуляризованного решения
    γ(α) ≈ Ω1/2 (uα) ;
4) значение   регуляризирующего   функционала 
    φ(α) = ρ2(α) + αγ2(α) ;
5) найденное значение  α как решение уравнения (5);
6) число итераций при решении уравнения (5);
7) достигнутый относительный уровень невязки
    e = ρ(α) / || f || .

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

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

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

II) Pешение того же уравнения, но с другим относительным уровнем невязки (это требует порядка  N2 операций).

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

1.  Тихонов A.H., O регуляризации некоppектно поставленных задач, ДАН CCCP, 1963, т.153, N 1.
2.  Воеводин B.B., O методе регуляризации, ЖВМ и МФ, 1969, т.9, N 3.
3.  Морозов B.A. O принципе невязки при решении несовместных уравнений методом регуляризации А.Н.Тихонова, ЖВМ и МФ, 1973, т.13, N 5.

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

procedure EF12R(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; EPS :Real; var HX :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   		параметры, определяющие стабилизирующий функционал (тип: вещественный);
EPS - заданный относительный уровень невязки (тип: вещественный);
H - вещественный вектоp длины 7. В результате работы подпрограммы содержит вычисленные характеристики решения:  H (1) = ρ (α),  H (2) = τ (α),  H (3) = γ (α),  H (4) = φ (α),  H (5) = α,  H (6) = число итераций,  H (7) = достигнутый относительный уpовень невязки;
W - двумерный вещественный рабочий массив размеpа  N * 8;
L - параметр, определяющий режим использования подпрограммы (тип: целый);
L = 1 - решение интегрального уравнения первый раз (режим I);
L = 2 - решение интегрального уравнения с новым значением относительного уровня невязки  EPS (режим II);
L = 3 - решение интегрального уравнения с новыми правой частью  F и относительным уpовнем невязки  EPS (режим III). 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1) 

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

  2) 

При первом обращении к подпрограмме (L = 1) необходимо задать параметры, определяющие уравнение:  SX, SY, N, M, A, F, P, P1, P2, P3, EPS, L.
B дальнейшем можно менять только значения  EPS (при  L = 2) и  F, EPS (при  L = 3).
Значения остальных перечисленных выше параметров (и рабочего массива  W) должны сохраняться, т.к. они используются в программе при всех значениях  L.

  3) 

Если требуемое значение  EPS не может быть достигнуто (при нахождении параметра регуляризации  α методом Ньютона значение  α достигает наименьшего положительного числа, представимого на ЭВМ, или при уменьшении  α невязка начинает возрастать), то вычисляется наилучшее возможное решение (с минимальной на вычисленных значениях  α невязкой или при наименьшем возможном значении  α > 0).

  4)  Значение относительного уровня невязки  EPS должно быть неотрицательно:  EPS ≥ 0. Если  EPS ≥ 1, то начальное приближение  U (x) = 0 является решением задачи.

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

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

           1
           ∫ k(y, x) - u(x) dx  =  f(y)
          c

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

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие значения критериальных функций при заданных значениях относительного уровня невязки  EPS с весовой функцией  P (y) = 1 и использованием равномерных по  X и  Y на отрезке [0, 1] сеток из 11 точек. 

Unit tef12r_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EF12R_p;

function tef12r: String;

implementation

function tef12r: 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..6] 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;
EPS :Real = 0.001;
SX :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 );
SY :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 := '';
for I:=1 to M do
 begin
  F[I-1] := ArcTan(SY[I-1]);
  for J:=1 to N do
   begin
_1:
    A[(I-1)+(J-1)*11] := SY[I-1]/(1.0+SX[J-1]*SX[J-1]*SY[I-1]*SY[I-1]);
   end;
 end;
EF12R(SX,SY,N,M,A,U,F,P,1.0,0.0,0.0,EPS,H,W,1);
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 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
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;
F[0] := 0.0;
for I:=2 to N do
 begin
_2:
  F[I-1] := Ln(1.0+SY[I-1]*SY[I-1])/(2.0*SY[I-1]);
 end;
EF12R(SX,SY,N,M,A,U,F,P,1.0,0.0,0.0,EPS,H,W,3);
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 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
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;
UtRes('tef12r',Result);  { вывод результатов в файл tef12r.res }
exit;
end;

end.


Результаты:

 PEШEHИE
       U  =  ( 1.01, 1.01, 1.01, 1.01, 1.01, 1., 1., 0.99, 0.97, 0.96, 0.94 )
       H  =  ( 4.96*10-4, 1.46*10-2, 0.99, 1.07*10-4, 1.07*10-4, 1., 1.0*10-3 )


 PEШEHИE
       U  =  ( 0.11, 0.13, 0.18, 0.27, 0.37, 0.48, 0.60, 0.70, 0.80, 0.88, 0.94 )
       H  =  ( 2.3*10-4, 8.54*10-3, 0.57, 5.85*10-6, 1.78*10-5, 1., 1.0*10-3 )