Текст подпрограммы и версий
qs82r_c.zip  qs82d_c.zip  		Тексты тестовых примеров
tqs82r_c.zip  tqs82d_c.zip 

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

Назначение

Вычисление определенного N - кратного (2 ≤ N ≤ 15) интеграла по прямоугольному параллелепипеду методом Гаусса.

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

qs82r_c вычисляет N - кратный интеграл по квадратурной формуле 

        b1 b2    bN
       ∫   ∫  ...  ∫  f (x1, x2, ... , xN) dx1 dx2 ... dxN  ≈
     a1 a2    aN
       N                      k1      k2       kN
 ≈ ( ∏ ( bi - ai )/2 )  ∑      ∑  ...  ∑  cn1 cn2 ... cnN *
      i=1                    n1=1  n2=1   nN=1
 * f (  ( b1 - a1 )/2 xn1 + ( a1 + b1 )/2 ,  ( b2 - a2 )/2 xn2 + ( a2 + b2 )/2 ,
         ... , ( bN - aN )/2 xnN + ( aN + bN )/2  ) ,

где 2 ≤ ki ≤ 128, а xni и cni ( ni = 1, 2, ..., ki ) соответственно узлы и веса квадратурной формулы Гаусса для отрезка [-1, 1].

Я.М.Жилейкин, А.Г.Симакин. Набор стандартных программ приближенного вычисления многократных интегралов с помощью квадратурных формул Гаусса. Сб. "Численный анализ на ФОРТРАНе", вып. 19, Изд-во МГУ, 1977.

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

    int qs82r_c (real *rint, integer *n, real *a, real *b, R_fp f,
            integer *np, real *x, real *w)

   Параметры
rint - вещественная переменная, содержащая вычисленное значение интеграла;
n - заданная кратность интегрирования (тип: целый);
a, b - вещественные векторы длины n, задающие соответвенно нижние и верхние пределы интегрирования;
f - имя вещественной подпрограммы - функции, вычисляющей подинтегральную функцию f (x);
np - целый вектоp длины n, задающий число узлов интегрирования по каждой переменной:  np (i) = ki,  i = 1, 2, ..., n;
x, w - вещественные двумерные массивы размера n на 64, используемые в подпрограмме как рабочие; 
   Версии
qs82d_c - вычисление с удвоенной точностью определенного N - кратного (2 ≤ N ≤ 15) интеграла по прямоугольному параллелепипеду методом Гаусса. 
   Вызываемые подпрограммы
qs80r_c - подпрограмма вычисления узлов и весов Гаусса - Лежандра.
qs80d_c - подпрограмма вычисления с удвоенной точностью узлов и весов Гаусса - Лежандра. 
   Замечания по использованию
 

Компоненты вектоpа np в пределах от 2 до 128.

Первый оператор подпрограммы - функции f должен иметь вид
      float f_c(float *x) ,
где x - вещественный вектоp длины n, являющийся аргументом подинтегральной функции.

В подпрограмме qs82d_c параметры rint, a, b, f, x, w имеют тип double.

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

int main(void)
{
    /* Local variables */
    static float rint;
    extern int qs82r_c(float *, int *, float *, float *,
                       R_fp, int *, float *, float *);
    static float a[5], b[5];
    extern float f_c();
    static int i__, n;
    int i__1;
    static float w[320]  /* was [5][64] */,
                 x[320]  /* was [5][64] */;
    static int np[5];

    n = 5;
    i__1 = n;
    for (i__ = 1; i__ <= i__1; ++i__) {
        np[i__ - 1] = 10;
        a[i__ - 1] = -3.f;
/* l1: */
        b[i__ - 1] = 3.f;
    }
    qs82r_c(&rint, &n, a, b, (R_fp)f_c, np, x, w);

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

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

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

    /* Local variables */
    static int i__;
    static float r__, s;

    /* Parameter adjustments */
    --x;

    /* Function Body */
    s = 0.f;
    for (i__ = 1; i__ <= 5; ++i__) {
        r__ = x[i__];
/* l1: */
        s += r__ * r__;
    }
    ret_val = (float)exp((float)(s * -.5f));
    return ret_val;
} /* f_c */


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