Текст подпрограммы и версий
qsfpr_c.zip  qsfpd_c.zip 
Тексты тестовых примеров
tqsfpr_c.zip  tqsfpd_c.zip 

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

Назначение

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

            B
     I =  ∫  [ cos(tx) Re F(x) - sin(tx) Im F(x) ] dx ,
          A 

где Re F (x) и  Im F (x) - действительная и мнимая части комплексной функции F(x) действительного переменного  x, на основе квадратурной формулы Филона.

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

Подпрограмма qsfpr_c вычисляет интеграл I с заданной абсолютной или относительной точностью на основе квадратурной формулы Филона [1], в которой плотность подынтегральной функции F (x) вычисляется в девяти равноотстоящих точках каждого подинтервала.

При делении подинтервала узлы и значения функции F (x) запоминаются для последующего использования [2].

1.  Справочник по специальным функциям. Под ред. М.А.Абрамовица и И.Стиган. - M.: Hаука, 1979, 832 с.
2.  Форсайт Дж., Малькольм M., Моулер K. Машинные методы математических вычислений. - M.: Мир, 1980, 280 с.

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

    int qsfpr_c(C_fp fun, real *a, real *b, real *abserr,
        real *relerr, real *t, real *result, real *errest, integer *nofun,
        real *flag)

Параметры

fun - имя комплексной подпрограммы - функции, вычисляющей плотность F (x) подынтегральной функции (должна быть описана в основной программе в операторе external);
a, b - заданные нижний и верхний пределы интегрирования (тип: вещественный);
abserr - заданная абсолютная погрешность вычисления интеграла (тип: вещественный);
relerr - заданная относительная погрешность вычисления интеграла (тип: вещественный);
t - заданное значение параметра  t (тип: вещественный);
result - приближение к интегралу, удовлетворяющее менее жесткой из двух границ погрешности (тип: вещественный);
errest - оценка погрешности вычисленного значения интеграла (тип: вещественный);
nofun - число вычисленных значений плотности подынтегральной функции F (x), потребовавшихся для получения результата (тип: целый);
flag - индикатор надежности:
если flag = 0, то result, вероятно, удовлетворяет заданной погрешности,
если flag = xxx.yyy, то xxx - число интегралов, на которых не удалось достигнуть заданной точности, а 0.yyy = (b - x0)/(b - a) - граница отрезка [a, x0], на котоpом был израсходован лимит обращений к подпрограмме, вычисляющей функцию fun.

Версии

qsfpd_c - вычисление этих же интегралов в режиме удвоенной точности.

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

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

  1. 

Для подпрограммы qsfpd_c параметры a, b, t, er, ea, rint, flag, errest, а также параметр x, из подпрограммы  fun имеют тип double;

  2.  Для подпрограммы qsfpd_c функция fun должна быть описана как комплексная двойной точности (см. пример использования)

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

              π                                 π     
          (   ∫  sin10x F(x) dx    =   ∫  sin10x sin x/2 dx  )
              0                                0

int main(void)
{
    /* Local variables */
    static float flag__, rint, a, b, t;
    static int nofun;
    extern int qsfpr_c(C_fp, float *, float *, float *, float *,
                       float *, float *, float *, int *, float *);
    static float ea, er, pi, errest;
    extern /* Complex */ VOID fun_c();

    pi = 3.1415926f;
    er = 1e-5f;
    ea = 0.f;
    a = 0.f;
    t = 10.f;
    b = pi;
    qsfpr_c((C_fp)fun_c, &a, &b, &ea, &er, &t, &rint, &errest, &nofun, &flag__);

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

/* Complex */ VOID fun_c(complex * ret_val, float *x)
{
    /* System generated locals */
    float r__1;
    complex q__1, q__2;

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

    /* Local variables */
    static complex zi;

    zi.r = 0.f, zi.i = 1.f;
    q__2.r = -zi.r, q__2.i = -zi.i;
    r__1 = (float)sin((float)(*x / 2.f));
    q__1.r = r__1 * q__2.r, q__1.i = r__1 * q__2.i;
    ret_val->r = q__1.r,  ret_val->i = q__1.i;
    return ;
} /* fun_c */


Результаты:

      rint  =  -1.0025063 * 10-1
      
      nofun  =  417 
      flag__  =  0
      errest  =  2.3 * 10-10