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

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

Назначение

Вычисление определеного двухкратного интеграла по прямоугольнику.

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

qs40r_c вычисляет s - кратный (s = 2) интеграл по s - меpному прямоугольному параллепипеду с заданной абслютной погрешностью. Заменой переменных данный s - мерный прямоугольный параллепипед переводится в единичный куб G, и интеграл вычисляется по этому кубу. B основу метода положены кубатурные формулы различной точности: 

     Q1(s), Q2(s), Q3(s) .

Куб G разбивается на 2s равных кубов Gk (k = 1, 2, ..., 2s). По каждому из Gk от подинтегральной функции вычисляются Q1 (s) и Q2 (s) и проверяется условие Is*| Q2 (s) - Q1 (s) | > E , где Is - якобиан исходной замены переменных, E - абсолютная погрешность вычислений. Если это условие не выполнено, то за приближенное значение интеграла по Gk принимается Q3(s), если же для Gk это условие выполняется, то весь процесс снова повторяется для этого Gk.

Если при дроблении кубов наступит момент, когда об'ем какого - либо куба будет машинным нулем, то интеграл для этого куба полагается равным нулю и вычисления продолжаются для следующего куба Gk.

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

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

    int qs40r_c (real *rint, real *a1, real *a2, R_fp f,
            real *eps, integer *ierr)

   Параметры
rint - вычисленное значение интеграла (тип: вещественный);
a1 - вещественный вектоp длины 2, в котоpом задаются нижние пределы интегрирования;
a2 - вещественный вектоp длины 2, в котоpом задаются верхние пределы интегрирования;
f - имя вещественной подпрограммы - функции вычисления подинтегральной функции;
eps - заданная абсолютная погрешность вычисления интеграла (тип: вещественный);
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы (см. "Замечания по использованию"). 
   Версии:  нет 
   Вызываемые подпрограммы
utqs10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qs40r_c. 
   Замечания по использованию
 

Пользователь должен задавать значение eps абсолютной погрешности вычисления интеграла на шаге. Для гладких функций полная погрешность приближенного значения интеграла будет  o (eps). В случае не гладких или очень хоpоших функций полная погрешность может значительно отличаться от величины eps. Для вычисления истинной погрешности можно просчитать интеграл для различных значений eps.

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

При работе подпрограммы может произойти следующее событие (назовем это событие V): при дроблении кубов об'ем какого - либо куба будет машинным нулем. в этом случае работа программы не прекращается: интеграл для этого куба полагается равным нулю, вычисления продолжаются для следующего куба и переменной ierr присваивается значение 65. Если же при работе программы событие V не произошло ни разу, то переменной ierr присваивается значение нуль. Таким образом, знание на выходе программы значения ierr дает возможность пользователю контpолировать точность вычисления интеграла.

Возможны случаи, когда аргумент функции за счет погрешности округления выйдет за границу области интегрирования. Поэтому, иногда, разумно доопределить функцию в окрестности граничной точки нулем, при этом достаточно взять радиус окрестности равным 10 - 9.

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

int main(void)
{
    /* Local variables */
    static int ierr;
    extern int qs40r_c(float *, float *, float *, R_fp, float *, int *);
    extern float f_c();
    static float a1[2], a2[2], eps, int__;

    a1[0] = 0.f;
    a1[1] = 0.f;
    a2[0] = 1.f;
    a2[1] = 2.f;
    eps = .001f;
    qs40r_c(&int__, a1, a2, (R_fp)f_c, &eps, &ierr);

    printf("\n %16.7e %16.7e \n",a1[0],a1[1]);
    printf("\n %16.7e %16.7e \n",a2[0],a2[1]);
    printf("\n %16.7e \n",int__);
    printf("\n %5i \n",ierr);
    return 0;
} /* main */

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

    /* Parameter adjustments */
    --x;

    /* Function Body */
    ret_val = x[1] * 12 * x[2] / (x[1] * x[1] + 1) * (x[1] * x[1] + 1) *
             (x[2] * x[2] + 1) * (x[2] * x[2] + 1);
    return ret_val;
} /* f_c */


Результаты:   int__  =  124. ,   ierr  =  0