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

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

Назначение

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

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

qs12r_c вычисляет интеграл с заданной абсолютной погрешностью на отрезке интегрирования по формуле 

               B
               ∫ f (x) dx ≈
              A
                                          N
                          ≈ (B-A)   ∑ Ci f ( (B-A) yi / 2 + (A+B) / 2 ) / 2 ,
                                         i=1

где yi и Ci соответственно узлы и веса Гаусса для отрезка [ - 1, 1].

Л.Г.Васильева. Набор стандартных программ численного интегрирования с фиксированным распределением узлов. Сб."Численный анализ на ФОРТРАНе", вып. 8, Изд-во МГУ, 1974.

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

    int qs12r_c (real *rint, real *a, real *b, R_fp f, real *e,
            integer *n, integer *ierr)

   Параметры
rint - вещественная переменная, содержащая вычисленное значение интеграла;
a, b - заданные нижний и верхний пределы интегрирования (тип: вещественный);
f - имя вещественной подпрограммы - функции, вычисляющей подинтегральную функцию f (x);
e - заданная абсолютная погрешность вычисления интеграла (тип: вещественный);
n - целая переменная, задающая число узлов интегрирования;
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
ierr = 65 - когда заданная точность не может быть достигнута при максимально возможном числе узлов интегрирования; значение n полагается равным - 1;
ierr = 66 - когда заданное число узлов интегрирования не принадлежит списку возможных узлов. 
   Версии:  нет 
   Вызываемые подпрограммы
utqs10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qs12r_c. 
   Замечания по использованию
 

B процессе вычисления число узлов интегрирования n может принимать только следующие значения:

(1)       2, 4, 6, 8, 10, 12, 16, 24, 32, 48, 64, 96.

Проводится несколько последовательных просчетов, причем для каждого последующего просчета берется очередное число узлов из списка (1).

Если абсолютная величина разности двух последовательных просчетов интеграла не превосходит e, то счет заканчивается, и за значение интеграла принимается результат последнего просчета. Если список (1) исчерпан и заданная точность не достигнута, то n полагается равным - 1.

При e ≤ 0 происходит только один просчет с заданным n.

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

int main(void)
{
    /* Local variables */
    static int ierr;
    extern int qs12r_c(float *, float *, float *, R_fp, float *, int *, int *);
    static float rint, a, b, e;
    extern float f_c();
    static int n;

    a = 1.f;
    b = 20.085536923186663f;
    e = 1e-5f;
    n = 2;
    qs12r_c(&rint, &a, &b, (R_fp)f_c, &e, &n, &ierr);

    printf("\n %16.7e %5i \n",rint,n);
    return 0;
} /* main */

float f_c(float *x)
{
    /* System generated locals */
    float ret_val;

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

    ret_val = 1.f / (*x * (float)sqrt((float)((float)log(*x) + 1.f)));
    return ret_val;
} /* f_c */


Результаты:

      rint  =  1.9999999999 ;
      n  =  32