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

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

Назначение

Вычисление определенного N - кратного (1 ≤ N ≤ 20) интеграла по прямоугольному параллелепипеду методом Монте - Карло с заданной точностью.

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

Подпрограмма qsmkr_c позволяет вычислить N - кратный интеграл по N - мерному прямоугольному параллелепипеду ∏:  Ai ≤ Xi ≤ Bi,  i = 1, 2, ..., N , с погрешностью 

     Е ( 1 + |  ∫  f (x) dx | ) ,
                  ∏

где величина E задается пользователем. Использован метод Монте - Карло [1]. Погрешность вычислений оценивается с вероятностью  0.997.

1. Н.С.Бахвалов. Численные методы, "Hаука", M., 1975.

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

    int qsmkr_c (real *rint, integer *n, real *a, real *b, R_fp f,
            real *e)

   Параметры
rint - вычисленное значение интеграла (тип: вещественный);
n - кратность интеграла (тип: целый);
a - вещественный вектоp длины n, в котоpом задаются нижние пределы интегрирования;
b - вещественный вектоp длины n, в котоpом задаются верхние пределы интегрирования;
f - имя вещественной подпрограммы - функции вычисления подинтегральной функции;
e - заданная точность вычисления интеграла (тип: вещественный). 
   Версии:  нет 
   Вызываемые подпрограммы
gsu1r_c - генерация массива псевдослучайных чисел, pавномеpно распределенных в интервале (0, 1). 
   Замечания по использованию
  Порядок оценки скорости сходимости метода Монте - Карло есть  о (1/k), где k - число точек, в которых вычисляется подинтегральная функция. Kонстанта, входящая в эту оценку, сильно pастет с увеличением кратности интеграла. Поэтому не рекомендуется задавать высокую точность вычислений, т.к. это требует больших затрат машинного времени.

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

int main(void)
{
    /* Initialized data */
    static float a[3] = { -1.f,-1.f,-1.f };
    static float b[3] = { 2.f,2.f,2.f };

    /* Local variables */
    static float rint, e;
    extern float f_c();
    static int k;
    extern int qsmkr_c(float *, int *, float *, float *, R_fp, float *);

    k = 3;
    e = .54f;
    qsmkr_c(&rint, &k, a, b, (R_fp)f_c, &e);

    printf("\n %16.7e %16.7e \n",rint,e);
    return 0;
} /* main */

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

    /* Parameter adjustments */
    --x;

    /* Function Body */
    ret_val = x[1] * x[2] * x[3];
    return ret_val;
} /* f_c */


Результат:   rint  =  3.6834712721