Подпрограмма:  ef12r_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 

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

          ρ(α)  =  e || f ||p  ,
   где
                        d           b
       ρ(α)  =  (  ∫ p(y) [  ∫ k(y, x) ua(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  +  a u T C u .
         i=1               j=1 

Здесь  p_i = p (y_i),  u_j = u (x_j),  f_i= f (y_i),  k_i j = k (y_i, x_j);  a =x₁ < x₂ < ... < x_N = b;  c = y₁ < y₂ < ... < y_M = d - сетки узлов по  x и  y,  σ_j' и  σ_i'' - коэффициенты квадратурных формул по  x и  y соответственно; C - матрица квадратичной формы, аппроксимирующей стабилизирующий функционал  Ω (u);  u = (u₁, u₂, ..., u_N) - вектор значений искомого решения.

Пусть f = ( f₁ √σ₁'',  f₂ √σ₂'', ..., f_M √σ_M'' ) и  A = (a_i j),  a_i j = k_i j p_i σ_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) Решение интегрального уравнения для фиксированный правой части и относительного уровня невязки (это требует порядка  MN² операций).

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

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

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

    int ef12r_c (real *sx, real *sy, integer *n, integer *m,
            real *a, real *u, real *f, real *p, real *p1, real *p2, real *p3,
            real *eps, real *hx, real *w, integer *l)

   Параметры 

   Версии:  нет 

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

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

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

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

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

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

Ниже приводится фрагмент программы, вычисляющей решения и соответствующие значения критериальных функций при заданных значениях относительного уровня невязки  eps с весовой функцией  p (y) = 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 eps = .001f;
    static float sx[11] = { 0.f,.1f,.2f,.3f,.4f,.5f,.6f,.7f,.8f,.9f,1.f };
    static float sy[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 ef12r_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__[7];
    static int 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(sy[i__ - 1]);
	i__2 = n;
	for (j = 1; j <= i__2; ++j) {
/* l1: */
	    a_ref(i__, j) = sy[i__ - 1] / (sx[j - 1] * sx[j - 1] * sy[i__ - 1]
		     * sy[i__ - 1] + 1.f);
	}
    }
    ef12r_c(sx, sy, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &eps, h__, w,
           &c__1);
           
    for (j = 0; j <= 6; j += 3) {       
         printf("\n %16.7e %16.7e %16.7e ",u[j], u[j+1], u[j+2]);
    }     
    printf("\n %16.7e %16.7e \n\n",u[9], u[10]);
    for (j = 1; j <= 7; ++j) {
         printf("\n %16.7e \n",h__[j]);
    }     
    f[0] = 0.f;
    i__2 = n;
    for (i__ = 2; i__ <= i__2; ++i__) {
/* l2: */
	f[i__ - 1] = (float)log((float)(sy[i__ - 1] * sy[i__ - 1] + 1.f)) /
	        (sy[i__ - 1] * 2.f);
    }
    ef12r_c(sx, sy, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &eps, h__, w,
           &c__3);
 
    for (j = 0; j <= 6; j += 3) {
         printf("\n %16.7e %16.7e %16.7e ",u[j], u[j+1], u[j+2]);
    }
    printf("\n %16.7e %16.7e \n\n",u[9], u[10]);
    for (j = 1; j <= 7; ++j) {
         printf("\n %16.7e \n",h__[j]);
    }
    return 0;
} /* main */


Результаты:

       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 )



       f[0] = 0.f;
       i__2 = n;
       for (i__ = 2; i__ <= i__2; ++i__) {
          f[i__ - 1] = (float)log((float)(sy[i__ - 1] * sy[i__ - 1] + 1.f)) /
	        (sy[i__ - 1] * 2.f);
       }
       ef12r_c(sx, sy, &n, &m, a, u, f, p, &c_b2, &c_b3, &c_b3, &eps, h__, w,
           &c__3);

Результаты:

       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 )