БЧА НИВЦ МГУ. QSFPR_C. Вычисление интегралов специального вида

Документ взят из кэша поисковой машины. Адрес оригинального документа : http://num-anal.srcc.msu.su/lib_na/cat/q_htm_c/qsfpr_c.htm
Дата изменения: Tue Apr 28 08:25:55 2015
Дата индексирования: Sun Apr 10 00:58:41 2016
Кодировка: Windows-1251

Текст подпрограммы и версий
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 - x₀)/(b - a) - граница отрезка [a, x₀], на кото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  )
              ⁰                                ⁰

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