Текст подпрограммы и версий
qs81r_c.zip  qs81d_c.zip  		Тексты тестовых примеров
tqs81r_c.zip  tqs81d_c.zip 

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

Назначение

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

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

qs81r_c вычисляет N - кратный интеграл с заданной абсолютной погрешностью E по квадратурной формуле 

        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  )  =  I ( k1, k2, ..., kN ) .

Для программным образом выбираемой последовательности вектоpов Km = (k1 m, k2 m, ..., kN m) вычисляются интегральные суммы I (Km) до тех пор, пока не будет найден вектоp Km0 с компонентами  ki m0 < 128, такой, что 

   |  I ( k1m0, k2m0, ..., kNm0) - I ( k1m0 + 1 , k2m0 + 1 , ..., kNm0 + 1)  |  <  E ,
 где E - заданная величина.

Значение интегральной суммы  I ( k1m0, k2m0, ..., kNm0 )
берется в качестве приближенного значения интеграла. B случае если такого вектоpа нет, счет прекращается, значение интеграла неопределено.

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

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

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

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

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

Если заданная точность достигнута, то при выходе из подпрограммы j - ый элемент вектоpа np pавен числу узлов интегрирования по j - ой переменной; B противном случае хотя бы один из элементов вектоpа np будет pавен 129. Последнее указывает на недостаточность узлов интегрирования по соответствующей переменной, при этом значение rint не определено, и подпрограммы utqs11_c, utqs13_c выдают сообщение: "заданная точность не может быть достигнута".

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

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

int main(void)
{
    /* Local variables */
    static int ierr;
    static float rint;
    extern int qs81r_c(float *, int *, float *, float *,
                       R_fp, float *, int *, float *, float *, int *);
    static float a[5], b[5], e;
    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__) {
        a[i__ - 1] = -3.f;
/* l1: */
        b[i__ - 1] = 3.f;
    }
    e = .1f;
    qs81r_c(&rint, &n, a, b, (R_fp)f_c, &e, np, x, w, &ierr);

    printf("\n %16.6e %16.7e &16.7e \n",a[0],a[1],a[2]);
    printf("\n %16.7e %16.7e \n",a[3],a[4]);
    printf("\n %16.7e %16.7e &16.7e \n",b[0],b[1],b[2]);
    printf("\n %16.7e %16.7e \n",b[3],b[4]);
    printf("\n %16.7e \n",rint);
    printf("\n %5i %5i %5i %5i %5i \n",np[0],np[1],np[2],np[3],np[4]);
    printf("\n %5i \n",ierr);
    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.6291

      np    =   (10, 10, 10, 10, 10)
      ierr  =  0