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

Назначение

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

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

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

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

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

(1)     || K D u - f ||p2  =   inf   || K D v - f ||p2 ,
                                      v∈G 

где  K - известная матрица порядка  M * N с элементами  k_i j , которая в зависимости от ядра  k (x, y) и способа его аппроксимации может быть матрицей общего вида, либо симметричной, циркулянтной, теплицевой и т.п.,

f = ( f₁, ..., f_M ) - известный  M -  мерный вектор ,

u = ( u₁, ..., u_N ) - искомый  N - мерный вектор приближенного решения,

D - диагональная матрица коэффициентов квадратурной формулы  d_j ,  j = 1, ..., N,

                        M 
       || Z ||p2  =  ∑  pi zi2 ,      pi > 0  -  весовые коэффициенты.
                       i=1 

Линейные ограничения  G представляют собой одно из множеств  G₁,  G₂ или  G₁ ∩ G₂, где

 G1 = { v = (v1, ..., vN) :  aj ≤ vj ≤ bj ,  j = 1, ..., N } ,
 G2 = { v = (v1, ..., vN) :  mj (vj + 1 - vj) ≥ 0 ,  j = 1, ..., N - 1 } , 

a = (a₁, ..., a_N),  b = (b₁, ..., b_N) - заданные  N - мерные векторы нижних и верхних ограничений в  G₁; параметры  m_j определяются условиями кусочной монотонности и принимают значения 1, 0 или -1.

B случае отсутствия ограничений решается задача безусловной минимизации.

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

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

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

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

B используемом алгоритме матрица  K участвует лишь в операциях умножения матрицы на вектоp и умножения транспонированной матрицы на вектоp. Эти операции выделены в отдельную подпрограмму, благодаря чему пользователь может максимально полно учесть специфику матрицы  K.

B самой подпрограмме ef01r_c матрица  K не используется, поэтому способ ее представления не регламентируется.

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

    int ef01r_c (integer *m, integer *n, real *a, real *f,
            integer *ng, real *au, real *bu, integer *mu, real *p, real *q,
            real *d, real *r, real *e, integer *isp, real *u, real *fu,
            integer *kn, real *b, S_fp op)

   Параметры 

                n
      wi  =  ∑   ki j vj ,    i = 1, ..., m ,
               j=1 

                m
      wj  =  ∑   kj i vi ,    j = 1, ..., n .
               i=1 

   Версии:  нет 

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

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

                             xi+1    yj+1
      ki j  =  1/h1x       ∫       ∫     k (x, y)  dx dy ,   dj ≡ 1 .
                              xi      yj 

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

Матрица  K имеет общий вид

int main(void)
{
    /* Initialized data */
    static float a[9] /* was [3][3] */ = { 3.f,1.f,5.f,1.f,-3.f,2.f,1.f,2.f,
            1.f };
    static float bu[3] = { 10.f,10.f,10.f };
    static float f[3] = { 4.f,3.f,4.f };
    static int mu[3] = { 1,1,0 };
    static float p[3] = { 1.f,1.f,1.f };
    static float q[3] = { 1.f,1.f,1.f };
    static float d__[3] = { 1.f,1.f,1.f };
    static float r__[4] = { 0.f,0.f,0.f,0.f };
    static float u[3] = { 0.f,0.f,0.f };
    static float au[3] = { -10.f,-10.f,-10.f };

    /* Local variables */
    extern int ef01r_c(int *, int *, float *, float *, int *, float *,
            float *, int *, float *, float *, float *, float *, float *,
            int *, float *, float *, int *, float *, U_fp);
    extern int ef01r1_c();
    static float b[26], e;
    static int m, n, ng, kn;
    static float fu;
    static int isp;

    m = 3;
    n = 3;
    ng = 2;
    e = 1e-7f;
    isp = 20;
l1:
    ef01r_c(&m, &n, a, f, &ng, au, bu, mu, p, q, d__, r__, &e, &isp, u, &fu,
            &kn, b, (U_fp)ef01r1_c);

        printf("\n %13.2e %13.2e %13.2e \n", u[0], u[1], u[2]);
        printf("\n %17.6e \n", fu);
    if (kn != 0) { goto l2;}
       else { goto l1;}
l2:
    return 0;
} /* main */


Результат:    u  =  (-1., 2., 5.)