Текст подпрограммы и версий
qs15r_c.zip  qs15d_c.zip  		Тексты тестовых примеров
tqs15r_c.zip  tqs15d_c.zip 

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

Назначение

Вычисление определенного интеграла от векторной функции по формулам Гаусса.

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

Подпрограмма qs15r_c вычисляет интеграл
           B
           ∫ F (x) dx
          A
от векторной функции   F (x) = ( f1 (x), ... , fN (x) )
с покомпонентной погрешностью
                         B
      ЕРS ( 1 + |  ∫  fi (x) dx | )
                        A

по квадратурным формулам Гаусса с разбиением основного отрезка интегрирования на частичные.

Точность может достигаться как за счет удвоения числа частичных отрезков разбиения при фиксированном числе узлов квадратурной формулы, применяемой на каждом частичном отрезке, так и за счет увеличения числа узлов квадратурной формулы при фиксированном числе отрезков разбиения.

Я.М.Жилейкин, М.В.Соколовский. Набор стандартных программ для вычисления интегралов от векторных функций, Сб. "Методы и алгоритмы в численном анализе", Изд-во МГУ, 1981 г.

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

    int qs15r_c (real *rint, real *a, real *b, S_fp f, integer *n,
             real *eps, integer *k, integer *ng, integer *l, integer *ind,
             real *r, integer *ierr)

   Параметры
rint - вещественный вектоp длины n, содержащий вычисленные значения интегралов;
a, b - заданные нижний и верхний пределы интегрирования (тип: вещественный);
f - имя подпрограммы, вычисляющей подинтегральные функции;
n - размерность вектоpа - функции F (x) (тип: целый);
eps - заданная меpа погрешности вычисления интеграла (тип: вещественный);
k - целая переменная, задающая начальное число частичных отрезков разбиения;
ng - целая переменная, задающая начальное число узлов квадратурной формулы;
l - задает режим работы подпрограммы (тип: целый); при этом, если:
l = 0 - то заданная точность достигается за счет увеличения числа частичных отрезков разбиения при заданном ng;
l ≥ 1 - то заданная точность достигается за счет увеличения числа узлов квадратурной формулы при заданном k;
ind - целый вектоp длины n, каждая компонента которого pавна при:
l = 0 - числу частичных отрезков разбиения, при котоpом достигается заданная точность по этой компоненте;
l ≥ 1 - числу узлов квадpатуной фоpмулы, пpи котоpом достигается заданная точность по этой компоненте;
r - вещественный рабочий вектоp длины 2n;
ierr - целая переменная, служащая для диагностических сообщений:
ierr = 65 - когда заданная точность не может быть достигнута при максимально возможном числе узлов квадратурной формулы (96);
ierr = 66 - если заданное число узлов квадратурной формулы не является допустимым;
ierr = 67 - если заданная точность не может быть достигнута при максимально возможном числе (1048576) частичных отрезков разбиения. 
   Версии
qs15d_c - вычисление с удвоенной точностью определенного интеграла от векторной функции по формулам Гаусса. 
   Вызываемые подпрограммы
utqs11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qs15r_c.
utqs13_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qs15d_c. 
   Замечания по использованию
  1. 

Число узлов квадратурной формулы может принимать лишь следующие допустимые значения 2, 4, 6, 8, 10, 12, 16, 24, 32, 48, 64, 96.

  2. 

Первый оператор подпрограммы f должен иметь вид:
int f(float *x, float *y, int *i)
Здесь:
x - аргумент функции f (тип: вещественный);
y - вещественный вектоp длины n вычисленных значений функции f;
i - целый вектоp длины n; если его i - ая компонента отлична от нуля, то интеграл от этой компоненты функции f не вычисляется.

  3. 

При eps ≤ 0 происходит только один просчет при заданных значениях k и ng.

  4. 

Вычисление интеграла от отдельной компоненты прекращается как только по этой компоненте достигается заданная точность.

  5. 

Если для каких - либо компонент векторной функции заданная погрешность интегрирования не может быть достигнута, то соответствующие компоненты вектоpа ind полагаются равными нулю.

  6. 

Погрешность определяется по модулю разности двух просчетов по m и 2*m частичным отрезкам разбиения при l = 0 либо по двум просчетам с различным, но следующим друг за другом, согласно указанному выше списку, числом узлов квадратурных формул при l ≥ 1.

  7.  В подпрограмме qs15d_c параметры rint, a, b, f, e, eps, r имеют тип double.

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

int main(void)
{
    /* Local variables */
    static int ierr;
    extern int qs15r_c(float *, float *, float *, U_fp, int *, float *,
                       int *, int *, int *, int *, float *, int *);
    static float rint1[5], rint2[5];
    static int i__;
    static float r__[10];
    extern int fun_c();
    static int ind1[5], ind2[5];

    qs15r_c(rint1, &c_b1, &c_b2, (U_fp)fun_c, &c__5, &c_b4, &c__2, &c__4,
            &c__0, ind1, r__, &ierr);
    qs15r_c(rint2, &c_b1, &c_b2, (U_fp)fun_c, &c__5, &c_b4, &c__1, &c__2,
            &c__1, ind2, r__, &ierr);

    for (i__ = 1; i__ <= 5; ++i__) {
        printf("\n %16.7e %2i %16.7e %2i \n",
        rint1[i__-1],ind1[i__-1],rint2[i__-1],ind2[i__-1]);
/* l1: */
    }
    return 0;
} /* main */

int fun_c(float *x, float *y, int *i__)
{
    /* System generated locals */
    float r__1;

    /* Local variables */
    static int j;

    /* Parameter adjustments */
    --y;

    /* Function Body */
/* Computing 6th power */
    r__1 = *x, r__1 *= r__1;
    y[1] = r__1 * (r__1 * r__1);
    for (j = 1; j <= 4; ++j) {
/* l1: */
        y[j + 1] = y[j] * *x;
    }
    return 0;
} /* fun_c */


Результаты:

               rint1                      ind1

      1.42857142855-01               2
      1.24999999998-01               2
      1.11111111108-01                 4
      9.99999999985-01               8
      9.09090909076-02               8

               rint2                      ind2

      1.42857142855-01               4
      1.24999999998-01               4
      1.11111111109-01                 6
      9.99999999985-01               6
      9.09090909078-02               6

      ierr  =  0