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

Подпрограмма:  QS40R (модуль QS40R_p)

Назначение

Вычисление определеного двухкратного интеграла по прямоугольнику.

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

QS40R вычисляет s - кратный (s = 2) интеграл по s - меpному прямоугольному параллепипеду с заданной абслютной погрешностью. Заменой переменных данный s - мерный прямоугольный параллепипед переводится в единичный куб G, и интеграл вычисляется по этому кубу. B основу метода положены кубатурные формулы различной точности: 

     Q1(s), Q2(s), Q3(s) .

Куб G разбивается на 2s равных кубов Gk (k = 1, 2, ..., 2s). По каждому из Gk от подинтегральной функции вычисляются Q1 (s) и Q2 (s) и проверяется условие Is*| Q2 (s) - Q1 (s) | > E , где Is - якобиан исходной замены переменных, E - абсолютная погрешность вычислений. Если это условие не выполнено, то за приближенное значение интеграла по Gk принимается Q3(s), если же для Gk это условие выполняется, то весь процесс снова повторяется для этого Gk.

Если при дроблении кубов наступит момент, когда об'ем какого - либо куба будет машинным нулем, то интеграл для этого куба полагается равным нулю и вычисления продолжаются для следующего куба Gk.

Е.А.Лапшин. Набор стандартных программ приближенного вычисления двух - и тpех - кратных интегралов. Сб."Численный анализ на ФОРТРАНе", вып. 8, Изд-во МГУ, 1974.

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

procedure QS40R(var RINT :Real; var A1 :Array of Real;
                var A2 :Array of Real; F :Func_F1A; EPS :Real;
                var IERR :Integer);


   Параметры
RINT - вычисленное значение интеграла (тип: вещественный);
A1 - вещественный вектоp длины 2, в котоpом задаются нижние пределы интегрирования;
A2 - вещественный вектоp длины 2, в котоpом задаются верхние пределы интегрирования;
F - имя вещественной подпрограммы - функции вычисления подинтегральной функции;
EPS - заданная абсолютная погрешность вычисления интеграла (тип: вещественный);
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы (см. "Замечания по использованию"). 
   Версии:  нет 
   Вызываемые подпрограммы
UTQS10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QS40R. 
   Замечания по использованию
 

Пользователь должен задавать значение EPS абсолютной погрешности вычисления интеграла на шаге. Для гладких функций полная погрешность приближенного значения интеграла будет  o (EPS). В случае не гладких или очень хоpоших функций полная погрешность может значительно отличаться от величины EPS. Для вычисления истинной погрешности можно просчитать интеграл для различных значений EPS.

Подпрограмма ориентирована на вычисление интегралов от функций, определенных всюду в области интегрирования за исключением конечного числа точек. В точках, где функция не определена (например, обращается в бесконечность) ее необходимо доопределить каким - либо значением (например, положить равной нулю).

При работе подпрограммы может произойти следующее событие (назовем это событие V): при дроблении кубов об'ем какого - либо куба будет машинным нулем. в этом случае работа программы не прекращается: интеграл для этого куба полагается равным нулю, вычисления продолжаются для следующего куба и переменной IERR присваивается значение 65. Если же при работе программы событие V не произошло ни разу, то переменной IERR присваивается значение нуль. Таким образом, знание на выходе программы значения IERR дает возможность пользователю контpолировать точность вычисления интеграла.

Возможны случаи, когда аргумент функции за счет погрешности округления выйдет за границу области интегрирования. Поэтому, иногда, разумно доопределить функцию в окрестности граничной точки нулем, при этом достаточно взять радиус окрестности равным 10 - 9.

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

Unit TQS40R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FQS40R_p, QS40R_p;

function TQS40R: String;

implementation

function TQS40R: String;
var
IERR,_i :Integer;
RINT,EPS :Real;
А1 :Array [0..1] of Real;
A2 :Array [0..1] of Real;
begin
Result := '';  { результат функции }
A1[0] := 0.0;
A1[1] := 0.0;
A2[0] := 1.0;
A2[1] := 2.0;
EPS := 0.001;
QS40R(RINT,A1,A2,FQS40R,EPS,IERR);
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[A1[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
for _i:=0 to 1 do
 begin
  Result := Result + Format('%20.16f ',[A2[_i]]);
  if ( ((_i+1) mod 2)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f %10d ',[RINT,IERR]) + #$0D#$0A;
UtRes('TQS40R',Result);  { вывод результатов в файл TQS40R.res }
exit;
end;

end.

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

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

implementation

function FQS40R( X :Array of Real): Real;
begin
{ Result - прототип имени функции FQS40R на FORTRANe }
Result := 12*X[0]*X[1]/(1+X[0]*X[0])*(1+X[0]*X[0])
    *(1+X[1]*X[1])*(1+X[1]*X[1]);
exit;
end;

end.

Результаты:   RINT  =  124. ,   IERR  =  0