Текст подпрограммы и версий
qsf6r_c.zip  qsf6d_c.zip  		Тексты тестовых примеров
tqsf6r_c.zip  tqsf6d_c.zip 

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

Назначение 

 Вычисление интегралов вида
                  B
(1)             ∫  f (x) sin( ρ x + φ ) dx ,
                A
                  B
(2)             ∫  f (x) cos( ρ x + φ ) dx
                A

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

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

qsf6r_c вычисляет интегралы (1) и (2) для A и B таких, что | B - A | ≥ 10 соответственно с погрешностями 

                     B
      E ( 1 + |  ∫  f (x) sin( ρ x + φ) dx | )   и
                   A     
                     B                               
      E ( 1 + |  ∫  f (x) cos( ρ x + φ) dx | ) ,
                   A
где E задается пользователем.

Если нижний (верхний) предел интегрирования pавен бесконечности, то его следует задавать близким к минимальному (максимальному) числу, представимому в машине.

Метод вычисления интеграла по большому отрезку состоит в том, что интеграл последовательно вычисляется по частичным отрезкам, таким, что длина каждого следующего в 2 раза больше предыдущего. Счет интегралов прекращается, если значение интеграла на очередном частичном отрезке становится достаточно малым.

Hа каждом частичном отрезке интеграл вычисляется с автоматическим выбором шага по формулам интерполяционного типа, точным для многочленов 5 - ой степени с весом  exp i (ρ x + φ) .

За приближенное значение интеграла на всем отрезке принимается сумма вычисленных интегралов по всем частичным отрезкам.

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

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

    int qsf6r_c (real *rint1, real *rint2, real *a, real *b,
        R_fp f, real *ro, real *fi, real *e, integer *l1, integer *l2,
            real *xa)

   Параметры
rint1 -
rint2   		вещественные переменные, содержащие вычисленные значения интегралов (1) и (2);
a, b - заданные нижний и верхний пределы интегрирования (тип: вещественный);
f - имя вещественной подпрограммы - функции, вычисляющей подинтегральную функцию  f (x);
ro - заданное значение параметра  ρ (тип: вещественный);
fi - заданное значение параметра  φ (тип: вещественный);
e - заданная точность вычисления интеграла (тип: вещественный);
l1, l2 - задают режим работы подпрограммы (тип: целый); при этом, если:
 

l1 = 1 и l2 = 0, то вычисляется интеграл (1);

l1 = 0 и l2 = 1, то вычисляется интеграл (2);

l1 = 1 и l2 = 1, то вычисляются интегралы (1) и (2);
xa - вещественная переменная, значение которой на выходе из подпрограммы pавно либо b, если заданная точность достигнута, либо значению точки "особенности" подинтегральной функции. 
   Версии
qsf6d_c - вычисление с удвоенной точностью интегралов вида (1) и (2) для больших отрезков интегрирования по формулам интерполяционного типа 5 - ой степени с гарантированной точностью. 
   Вызываемые подпрограммы
utqs11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qsf6r_c;
utqs13_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qsf6d_c. 
   Замечания по использованию
 

Могут встретиться случаи, когда подинтегральная функция имеет особенность, в окрестности которой шаг интегрирования становится недопустимо малым. Тогда значение xa полагается равным координате такой "особенности", значения rint1 и rint2 не определены, а utqs11_c и utqs13_c выдают диагностическое сообщение: "заданная точность не может быть достигнута" и печатают значение xa.

В подпрограмме qsf6d_c параметры rint1, rint2, a, b, f, ro, fi, e, xa имеют тип double.

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

int main(void)
{
    /* Local variables */
    static float a, b;
    extern int qsf6r_c(float *, float *, float *, float *, R_fp,
                       float *, float *, float *, int *, int *, float *);
    static float e;
    extern float f_c();
    static float s, i1, i2;
    static int l1, l2;
    static float fi, xa, ro;

    e = 1e-6f;
    a = 0.f;
    b = 1e4f;
    s = 0.f;
    ro = .01f;
    fi = 0.f;
    l1 = 1;
    l2 = 0;
    qsf6r_c(&i1, &i2, &a, &b, (R_fp)f_c, &ro, &fi, &e, &l1, &l2, &xa);

    printf("\n %16.7e %16.7e \n",i1,i2);
    printf("\n %16.7e %16.7e \n",xa,b);
    return 0;
} /* main */

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

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

    if (*x <= 0.f) {
        goto l2;
    } else {
        goto l4;
    }
l2:
    ret_val = 0.f;
    goto l3;
l4:
    ret_val = 1.f / (float)sqrt(*x);
l3:
    return ret_val;
} /* f_c */


Результаты:

      i1  =  11.673418089
      ха  =  10000