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

Назначение

Решение одномерного интегрального уравнения Винера - Хопфа первого рода с симметричным неотрицательно определенным ядром методом регуляризации первого порядка без выбора параметра регуляризации.

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

Приближенное решение интегpaльнoго уравнения Винера - Хопфа первого рода

                    ∞   
(1)    Au  ≡   ∫    K(x - z) u(t) dt  =  f(x)  ,
                    0
          x ∈ [0, +∞) ,    K(x - t)  =  K(t - x)  , 

с гладким неотрицательно определенным ядром K осуществляется методом упрощенной однопараметрической регуляризации А.Н.Тихонова [1] путем сведения задачи к минимизации по  u  сглаживающего функционала

             ∞   ∞
(2)        ∫     ∫   K(x - t) u(x) u(t) dx dt  -
           0     0
                           ∞                             ∞
                    -  2  ∫   f(x) u(x) dx  +  α ∫   u ' 2(x) dx  ,
                          0                              0
   где  α > 0  -  параметр регуляризации ,
      ∞
      ∫   u ' 2(x) dx  -  стабилизирующий функционал. 
     0 

Считается, что  K (t)→0,  u (t)→0,  f (t)→0 (достаточно быстро) при  t→∞, так что функции  K (t),  f (t) можно считать заданными на конечном отрезке [0, T].

Тогда вместо (2) будем иметь:

             T    T
(3)        ∫     ∫   K(x - t) u(x) u(t) dx dt  -
           0     0
                            T                             T
                    -  2  ∫   f(x) u(x) dx  +  α ∫   u ' 2(x) dx  ,
                          0                              0 

Для дискретизации первого и второго слагаемого (3) используется квадратурная формула трапеций на равномерной сетке  0 = x₁ < x₂ < ... < x_N = T ,  x_i + 1 - x_i = δx , а при аппроксимации третьего слагаемого - разностные отношения:

     du/dx (xi)  =  ( u(xi+1) - u(xi) ) / δx ,        i = 1, 2, ..., N-1 . 

B предположении  du/dx (x_N) = 0 дискретный аналог (3) имеет вид:

(4)    (D K D u, u - 2(D f, u) + α (C u, u).

   Здесь: 

D = diag {d₁₁, d₂₂, ..., d_NN},  причем  d₁₁ = d_NN = δx /2 ,  d_i i = δx ,  i = 2, 3, ..., N-1 ;

K - симметричная теплицева матрица, определенная элементами первой строки  K_1 i = K (x_i) ;

u = (u₁, u₂, ..., u_N) ,  u_i = u (x_i) ,

f = (f₁, f₂, ..., f_N) ,  f_i = f (x_i) ,  i = 1, 2, ..., N ,

C - трехдиагональная матрица, аппроксимирующая стабилизирующий оператор.

Задача минимизации (4) по  u эквивалентна решению системы линейных алгебраических уравнений:

(5)     D K D u + α C u  =  D f . 

При ее решении используется схема, предложенная в [2], в основе которой лежит существенное использование теплицевости матрицы  K и "квазитеплицевости" матрицы  C.

Алгоритм, описанный в [3], был адаптирован для решения систем уравнений с симметричной теплицевой матрицей. Эта адаптация дает почти двойной выигрыш по времени.

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

    int ec21r_c (real *a, real *x, real *ba, real *b, real *c,
            integer *n, real *t, integer *j1, real *alfa)

   Параметры 

   Версии:  нет 

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

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

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

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

    k(t) = exp( -| t | )
    и правой частью  f(t) = 2/5 exp(-t) ( cos(t) + 2 sin(t) )
    ( точное решение  u(t) = exp(-t) cos(t) ). 

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

int main(void)
{
    /* Builtin functions */
    double exp(double), cos(double), sin(double);

    /* Local variables */
    static float alfa;
    extern int ec21r_c(float *, float *, float *, float *, float *, int *,
                       float *, int *, float *);
    static float uton[200], a[200], b[200], c__[200], h__;
    static int i__, n;
    static float t, x[200];
    static int j1;
    static float ba[200];
    int i__1;

    t = 10.f;
    n = 200;
    h__ = t / (float) (n - 1);
    j1 = 0;
    alfa = 1.0000000000000001e-9f;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        uton[i__ - 1] = (float)exp((float)((1.f - i__) * h__)) *
                (float)cos((float)((i__ - 1.f) * h__));
/* l4: */
    }
    printf("\n %13.4e %13.4e %13.4e %13.4e ",
           uton[0], uton[9], uton[39], uton[69]);
    printf("\n %13.4e %13.4e %13.4e %13.4e \n",
           uton[99], uton[129], uton[159], uton[189]);

    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        a[i__ - 1] = (float)exp((float)((1.f - i__) * h__));
        ba[i__ - 1] = (float)exp((float)((1.f - i__) * h__)) *
                .40000000000000002f * ((float)cos((float)((i__ - 1.f) * h__)) +
                (float)sin((float)((i__ - 1.f) * h__)) * 2.f);
/* l1: */
    }
    ec21r_c(a, x, ba, b, c__, &n, &t, &j1, &alfa);

    printf("\n %13.4e %13.4e %13.4e %13.4e ",
           x[0], x[9], x[39], x[69]);
    printf("\n %13.4e %13.4e %13.4e %13.4e \n",
           x[99], x[129], x[159], x[189]);
    return 0;
} /* main */


Результаты:

       x(1)      =   0.9830 ,
       x(10)    =   0.5722 ,
       x(40)    =  -0.0534 ,
       x(70)    =  -0.0296 ,
       x(100)  =   0.0018 ,
       x(130)  =   0.0015 ,
       x(160)  =  -0.0000 ,
       x(190)  =  -0.0001 .