Текст подпрограммы и версий ( Фортран )
qs81r.zip  qs81d.zip  		Тексты тестовых примеров ( Фортран )
tqs81r.zip  tqs81d.zip 
Текст подпрограммы и версий ( Си )
qs81r_c.zip  qs81d_c.zip  		Тексты тестовых примеров ( Си )
tqs81r_c.zip  tqs81d_c.zip 
Текст подпрограммы и версий ( Паскаль )
qs81r_p.zip  qs81e_p.zip  		Тексты тестовых примеров ( Паскаль )
tqs81r_p.zip  tqs81e_p.zip 

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

Назначение

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

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

QS81R вычисляет 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.

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

    SUBROUTINE  QS81R (RINT, N, A, B, F, E, NP, X, W, IERR) 
   Параметры
RINT - вещественная переменная, содержащая вычисленное значение интеграла;
N - заданная кратность интегрирования (тип: целый);
A, B - вещественные векторы длины N, задающие соответственно нижние и верхние пределы интегрирования;
F - имя вещественной подпрограммы - функции, вычисляющей подинтегральную функцию f (x);
E - заданная абсолютная погрешность вычисления интеграла (тип: вещественный);
NP - целый вектоp длины N, используемый в подпрограмме как рабочий;
X, W - вещественные двумерные массивы размера N на 64, используемые в подпрограмме как рабочие;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
IERR=65 - когда заданная точность не может быть достигнута по крайней меpе по одному направлению при максимально возможном числе узлов. 
   Версии
QS81D - вычисление с удвоенной точностью определенного N - кратного (2 ≤ N ≤ 15) интеграла по прямоугольному параллепипеду методом Гаусса c заданной абсолютной погрешностью. 
   Вызываемые подпрограммы
QS80R - подпрограмма вычисления узлов и весов Гаусса - Лежандра.
QS80D - подпрограмма вычисления с удвоенной точностью узлов и весов Гаусса - Лежандра.
UTQS11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QS81R.
UTQS13 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QS81D. 
   Замечания по использованию
 

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

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

  В подпрограмме QS81D параметры RINT, A, B, F, E, X, W имеют тип DOUBLE PRECISION.

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

      FUNCTION  F(X)
      DIMENSION  X(5)
      S = 0.
      DO 1 I = 1, 5
      R = X(I)
   1 S = S + R*R
      F = EXP(-0.5*S)
      RETURN
      END

      EXTERNAL  F
      DIMENSION  A(5), B(5), NP(5), X(5, 64), W(5, 64)
      E = 0.001
      N = 5
      DO 1 I = 1, N
      A (I) = -3.
   1 B (I) = 3.
      CALL  QS81R (RINT, N, A, B, F, E, NP, X, W, IERR)

Результаты:

      RINT  =  97.6291

      NP    =   (10, 10, 10, 10, 10)
      IERR  =  0