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

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

Назначение

Решение интегрального уравнения Фредгольма 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

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

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

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

1)  Невязка
                      d           b
     ρ(α)  =  (  ∫ P(y) [  ∫  K(y, x) uα(x) dx - f(y) ]2 dy )1/2 ;
                     c           a
2)  Функция  "чувствительности"
     τ(α)  =  α Ω1/2 (duα / dα) ;

3)  Hоpма решения
      γ(α)  =  Ω1/2 (uα ) ;

4)  Значение регуляризирующего функционала
      φ(α)  =  ρ2(α) + αγ2(α) ;

5)  Функция отношения
        ψ(α)  =  ρ2(α) / ρ(α) ,
где
                        d           b
     ρ2(α)  =  {  ∫ P(y) [  ∫  K(y, x) ( α d uα(x)/dα - uα(x) ) dx - f(y) ]2 dy }1/2 .
                       c           a

Здесь  uα = uα (x) - функция, минимизирующая (2).

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

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

1)  Выполнить только приведение задачи к каноническому виду, не решая интегрального уравнения (это требует порядка  M * N2 операций).

При работе с задачей, приведенной к каноническому виду, подпрограмма позволяет:

2)  Находить приближенное решение интегрального уравнения и значения критериальных функций при заданном значении паpаметpа регуляризации (это требует порядка  N2 операций).
3)  При заданном значении  α находить только значения критериальных функций, не находя решения (это требует порядка N операций).
4)  Если необходимо решить уравнение (1) с другой правой частью, но с тем же ядром, то подпрограмма позволяет выполнить необходимое преобразование новой задачи к "канонической" за число операций порядка  M * N.

1. Тихонов A.H., O регуляризации некоppектно поставленных задач, ДАН CCCP, 1963, т.153, N 1.

2. Воеводин B.B., O методе регуляризации, ЖВМ и МФ, 1969, т.9, N 3.

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

    int ef13r_c(real *sx, real *sy, integer *n, integer *m,
        real *a, real *u, real *f, real *p, real *p1, real *p2, real *p3,
        real *q, real *h, real *w, integer *l)

   Параметры
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   		параметры, определяющие стабилизирующий функционал (тип: вещественный);
q - заданный параметр регуляризации (тип: вещественный);
h - вещественный вектоp длины 50 ; в результате работы подпрограммы содержит вычисленные значения критериальных функций при заданном значении  q:  h (1) = ρ (α),  h (2) = τ (α),  h (3) = γ (α),  h (4) = φ (α),  h (5) = ρ1 (α);
w - двумерный вещественный рабочий массив размеpа  n * 8;
l - параметр, определяющий режим использования подпрограммы (тип: целый);
l = 1 - выполняется сведение задачи к канонической;
l = 2 - вычисляется решение и значения критериальных функций;
l = 3 - вычисляются только значения критериальных функций;
l = 4 - выполняются дополнительные преобразования, приводящие задачу с новой правой частью к канонической: новая правая часть должна быть расположена на месте f. 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1) 

Подпрограмма ef13r_c обращается к вспомогательным подпрограммам с именами: effr_c, efrs_c, efus_c, efuf_c, efsl_c, efud_c, efrd_c.

  2)  При работе подпрограммы в режиме, выполняющем дополнительные преобразования, приводящие задачу с новой правой частью к "канонической" (p;4). Значения остальных перечисленных выше параметров (и рабочего массива  w) должны сохраняться, т.к. они используются в программе при всех значениях  l.

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

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

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

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

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие знначения критериальных функций при фиксированном параметре регуляризации, с весовой функцией  p ≡ 1 и использованием равномерных по  x и  y на отрезке [0,1] сеток из 11 точек. 

int main(void)
{
    /* Initialized data */
    static float p[11] = { 1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f };
    static int n = 11;
    static int m = 11;
    static float q = 1e-6f;
    static float x[11] = { 0.f,.1f,.2f,.3f,.4f,.5f,.6f,.7f,.8f,.9f,1.f };
    static float y[11] = { 0.f,.1f,.2f,.3f,.4f,.5f,.6f,.7f,.8f,.9f,1.f };

    /* Builtin functions */
    double atan(double), log(double);

    /* Local variables */
    extern  int ef13r_c(float *, float *, int *, int *, float *, float *,
            float *, float *, float *, float *, float *, float *,
            float *, float *, int *);
    static float a[121]  /* was [11][11] */, f[11], h__[50];
    static int i__, i, j;
    static float u[11], w[88]  /* was [11][8] */;
    int i__1, i__2;

#define a_ref(a_1,a_2) a[(a_2)*11 + a_1 - 12]

    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        f[i__ - 1] = (float)atan(x[i__ - 1]);
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
/* l1: */
            a_ref(i__, j) = y[i__ - 1] / (x[j - 1] * x[j - 1] * y[i__ - 1] * 
                    y[i__ - 1] + 1.f);
        }
    }
    ef13r_c(x, y, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &q, h__, w, &c__1);
    ef13r_c(x, y, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &q, h__, w, &c__2);

    for (i = 0; i <= 6; i+= 3) {
       printf("\n %16.7e %16.7e %16.7e ",u[i], u[i+1], u[i+2]);
    }
    printf("\n %16.7e %16.7e \n",u[9], u[10]);
    for (i = 1; i <= 5; ++i) {
        printf("\n %16.7e \n",h__[i-1]);
    }
    f[0] = 0.f;
    i__2 = m;
    for (i__ = 2; i__ <= i__2; ++i__) {
/* l2: */
        f[i__ - 1] = (float)log((float)(y[i__ - 1] * y[i__ - 1] + 1.f)) /
                (y[i__ - 1] * 2.f);
    }
    q = 1e-6f;
    ef13r_c(x, y, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &q, h__, w, &c__4);
    ef13r_c(x, y, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &q, h__, w, &c__2);

    for (i = 0; i <= 6; i+= 3) {
       printf("\n %16.7e %16.7e %16.7e ",u[i], u[i+1], u[i+2]);
    }
    printf("\n %16.7e %16.7e \n",u[9], u[10]);
    for (i = 1; i <= 5; ++i) {
        printf("\n %16.7e \n",h__[i-1]);
    }
    return 0;
} /* main */


Результаты:

       u  =  ( 0.99, 0.99, 0.99, 1.0, 1, 0, 1.0, 1.0, 1.0, 0.99, 0.98, 0.96 )
       h  =  ( 5.6*10-6,  5.7*10-4,  1.0,  1.0*10-6,  1.7*105 )



       for (i__ = 2; i__ <= i__2; ++i__) {
           f[i__ - 1] = (float)log((float)(y[i__ - 1] * y[i__ - 1] + 1.f)) /
                   (y[i__ - 1] * 2.f);
       }
       q = 1e-6f;
       ef13r_c(x, y, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &q, h__, w, &c__4);
       ef13r_c(x, y, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &q, h__, w, &c__2);

Результаты:

       u    = ( 0.04, 0.10, 0.19, 0.296, 0.39, 0.48, 0.59, 0.72, 0.85, 0.92, 0.84 )
       h__  = ( 9.8*10-8, 1.7*10-2, 0.577, 3.*10-12, 4.6*10-6 )