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

Подпрограмма:  EF04R (модуль EF04R_p)

Назначение

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

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

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

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

с гладким ядром и правой частью  f (x), заданной с погрешностью  δ (x), при априорно известной информации об ограниченности первой призводной исходного решения  u (y) и об участках знакопостоянства его кривизны сводится методом поточной невязки с дескриптивной регуляризацией [1,2] к отысканию общей точки  u (y) множеств  Uδ , Uν , Uμ : 

(2)     Uδ = { u:  | Au - f(x) | ≤ δ(x) }  ,
(3)     Uν = { u:  v1(y) ≤ u'(y) ≤ v2(y)   ∀y ∈ [a, b] }  ,
(4)     Uμ = { u:  μ(y) u''(y) ≥ 0   ∀y ∈ [a, b] ,  μ(y) = sign u''(y) }

т.е. к выделению из множества  Uδ формальных решений уравнения (1) приближенного решения  u (y), удовлетворяющего условиям (3), (4). v1 (y) и  v2 (y) - заданные на отрезке [a, b] функции.

Дискретный аналог задачи (2) - (4) состоит в определении вектоpа  u = (u1, ..., uN)∈EN, удовлетворяющего системе линейных неpавенств: 

(5)     | (K Du)i - fi | ≤ δi ,      i = 1, 2, ..., M  ,
(6)     v1 j ≤ (uj+1 - uj) / hj ≤ v2 j ,      j = 1, 2, ..., N-1  ,
(7)     μj ( (uj+1 - uj) / hj - (uj - uj-1) / hj-1 ) ≥ 0 ,      j = 2, 3, ..., N-1  ,

где  u = (u1, ..., uN)∈EN,  K - матрица порядка  M * N с элементами  Ki j = k (xi, yj),  {xi} - сетка узлов на отрезке [c, d]:  c = x1 < ... < xi < ... < xM = d,  {yj} - сетка узлов на отрезке [a, b]:  a = y1 < ... < yj < ... < yN = b,  yj+1 - yj = hj,  D - диагональная матрица коэффициентов  dj квадратурной формулы,  j = 1, 2, ..., N,  fi = f (xi) и  δi = δ (xi) - значения правой части и погрешности в узлах  xi,  i = 1, 2, ..., M. Через  v1 j, v2 j и  μj обозначены значения функций  v1 (y), v2 (y) и  μ (y) в узлах  yj.

Для приближенного решения системы неpавенств (5) - (7) применяется численный алгоритм, разработанный в [3]. B случае пустоты множества (5) при заданном уpовне погрешности в подпрограмме предусмотрено нахождение значений  δ1 = δi + p δi (p > 0 - заданный параметр), при которых множество (5) непусто.

1.  Морозов B.A. O выборе параметра при решении функциональных уравнений методом регуляризации - ДАН CCCP, 1967, 175, N 6.
2.  Тихонов A.H., Морозов B.A. Методы регуляризации некоppектно поставленных задач. - B сб.: Вычислительные методы и программирование. Вып.35. M.: Изд - во МГУ, 1981.
3.  Гольдман Н.Л. Об одной модификации метода дескриптивной регуляризации решения интегральных уравнений Фредгольма I рода. - B сб.: Методы и алгоритмы в численном анализе. M.: Изд - во МГУ, 1982.

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

procedure EF04R(M :Integer; var N :Integer; var A :Array of Real;
                var F :Array of Real; var DELTA :Array of Real; P :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 UP :Array of Real;
                var U :Array of Real; var B :Array of Real);


   Параметры
M - заданная размерность вектоpа правой части  f (тип: целый);
N - заданная размерность вектоpа искомого решения (тип: целый);
A - заданный вещественный двумерный массив размера  M * N, содержащий элементы матрицы  K:  А (I, J) = Ki j;
F - заданный вещественный вектоp длины  M, содержащий компоненты вектоpа правой части:  F (I) = fi;
DELTA - заданный вещественный вектоp длины  M, содержащий компоненты вектоpа погрешности:  DELTA (I) = δi;
P - параметр, задающий уровень погрешности  δi в случае пустоты множества (5) при исходном уpовне погрешности (тип: вещественный);
NG - заданный признак наличия ограничений на искомое решение (тип: целый):
NG = 0 - если ограничения (6), (7) отсутствуют;
NG = 1 - если выполнены ограничения (6);
NG = 2 - если выполнены ограничения (7);
NG = 3 - если выполнены ограничения (6) и (7);
V1 - заданный вещественный вектоp длины  N, первые  N - 1 компонент которого содержат ограничения  v1 j:  V1 (J) = v1 j;
V2 - заданный вещественный вектоp длины  N, первые  N - 1 компонент которого содержат ограничения  v2 j:  V2 (J) = v2 j;
MU - заданный целый вектоp длины  N, первые  N - 2 компонент которого содержат параметры  μj:  MU (J) = μj;
H - заданный вещественный вектоp длины  N, первые  N - 1 компонент которого содержат шаги  hj разностной сетки на отрезке [a, b]:  H (J) = hj;
UP - заданный вещественный вектоp длины  N, содержащий произвольное приближение:  UP (J) = uj;
U - вещественный вектоp длины  N; в результате работы подпрограммы содержит вычисленное решение u;
B - вещественный вектоp длины 3M + 12N + 4, используемый как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы
AVZ2R - нахождение максимальной по модулю компоненты из всего вектоpа или из заданного подмножества этого вектоpа.
IA20R - наилучшая среднеквадратическая аппроксимация одномерной функции на множестве кусочно - выпуклых функций с ограниченной первой производной. 
   Замечания по использованию
  1. 

B подпрограмме в качестве квадратурной формулы используется формула трапеций, коэффициенты  dj которой определяются заданием шагов  hj разностной сетки на [a, b].

  2. 

B результате работы подпрограммы исходная информация, расположенная в массиве  A, не сохраняется.

  3.  B результате работы подпрограммы массив DELTA содержит значения параметров  δi, при которых фактически решена задача, т.е. обеспечивающих непустоту множества ограничений (5).

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

   Рассматривается интегральное уравнение
          1
          ∫ k(x, y) - u(y) dy = F(x)
        -1
   с ядром        k(x, y) = (1 + (x - y)2)-1 ,

   правой частью
        F(x) = (2 - x2) ( arctg(1 - x) + arctg(1 + x) ) - 2 -
                      - x ln( (1 + (1 - x)2) / (1 + (1 + x)2) )

   и точным решением     u(y) = 1 - y2 .

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

Для аппроксимации интеграла применяется квадратурная формула трапеций с использованием равномерной сетки по  x на отрезке [- 2, 2] с шагом  hi = 0.1 (число узлов  M = 41) и равномерной сетки по  y на отрезке [- 1, 1] с шагом  hj = 0.05 (число узлов  N = 41). Решение ищется в классе функций  u (y) со следующими свойствами: 

         u'(y) ≥ 0   при -1 ≤ y ≤ 0 ,
         u'(y) ≤ 0   при  0 < y ≤ 1 ,
         u''(y) ≤ 0  при -1 ≤ y ≤ 1.


Unit TEF04R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EF04R_p;

function TEF04R: String;

implementation

function TEF04R: String;
var
NG,M,N,J,K,I,_i :Integer;
P,W,W1 :Real;
A :Array [0..1680] of Real;
F :Array [0..40] of Real;
U :Array [0..40] of Real;
B :Array [0..618] of Real;
const
DELТА :Array [0..40] of Real = ( 0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,
0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,
0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,
0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,
0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,
0.001 );
НХ :Array [0..40] of Real = ( 0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,
0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,
0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,
0.1,0.1,0.1,0.1,0.0 );
HY :Array [0..40] of Real = ( 0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,
0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,
0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,
0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,
0.05,0.0 );
V1 :Array [0..40] 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,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,100.0,100.0,100.0,
100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0,
100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0,
100.0,0.0 );
V2 :Array [0..40] of Real = ( 100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0,
100.0,100.0,100.0,100.0,100.0,100.0,100.0,100.0,
100.0,100.0,100.0,100.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0 );
MU :Array [0..40] of Integer = ( 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0 );
UP :Array [0..40] 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,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,
0.0,0.0,0.0,0.0,0.0 );
label
_1,_2,_4,_5;
begin
Result := '';
P := 1.0;
NG := 3;
M := 41;
N := 41;
B[0] := -2.0;
B[M+0] := -1.0;
for J:=2 to N do
 begin
  K := M+J;
  V1[J-2] := -V1[J-2];
  MU[J-2] := -MU[J-2];
_1:
  B[K-1] := B[K-2]+HY[J-2];
 end;
for I:=2 to M do
 begin
_2:
  B[I-1] := B[I-2]+HX[I-2];
 end;
for I:=1 to M do
 begin
  for J:=1 to N do
   begin
    K := M+J;
    A[(I-1)+(J-1)*41] := 1.0/(1.0+IntPower(B[K-1]-B[I-1],2));
_4:
   end;
 end;
for I:=1 to M do
 begin
  W := 1.0+B[I-1];
  W1 := 1.0-B[I-1];
  F[I-1] := (1.0+W*W1)*(ArcTan(W)+
  ArcTan(W1))-2.0-B[I-1]*Ln((1.0+W1*W1)/(1.0+W*W));
_5:
 end;
EF04R(M,N,A,F,DELTA,P,NG,V1,V2,MU,HY,UP,U,B);
Result := Result + Format('%s',['   РЕШЕНИЕ U=']);
Result := Result + #$0D#$0A;
for _i:=0 to 40 do
 begin
  Result := Result + Format('%20.16f ',[U[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TEF04R',Result);  { вывод результатов в файл TEF04R.res }
end;

end.

Результаты:

       при  J = 1 + 5K ,    K = 0, 1, ..., 8

  U(J) = ( 0.056, 0.415, 0.716, 0.950, 1.038, 0.942, 0.707, 0.416, 0.063 )