Текст подпрограммы и версий
qsf9r_c.zip  qsf9d_c.zip  		Тексты тестовых примеров
tqsf9r_c.zip  tqsf9d_c.zip 

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

Назначение

Вычисление определенных интегралов 

                  B
(1)             ∫  F (x) sin( ρ x + φ ) dx ,
                A
                  B
(2)             ∫  F (x) cos( ρ x + φ ) dx
                A

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

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

Подпрограмма qsf9r_c вычисляет интегралы (1) и (2), или один из них от векторной функции F (x) = ( f1 (x), ..., fN (x) ) с покомпонентной погрешностью 

                     B
      E ( 1 + |  ∫  fi (x) sin( ρ x + φ) dx | )  или
                   A     
                     B                               
      E ( 1 + |  ∫  fi (x) cos( ρ x + φ) dx | )      
                   A

При ρ h ≥ 0.33 вычисление ведется по формулам интерполяционного типа пятой степени точности, а при ρ h < 0.33 по обобщенной квадратурной формуле Hьютона - Котеса шестой степени точности с шагом  h = (B - A)/M, где M - число частичных отрезков разбиения.

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

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

    int qsf9r_c (real *rint, real *a, real *b, S_fp f, real *ro,
            real *fi, integer *n, real *e, integer *k, integer *l1, integer *l2,
             integer *ind, real *r, integer *ierr)

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

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

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

l1 = 1 и l2 = 1, то вычисляются интегралы (1) и (2);
ind - целый вектоp длины n, каждая компонента котоpого pавна числу частичных отрезков разбиения, при котоpом достигается заданная точность по этой компоненте;
r - вещественный рабочий вектоp длины 3n, если вычисляются оба интеграла, и 2n, если вычисляется один из них;
ierr - целая переменная, служащая для диагностических сообщений:
ierr = 65 - если заданная точность не может быть достигнута при максимально возможном числе (1048576) частичных отрезков разбиения. 
   Версии
qsf9d_c - вычисление с удвоенной точностью определенных интегралов (1) и (2) от векторной функции по формулам интерполяционного типа пятой степени точности. 
   Вызываемые подпрограммы
utqs10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qsf9r_c.
utqs12_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qsf9d_c. 
   Замечания по использованию
  1. 

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

  2. 

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

  3. 

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

  4. 

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

  5. 

Погрешность определяется по модулю разности двух просчетов по m и 2m частичным отрезкам разбиения. Если заданная точность не достигнута, то значение m удваивается.

  6.  В подпрограмме qsf9d_c параметры rint, a, b, f, ro, fi, e, r имеют тип double.

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

int main(void)
{
    /* System generated locals */
    float r__1;

    /* Local variables */
    static int ierr;
    static float rint[10] /* was [5][2] */;
    extern int qsf9r_c(float *, float *, float *, U_fp, float *, float *,
                       int *, float *, int *, int *, int *, int *,
                       float *, int *);
    static int i__;
    static float r__[15], pi0;
    static int ind[5];
    extern int fun_c();

#define rint_ref(a_1,a_2) rint[(a_2)*5 + a_1 - 6]

    pi0 = 1.5707963268f;
    r__1 = -pi0;
    qsf9r_c(rint, &r__1, &pi0, (U_fp)fun_c, &c_b1, &c_b2, &c__5, &c_b4, &c__2,
            &c__1, &c__1, ind, r__, &ierr);

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


int fun_c(float *x, float *y, int *i__)
{
    static int j;

    /* Parameter adjustments */
    --y;

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


Результаты:

              rint_ref (i__, 1)                         rint_ref (i__, 2)               ind (i__)

      1.9999999618 e+00             -8.661 e-12                            2
     -2.676 e-11                             9.3480217805 e-01              2
      2.8044065499 e+00             -6.428 e-12                            2
     -5.268 e-11                             9.5850994403 e-01              2
      4.7925452374 e+00             -2.098 e-11                            4

      ierr  =  0

      i__  =   (1, 2, 3, 4, 5)