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

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

Назначение

Решение интегрального уравнения Фредгольма I рода методом поточной невязки с дескриптивной регуляризацией.

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

Приближенное решение интегpaльнoго уравнения Фредгольма I рода 

                     b   
(1)    Au  ≡   ∫    k(x, y) u(y) dy = f(x)  ,         c ≤ x ≤ d ,
                    a

с гладким ядром и правой частью  f (x), заданной с погрешностью  δ (x), при априорно известной информации об ограниченности первой призводной исходного решения  u (y) и об участках знакопостоянства его кривизны сводится методом поточной невязки с дескриптивной регуляризацией [1,2] к отысканию общей точки  u (y) множеств  Uδ , Uν , Uμ : 

(2)     Uδ = { u:  | Au - f(x) | ≤ δ(x) }  ,
(3)     Uν = { u:  v1(y) ≤ u'(y) ≤ v2(y)   ∀y ∈ [a, b] }  ,
(4)     Uμ = { u:  μ(y) u''(y) ≥ 0   ∀y ∈ [a, b] ,  μ(y) = sign u''(y) }

т.е. к выделению из множества  Uδ формальных решений уравнения (1) приближенного решения  u (y), удовлетворяющего условиям (3), (4). v1 (y) и  v2 (y) - заданные на отрезке [a, b] функции.

Дискретный аналог задачи (2) - (4) состоит в определении вектоpа  u = (u1, ..., uN)∈EN, удовлетворяющего системе линейных неpавенств: 

(5)     | (K Du)i - fi | ≤ δi ,      i = 1, 2, ..., M  ,
(6)     v1 j ≤ (uj+1 - uj) / hj ≤ v2 j ,      j = 1, 2, ..., N-1  ,
(7)     μj ( (uj+1 - uj) / hj - (uj - uj-1) / hj-1 ) ≥ 0 ,      j = 2, 3, ..., N-1  ,

где  u = (u1, ..., uN)∈EN,  K - матрица порядка  M * N с элементами  Ki j = k (xi, yj),  {xi} - сетка узлов на отрезке [c, d]:  c = x1 < ... < xi < ... < xM = d,  {yj} - сетка узлов на отрезке [a, b]:  a = y1 < ... < yj < ... < yN = b,  yj+1 - yj = hj,  D - диагональная матрица коэффициентов  dj квадратурной формулы,  j = 1, 2, ..., N,  fi = f (xi) и  δi = δ (xi) - значения правой части и погрешности в узлах  xi,  i = 1, 2, ..., M. Через  v1 j, v2 j и  μj обозначены значения функций  v1 (y), v2 (y) и  μ (y) в узлах  yj.

Для приближенного решения системы неpавенств (5) - (7) применяется численный алгоритм, разработанный в [3]. B случае пустоты множества (5) при заданном уpовне погрешности в подпрограмме предусмотрено нахождение значений  δ1 = δi + p δi (p > 0 - заданный параметр), при которых множество (5) непусто.

1.  Морозов B.A. O выборе параметра при решении функциональных уравнений методом регуляризации - ДАН CCCP, 1967, 175, N 6.
2.  Тихонов A.H., Морозов B.A. Методы регуляризации некоppектно поставленных задач. - B сб.: Вычислительные методы и программирование. Вып.35. M.: Изд - во МГУ, 1981.
3.  Гольдман Н.Л. Об одной модификации метода дескриптивной регуляризации решения интегральных уравнений Фредгольма I рода. - B сб.: Методы и алгоритмы в численном анализе. M.: Изд - во МГУ, 1982.

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

    int ef04r_c(integer *m, integer *n, real *a, real *f,
        real *delta, real *p, integer *ng, real *v1, real *v2, integer *mu,
        real *h, real *up, real *u, real *b)

   Параметры
m - заданная размерность вектоpа правой части  f (тип: целый);
n - заданная размерность вектоpа искомого решения (тип: целый);
a - заданный вещественный двумерный массив размера  m * n, содержащий элементы матрицы  K:  a (i, j) = ki j;
f - заданный вещественный вектоp длины  m, содержащий компоненты вектоpа правой части:  f (i) = fi;
delta - заданный вещественный вектоp длины  m, содержащий компоненты вектоpа погрешности:  delta (i) = δi;
p - параметр, задающий уровень погрешности  δi в случае пустоты множества (5) при исходном уpовне погрешности (тип: вещественный);
ng - заданный признак наличия ограничений на искомое решение (тип: целый):
ng = 0 - если ограничения (6), (7) отсутствуют;
ng = 1 - если выполнены ограничения (6);
ng = 2 - если выполнены ограничения (7);
ng = 3 - если выполнены ограничения (6) и (7);
v1 - заданный вещественный вектоp длины  n, первые  n - 1 компонент которого содержат ограничения  v1 j:  v1 (j) = v1 j;
v2 - заданный вещественный вектоp длины  n, первые  n - 1 компонент которого содержат ограничения  v2 j:  v2 (j) = v2 j;
mu - заданный целый вектоp длины  n, первые  n - 2 компонент которого содержат параметры  μj:  mu (j) = μj;
h - заданный вещественный вектоp длины  n, первые  n - 1 компонент которого содержат шаги  hj разностной сетки на отрезке [a, b]:  h (j) = hj;
up - заданный вещественный вектоp длины  n, содержащий произвольное приближение:  up (j) = uj;
u - вещественный вектоp длины  n; в результате работы подпрограммы содержит вычисленное решение u;
b - вещественный вектоp длины 3m + 12n + 4, используемый как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы
avz2r_c - нахождение максимальной по модулю компоненты из всего вектоpа или из заданного подмножества этого вектоpа.
ia20r_c - наилучшая среднеквадратическая аппроксимация одномерной функции на множестве кусочно - выпуклых функций с ограниченной первой производной. 
   Замечания по использованию
  1. 

B подпрограмме в качестве квадратурной формулы используется формула трапеций, коэффициенты  dj которой определяются заданием шагов  hj разностной сетки на [a, b].

  2. 

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

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

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

   Рассматривается интегральное уравнение
          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) Решение ищется в классе функций  u (y) со следующими свойствами: 

         u'(y) ≥ 0   при -1 ≤ y ≤ 0 ,
         u'(y) ≤ 0   при  0 < y ≤ 1 ,
         u''(y) ≤ 0  при -1 ≤ y ≤ 1.


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 hx[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 float hy[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 v1[41] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
            0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,100.f,100.f,100.f,100.f,100.f,
            100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,
            100.f,100.f,100.f,100.f,0.f };
    static float v2[41] = { 100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,
            100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,100.f,
            100.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
            0.f,0.f,0.f,0.f,0.f,0.f };
    static int mu[41] = { 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
            1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0 };
    static float up[41] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
            0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
            0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };

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

    /* Local variables */
    extern int ef04r_c(int *, int *, float *, float *,
            float *, float *, int *, float *, float *, int *, float *,
            float *, float *, float *);
    static float a[1681] /* was [41][41] */, b[619], f[41];
    static int i__, i, j, k, m, n;
    static float p, u[41], w, w1;
    static int ng;
    int i__1, i__2;
    float r__1;

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

    p = 1.f;
    ng = 3;
    m = 41;
    n = 41;
    b[0] = -2.f;
    b[m] = -1.f;
    i__1 = n;
    for (j = 2; j <= i__1; ++j) {
        k = m + j;
        v1[j - 2] = -v1[j - 2];
        mu[j - 2] = -mu[j - 2];
/* l1: */
        b[k - 1] = b[k - 2] + hy[j - 2];
    }
    i__1 = m;
    for (i__ = 2; i__ <= i__1; ++i__) {
/* l2: */
        b[i__ - 1] = b[i__ - 2] + hx[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);
/* l4: */
        }
    }
    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)));
/* l5: */
    }
    ef04r_c(&m, &n, a, f, delta, &p, &ng, v1, v2, mu, hy, up, u, b);

    for (i = 0; i <= 36; i+= 4) {
     printf("\n %16.7e %16.7e %16.7e %16.7e ",u[i], u[i+1], u[i+2], u[i+3]);
    }
    printf("\n %16.7e \n",u[40]);
    return 0;
} /* main */


Результаты:

       при  j = 1 + 5k ,    k = 0, 1, ..., 8

  u(j) = ( 0.056, 0.415, 0.716, 0.950, 1.038, 0.942, 0.707, 0.416, 0.063 )