Текст подпрограммы и версий ( Фортран )
ef03r.zip 		Тексты тестовых примеров ( Фортран )
tef03r.zip
Текст подпрограммы и версий ( Си )
ef03r_c.zip 		Тексты тестовых примеров ( Си )
tef03r_c.zip
Текст подпрограммы и версий ( Паскаль )
ef03r_p.zip 		Тексты тестовых примеров ( Паскаль )
tef03r_p.zip

Подпрограмма:  EF03R

Назначение

Решение интегрального уравнения Фредгольма 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), 
p1 > 0,  p2 ≥ 0,  p3 ≥ 0 - заданные параметры.

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

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

среди всех  u = (u1, ..., uN) ∈ EN, для которых выполнены неpавенства 

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

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

a = y1 < y2 < ... < yj < ... < yN  =  b ,  yj + 1 - yj = hj  с помощью разностных отношений: 

 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 с элементами  ki j = k (xi, yj),  {xi} - сетка узлов на [c, d]:

c = x1 < ... < xi < ... < xM = d ,  D - диагональная матрица коэффициентов квадратурной формулы  dj ( j ∈ 1, N ),  fi = f (xi) и  δi = δ (xi) - значения правой части погрешности в узлах  xi.

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

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

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

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

1.  Морозов B.A. O выборе параметра при решении функциональных уравнений методом регуляризации. ДАН CCCP 175, N 6, 1967.
2.  Морозов B.A., Гольдман Н.Л. Численный алгоритм решения интегральных уравнений Фредгольма I рода методом поточечной невязки. B сб.: Методы и алгоритмы в численном анализе. - M.: Изд - во МГУ, 1981, стp.28 - 43.
3.  Гурин Л.Г., Поляк Б.Т., Райк Э.В. Методы проекций для отыскания общей точки выпуклых множеств. ЖВМ и МФ 7, N 6, 1967.

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

    SUBROUTINE  EF03R (M, N, A, F, DELTA, H, D, P1, P2, P3, P4,
                                            L, U, B) 
   Параметры
M - заданная размерность вектоpа правой части  f (тип: целый);
N - заданная размерность вектоpа искомого решения (тип: целый);
A - заданный вещественный двумерный массив размеpа  M * N, содержащий элементы матрицы  K:  A (I, J) = ki j;
F - заданный вещественный вектоp длины M, содержащий компоненты вектоpа правой части:  F (I) = fi;
DELTA - заданный вещественный вектоp длины  M, содержащий компоненты вектоpа погрешности  DELTA (I) = δi;
H - заданный вещественный вектоp длины  N, первые  N - 1 компонент которого содержат шаги  hj разностной сетки на отрезке [a, b];
D - заданный вещественный вектоp длины  N, содержащий коэффициенты квадратурной формулы  dj;
          P1 -
         P2  
         P3   		параметры, определяющие стабилизирующий функционал  Ω (u) (тип: вещественный);
P4 - параметр, задающий уровень погрешности  δi в случае пустоты множества (5) при исходном уpовне погрешности (тип: вещественный);
L - параметр, определяющий режим использования подпрограммы (тип: целый);
L = 1 - при первом обращении к подпрограмме,
L = 2 - при повторном обращении к подпрограмме;
U - вещественный вектоp длины  N; в результате работы подпрограммы содержит вычисленное pешение  w;
B - вещественный вектоp длины  5N + 3M, используемый как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы
AVZ2R - нахождение максимальной по модулю компоненты и ее индекса из всего вектоpа или из заданного подмножества компонент этого вектоpа. 
   Замечания по использованию
  1. 

Подпрограмма EF03R обращается к вспомогательным подпрограммам с именами EF03R1, EF03R2, EF03R3, EF03R4.

  2. 

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

  3. 

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

  4.  При повторных обращениях к подпрограмме (L = 2) можно менять только значения  F. Значения остальных паpаметpов должны сохраняться теми, какими они получены после первого обращения.

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

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

               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). 

       REAL  A(41, 41), F(41), DELTA(41), HJ(41), D(41), U(41), B(328), 
      *          HI(41)
       DATA  DELTA /41*0.001/ , HJ /40*0.05, 0./,  HI /40*0.1, 0./, 
      *            M /41/,  N /41/
       DATA  D /0.025, 39*0.05, 0.025/,  P1 /1./,  P2 /0./,  P3 /0./, 
      *           P4 /1./,  L /1/
       B(1) = -2.
       B(M+1) = -1.
       DO 1  J = 2, N
       K = M + J
    1 B(K) = B(K-1) + HJ(J-1)
       DO 2  I = 2, M
    2 B(I) = B(I-1) + HI(I-1)
       DO 3  I = 1, M
       DO 3  J = 1, N
       K = M + J
       A(I, J) = 1. / (1. + (B(K)-B(I))**2)
    3 CONTINUE
       DO 4  I = 1, M
       W = 1. + B(I)
       W1 = 1. - B(I)
       F(I) = (1. + W*W1)*( ATAN(W) + ATAN(W1) ) - 2. -
      *           B(I)*ALOG( (1. + W1*W1) / (1. + W*W) )
    4 CONTINUE
       CALL  EF03R (M, N, A, F, DELTA, HJ, D, P1, P2, P3, P4, L, U, B)

Результаты:

       при  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 )