Текст подпрограммы и версий
qs81r_p.zip  qs81e_p.zip  		Тексты тестовых примеров
tqs81r_p.zip  tqs81e_p.zip 

Подпрограмма:  QS81R (модуль QS81R_p)

Назначение

Вычисление определенного 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.

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

procedure QS81R(var RINT :Real; N :Integer; var A :Array of Real;
                var B :Array of Real; F :Func_F1A; E :Real;
                var NP :Array of Integer; var X :Array of Real;
                var W :Array of Real; var IERR :Integer);


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

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

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

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

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

Unit TQS81R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FQS81R_p, QS81R_p;

function TQS81R: String;

implementation

function TQS81R: String;
var
N,I,_i,IERR :Integer;
E,RINT :Real;
A :Array [0..4] of Real;
B :Array [0..4] of Real;
NP :Array [0..4] of Integer;
X :Array [0..319] of Real;
W :Array [0..319] of Real;
label
_1;
begin
Result := '';  { результат функции }
N := 5;
FОR I:=1 to N do
 begin
  А[I-1] := -3.0;
_1:
  B[I-1] := 3.0;
 end;
E := 1.E-1;
QS81R(RINT,N,A,B,FQS81R,E,NP,X,W,IERR);
Result := Result + #$0D#$0A;
for _i:=0 to 4 do
 begin
  Result := Result + Format('%20.16f ',[A[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 4 do
 begin
  Result := Result + Format('%20.16f ',[B[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f ',[RINT]) + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 4 do
 begin
  Result := Result + Format('%10d ',[NP[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('  %5d ',[IERR]) + #$0D#$0A;
UtRes('TQS81R',Result);  { вывод результатов в файл TQS81R.res }
exit;
end;

end.

Unit FQS81R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

function FQS81R( X :Array of Real): Real;

implementation

function FQS81R( X :Array of Real): Real;
var
I :Integer;
S,R :Real;
label
_1;
begin
{ Result - прототип имени функции FQS81R на FORTRANe }
S := 0.0;
for I:=1 to 5 do
 begin
  R := X[I-1];
_1:
  S := S+R*R;
 end;
Result := Exp(-0.5*S);
exit;
end;

end.

Результаты:

      RINT  =  97.6291

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