Подпрограмма:  EF03R (модуль EF03R_p)

Назначение

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

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

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

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

с известными ядром  K (x, y) и правой частью  f (x) сводится методом поточечной невязки [1,2] к минимизации функционала

                      b
(2)    Ω(u)  =  ∫ ( p1 u2 + p2 (du/dy)2 + p3 (d2u/dy2) ) dy
                     a
при условии

(3)     | Au - f(x) | ≤ δ(x) , 

где  δ (x) - заданная погрешность правой части уравнения (1), 
p₁ > 0,  p₂ ≥ 0,  p₃ ≥ 0 - заданные параметры.

Дискретный аналог задачи (2), (3) состоит в определении вектоpа  u = (u₁, ..., u_N) ∈ E^N, обладающего наименьшей  C - нормой

(4)      || u ||c2  =  (C u, u)  =  min u (C u, u) 

среди всех  u = (u₁, ..., u_N) ∈ E^N, для которых выполнены неpавенства

(5)      | (K D u)i - fi | ≤ δi ,    i ∈ 1, M 

Здесь  C - ленточная положительно определенная матрица квадратичной формы, аппроксимирующей функционал  Ω (u) на сетке  { y_j } :

a = y₁ < y₂ < ... < y_j < ... < y_N  =  b ,  y_j + 1 - y_j = h_j  с помощью разностных отношений:

 du/dy [ в точке yj+1/2 ]  =  [ u(yj+1) - u(yj) ] / hj ,

 d2u/dy2 [ в точке yj ]  = [ du/dy [ в точке yj+1/2 ]  -
                                    -   du/dy [ в точке yj-1/2 ]  ] / 0.5(hj + hj-1) 

и квадратурной формулы трапеций.  K - известная матрица порядка  M * N с элементами  k_i j = k (x_i, y_j),  {x_i} - сетка узлов на [c, d]:

c = x₁ < ... < x_i < ... < x_M = d ,  D - диагональная матрица коэффициентов квадратурной формулы  d_j ( j ∈ 1, N ),  f_i = f (x_i) и  δ_i = δ (x_i) - значения правой части погрешности в узлах  x_i.

Для приближенного решения задачи (4), (5) в подпрограмме применяется метод последовательного проектирования [3], позволяющий определить элемент  w = (w₁, ..., w_n) ∈ E^N такой, что

     || u ||c ≤ || w ||c ≤ 2 || u ||c . 

В случае пустоты множества (5) при заданных  δ_i в подпрограмме находятся значения  δ_i = δ_i + p₄ δ_i (p₄ > 0 - заданный параметр), при котоpом множество (5) непусто.

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

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

procedure EF03R(var M :Integer; var N :Integer; var A :Array of Real;
                var F :Array of Real; var DELTA :Array of Real;
                var H :Array of Real; var D :Array of Real; var P1 :Real;
                var P2 :Real; var P3 :Real; P4 :Real; L :Integer;
                var U :Array of Real; var B :Array of Real);

   Параметры 

   Версии:  нет 

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

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

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

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

               1
               ∫  K(x, y) u(y) dy  =  F(x)
             -1 

с ядром  K (x, y) = ( 1 + (x - y)² )^- 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 - y².

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

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

Unit TEF03R_p;
interface
uses
SysUtils, Math, { Delphi }
LStruct, Lfunc, UtRes_p, EF03R_p;

function TEF03R: String;

implementation

function TEF03R: String;
var
J,K,I :Integer;
W,W1 :Real;
A :Array [0..1680] of Real;
F :Array [0..40] of Real;
U :Array [0..40] of Real;
B :Array [0..327] 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 );
HJ :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 );
HI :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 );
M :Integer = 41;
N :Integer = 41;
D :Array [0..40] of Real = ( 0.025,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.025 );
P1 :Real = 1.0;
P2 :Real = 0.0;
P3 :Real = 0.0;
P4 :Real = 1.0;
L :Integer = 1;
label
_1,_2,_3,_4,_6;
begin
Result := '';
B[0] := -2.0;
B[M+0] := -1.0;
for J:=2 to N do
 begin
  K := M+J;
_1:
  B[K-1] := B[K-2]+HJ[J-2];
 end;
for I:=2 to M do
 begin
_2:
  B[I-1] := B[I-2]+HI[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));
_3:
   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));
_4:
 end;
EF03R(M,N,A,F,DELTA,HJ,D,P1,P2,P3,P4,L,U,B);
J := 1;
WHILE ( J<=41 ) do
 begin
  I := 42-J;
  Result := Result + Format('%s',['     I=']);
  Result := Result + Format('%3d',[I]) + #$0D#$0A;
  Result := Result + Format('%s',['   U[I-1]=']);
  Result := Result + Format('%20.16f',[U[I-1]]) + #$0D#$0A;
_6:
  inc(J,5);
 end;
UtRes('TEF03R',Result);  { вывод результатов в файл TEF03R.res }
exit;
end;

end.

Результаты:

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

       U(J)  =  ( 0.104,  0.411,  0.715,  0.952,  1.041,  0.944,  0.707, 
                      0.411,  0.119 )