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

Назначение

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

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

    int ef03r_c(integer *m, integer *n, real *a, real *f,
        real *delta, real *h, real *d, real *p1, real *p2, real *p3,
        real *p4, integer *l, real *u, real *b)

   Параметры 

   Версии:  нет 

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

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

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

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

               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)

int main(void)
{
    /* Initialized data */
    static float delta[41] = { .001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,
            .001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,
            .001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,
            .001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f,.001f 
            };
    static float p4 = 1.f;
    static int l = 1;
    static float hj[41] = { .05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,
            .05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,
            .05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,
            .05f,.05f,.05f,.05f,0.f };
    static float hi[41] = { .1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,
            .1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,
            .1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,0.f };
    static int m = 41;
    static int n = 41;
    static float d__[41] = { .025f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,
            .05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,
            .05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,.05f,
            .05f,.05f,.05f,.05f,.05f,.025f };
    static float p1 = 1.f;
    static float p2 = 0.f;
    static float p3 = 0.f;

    /* System generated locals */
    int i__1, i__2;
    float r__1;

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

    /* Local variables */
    extern int ef03r_c(int *, int *, float *, float *, float *, float *,
             float *, float *, float *, float *, float *, int *,
             float *, float *);
    static float a[1681] /* was [41][41] */, b[328], f[41];
    static int i__, j, k;
    static float u[41], w, w1;

#define a_ref(a_1,a_2) a[(a_2)*41 + a_1 - 42]

    b[0] = -2.f;
    b[m] = -1.f;
    i__1 = n;
    for (j = 2; j <= i__1; ++j) {
        k = m + j;
/* l1: */
        b[k - 1] = b[k - 2] + hj[j - 2];
    }
    i__1 = m;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* l2: */
        b[i__ - 1] = b[i__ - 2] + hi[i__ - 2];
    }
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            k = m + j;
/* Computing 2nd power */
            r__1 = b[k - 1] - b[i__ - 1];
            a_ref(i__, j) = 1.f / (r__1 * r__1 + 1.f);
/* l3: */
        }
    }
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__) {
        w = b[i__ - 1] + 1.f;
        w1 = 1.f - b[i__ - 1];
        f[i__ - 1] = (w * w1 + 1.f) * ((float)atan(w) + (float)atan(w1)) -
                2.f - b[i__ - 1] * (float)log((float)((w1 * w1 + 1.f) /
                (w * w + 1.f)));
/* l4: */
    }
    ef03r_c(&m, &n, a, f, delta, hj, d__, &p1, &p2, &p3, &p4, &l, u, b);

    for (j = 1; j <= 41; j += 5) {
        i__ = 42 - j;
        printf("\n  i = %2i   u(i) = %10.3e \n",i__, u[i__ - 1]);
/* l6: */
    }
    return 0;
} /* main */


Результаты:

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