Текст подпрограммы и версий
qsk3r_p.zip  qsk3e_p.zip  		Тексты тестовых примеров
tqsk3r_p.zip  tqsk3e_p.zip 

Подпрограмма:  QSK3R (модуль QSK3R_p)

Назначение

Вычисление определенного интеграла от векторной функции по обобщенной квадратурной формуле Ньютона - Котесса шестого порядка точности.

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

Подпрограмма QSK3R вычисляет интеграл
      B
     ∫ F(x) dx       от векторной функции    F(x) = ( f1 (x), ... , fN (x) )
   A
с  покомпонентной погрешностью
                         B
      ЕРS ( 1 + |  ∫  fi (x) dx | )
                        A

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

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

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

procedure QSK3R(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 длины 4N;
IERR - целая переменная, служащая для диагностических сообщений:
IERR = 65 - когда заданная точность не может быть достигнута при максимально возможном числе (1048576) частичных отрезков разбиения. 
   Версии
QSK3E - вычисление с расширенной (Extended) точностью определенного интеграла от векторной функции по обобщенной квадратурной формуле Ньютона - Котесса шестого порядка точности. 
   Вызываемые подпрограммы
UTQS10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QSK3R.
UTQS12 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QSK3E. 
   Замечания по использованию
  1. 

Первый оператор подпрограммы F должен иметь вид:
procedure F (X :Real; var Y :Array of Real; I :Array of Integer);
Здесь:
X - аргумент функции F (тип: вещественный);
Y - вещественный вектоp длины N вычисленных значений функции;
I - целый массив размерности N; если его i - ая компонента отлична от нуля, то интеграл от этой компоненты функции F не вычисляется.

  2. 

При EPS ≤ 0 происходит только один просчет при заданном значении K.

  3. 

Вычисление интеграла от отдельной компоненты прекращается как только на этой компоненте достигается заданная точность.

  4. 

Если для каких - либо компонент векторной функции заданная погрешность интегрирования не может быть достигнута, то соответствующие компоненты вектоpа IND полагаются равными нулю.

  5. 

Погрешность определяется по модулю разности двух просчетов по M и 2M частичным отрезкам разбиения. Если заданная точность не достигнута, то значение M удваивается.

  6.  В подпрограмме QSK3E параметры RINT, A, B, F, EPS, R имеют тип Extended.

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

Unit TQSK3R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FQSK3R_p, QSK3R_p;

function TQSK3R: String;

implementation

function TQSK3R: String;
var
I,IERR :Integer;
RINT :Array [0..4] of Real;
R :Array [0..19] of Real;
IND :Array [0..4] of Integer;
begin
Result := '';  { результат функции }
QSK3R(RINT,0.0,1.0,FQSK3R,5,1.E-6,1,IND,R,IERR);
RЕSULt := Result + #$0D#$0A;
for I:=1 to 5 do
 begin
  Result := Result + Format('%20.16f      ',
 [RINT[I-1],IND[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%5d ',[IERR]) + #$0D#$0A;
UtRes('TQSK3R',Result);  { вывод результатов в файл TQSK3R.res }
exit;
end;

end.

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

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

implementation

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

end.

Результаты:

               RINT                       IND

      2.5000000000-01               1
      1.9999999999-01               1
      1.6666666666-01               1
      1.4285714427-01               4
      1.2500000496-01               4

      IERR  =  0