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

Назначение

Приближенное решение интегрального уравнения Фредгольма I рода в классе кусочно - выпуклых функций с ограниченной первой производной.

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

Решение интегрального уравнeния фредгольма I рода

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

при априорно известной информации об ограниченности первой производной искомого решения  u = u (y) и об участках знакопостоянства его кривизны сводится на основе метода квазирешений В.К.Иванова к задаче минимизации

                   d      b   
     min        ∫   (  ∫  k(x, y)  u(y) dy - f(x) )2 dx
    v, u (a)    c      a 

на множестве кусочно - монотонных функций  v (y) = u' (y), определенных на отрезке [a, b]

                                 y
          u(y)  =  u(a) + ∫ v(t) dt ,   a ≤ y ≤ b ,
                               a
 и удовлетворяющих ограничениям
        v1(y) ≤ v(y) ≤ v2(y) ,   a ≤ y ≤ b , 

v₁ (y) и  v₂ (y) - заданные на отрезке [a, b] функции. Значение  u (a) может считаться как известным, так и подлежащим определению.

Численная дискретизация этой вариационной задачи после аппроксимации соответствующих интегралов с помощью квадратурных формул на сетках  Ω_x = {x_i :  c = x₁ < ... < x_i < ... < x_M = d} и Ω_y = {y_j :  a = y₁ < ... < y_j < ... < y_N = b} приводит к задаче квадратичного программирования

                        M            N
  (1)      min      ∑  pi  (   ∑  dj ki j uj  -  fi  )2
            v, u1    i=1         j=1 

при условии, что вектоp  v = (v₁, ..., v_N - 1) удовлетворяет ограничениям

                               j -1
  (2)     uj  =  u1  +   ∑   he ve ,     j = 2, ..., N
                               e=1
        
  (3)     v1 j ≤ vj ≤ v2 j ,     j = 1, ..., N-1

  (4)     mj ( vj + 1 - vj ) ≥ 0 ,     j = 1, ..., N-2 , 

где  k_i j = k (x_i, y_j),  f_i = f (x_i),  d_j > 0 - коэффициенты квадратурной формулы трапеции на сетке  Ω_y ,  p_i > 0 - весовые коэффициенты,  h_j = y_j + 1 - y_j - шаг сетки  Ω_y ,  v_1j = v₁ (y_j),  v_2j = v₂ (y_j), параметры  m_j определяются условиями кусочной монотонности функции  v = u' и принимают значения 1, 0 или -1.

Численное решение задачи (1) - (4) основано на методе проекции сопряженных градиентов, который позволяет проводить итерационный процесс минимизации в конечномерном пространстве  W_2, ρ¹ с нормой

                                      N-1
    || v ||  =  ( || v || Q2  +   ∑   ρj (vj - vj-1)2 )1/2 ,
                                      j=2
                       N-1
    || v || Q  =  (   ∑   qj vj2 ) ,
                       j=1 

q_j > 0 и  ρ_j ≥ 0 - заданные числа. Если априори известно, что производная искомого решения является достаточно гладкой функцией, то задание параметров  ρ_j > 0, связанных с конкретной нормировкой  W_2, ρ¹ (при  ρ_j ≡ 0 эта ноpма v  совпадает с || v || _Q) позволяет увеличить точность и улучшить качество приближений.

Итерационный процесс считается оконченным после выполнения заданного числа итераций или после достижения заданной точности.

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

    int ef02r_c (integer *m, integer *n, real *a, real *f,
            integer *ng, real *v1, real *v2, integer *mu, real *h, real *p,
            real *q, real *r, real *e, integer *isp, integer *lu, real *u,
            real *fu, integer *kn, real *b)

   Параметры 

   Версии:  нет 

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

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

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

int main(void)
{
    /* Initialized data */
    static float v1[11] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };
    static int ng = 3;
    static int lu = 1;
    static int isp = 20;
    static float v2[11] = { 1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,
            0.f };
    static int mu[11] = { 1,1,1,1,1,-1,-1,-1,-1,0,0 };
    static float h__[11] = { .1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,0.f };
    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 float q[11] = { .05f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.05f };
    static float r__[12] = { 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 m = 11;
    static int n = 11;

    /* Builtin functions */
    double sin(double);

    /* Local variables */
    extern int ef02r_c(int *, int *, float *, float *, int *, float *,
            float *, int *, float *, float *, float *, float *, float *,
            int *, int *, float *, float *, int *, float *);
    static float a[121]  /* was [11][11] */, b[90], e, f[11];
    static int i__, i, j;
    static float s, u[11], y, sf;
    static int kn;
    static float fu;
    int i__1, i__2;
    float r__1;

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

    e = 1e-7f;
    i__1 = m;
    for (i__ = 1; i__ <= i__1; ++i__) {
        i__2 = n;
        for (j = 1; j <= i__2; ++j) {
            if (i__ - j != 0) {
                goto l2;
            } else {
                goto l1;
            }
l1:
            a_ref(i__, j) = 1.f / q[i__ - 1];
            goto l3;
l2:
            a_ref(i__, j) = 0.f;
l3:
            ;
        }
    }
    u[0] = -.45833333333333331f;
    i__2 = n;
    for (j = 2; j <= i__2; ++j) {
/* l4: */
        u[j - 1] = 0.f;
    }
    sf = .025f;
    y = -.5f;
    i__2 = m;
    for (i__ = 1; i__ <= i__2; ++i__) {
        s = sf * (float)sin((float)(i__ * 123.456f + 789.f));
/* Computing 3rd power */
        r__1 = y;
        f[i__ - 1] = y - r__1 * (r__1 * r__1) / 3.f + s;
/* l5: */
        y += h__[i__ - 1];
    }
    ef02r_c(&m, &n, a, f, &ng, v1, v2, mu, h__, p, q, r__, &e, &isp, &lu, u,
            &fu, &kn, b);

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


Результат:

       u  =  ( -4.58e-1,   -4.11e-1,   -3.16e-1,   -2.21e-1,   -1.26e-1, 
                  -3.13e-2,    7.87e-2,    1.73e-1,    2.68e-1,    3.64e-1,   4.58e-1 )