БЧА НИВЦ МГУ. QST2R_P. Вычисление однократных интегралов (с контролем точности)

Текст подпрограммы и версий
qst2r_p.zip qst2e_p.zip

Тексты тестовых примеров
tqst2r_p.zip tqst2e_p.zip

Подпрограмма: QST2R (модуль QST2R_p)

Назначение

Вычисление определенного интеграла от векторной функции по обобщенной квадратурной формуле трапеций.

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

 Подпрограмма  QST2R вычисляет интеграл 
                     _B
                    ∫  f (x) dx 
                   ^A 
от векторной функции      F(x)  =  { f _i (x) } ,   i =1, ... , N
с покомпонентной  погрешностью 
                            _B
         ЕРS ( 1 + |  ∫  f _i (x) dx | )
                           ^A

по обобщенной квадратурной формуле трапеций с шагом h = (B - A)/M, где M - число частичных отрезков разбиения.

Я.М.Жилейкин, М.В.Соколовский. Набор стандартных программ для вычисления интегралов от векторных функций, Сб. "Методы и алгоритмы в численном анализе", Изд - во МГУ, 1981.

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

procedure QST2R(var RINT :Array of Real; A :Real; B :Real;
                F :Proc_F3_A; N :Integer; EPS :Real; K :Integer;
                var IND :Array of Integer; var R :Array of Real;
                var IERR :Integer);

   Параметры

RINT -	вещественный вектоp длины N, содержащий вычисленные значения интегралов;
A, B -	заданные нижний и верхний пределы интегрирования (тип: вещественный);
F -	имя подпрограммы, вычисляющей подинтегральные функции;
N -	размерность вектоpа - функции F (x) (тип: целый);
EPS -	заданная меpа погрешности вычисления интеграла (тип: вещественный);
K -	целая переменная, задающая начальное число частичных отрезков разбиения;
IND -	целый вектоp длины N, каждая компонента которого pавна числу частичных отрезков разбиения, при котоpом достигается заданная точность;
R -	вещественный рабочий вектоp длины 2 N;
IERR -	целая переменная, служащая для диагностических сообщений:

IERR=65 -

когда заданная точность не может быть достигнута при максимально возможном числе (1048576) частичных отрезков разбиения.

   Версии

QST2E -

вычисление с расширенной (Extended) точностью определенного интеграла от векторной функции по обобщенной квадратурной формуле трапеций.

   Вызываемые подпрограммы

UTQS10 -	подпрограмма выдачи диагностических сообщений при работе подпрограммы QST2R.
UTQS12 -	подпрограмма выдачи диагностических сообщений при работе подпрограммы QST2E.

   Замечания по использованию

	1.	Первый оператор подпрограммы F должен иметь вид: procedure F (X :Real; var Y :Array of Real; I :Array of Integer); Здесь: X - аргумент функции F (тип: вещественный); Y - вещественный вектоp длины N вычисленных значений функции F; I - целый вектоp длины N; если его i - ая компонента отлична от нуля, то интеграл от этой компоненты функции F не вычисляется.
	2.	При EPS ≤ 0 происходит только один просчет при заданном значении K.
	3.	Вычисление интеграла от отдельной компоненты прекращается как только на этой компоненте достигается заданная точность.
	4.	Если для каких - либо компонент векторной функции заданная погрешность интегрирования не может быть достигнута, то соответствующие компоненты вектоpа IND полагаются равными нулю.
	5.	Погрешность определяется по модулю разности двух просчетов по M и 2 M частичным отрезкам разбиения. Если заданная точность не достигнута, то значение M удваивается.
	6.	В подпрограмме QST2E параметры RINT, A, B, F, EPS, R имеют тип Extended.

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

Unit TQST2R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FQST2R_p, QST2R_p;

function TQST2R: String;

implementation

function TQST2R: String;
var
_i,IERR :Integer;
RINT :Array [0..4] of Real;
R :Array [0..9] of Real;
IND :Array [0..4] of Integer;
begin
Result := '';  { результат функции }
QST2R(RINT,0.0,1.0,FQST2R,5,0.5E-5,1,IND,R,IERR);
RЕSUlt := Result + #$0D#$0A;
for _i:=0 to 4 do
 begin
  Result := Result + Format('%20.16f ',[RINT[_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('%4d ',[IND[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TQST2R',Result);  { вывод результатов в файл TQST2R.res }
exit;
end;

end.

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

procedure FQST2R(X :Real; var Y :Array of Real;
                 I :Array of Integer);

implementation

procedure FQST2R(X :Real; var Y :Array of Real;
                 I :Array of Integer);
var
J :Integer;
label
_1;
begin
Y[0] := 1.0;
for J:=1 to 4 do
 begin
_1:
  Y[J+0] := Y[J-1]*X;
 end;
end;

end.

Результаты:

                RINT                          IND
        1.0000000000 + 00                 1
        5.0000000000 - 01                  1
        3.3333349227 - 01              512
        2.5000023841 - 01              512
        2.0000031788 - 01              512

     IERR  =  0