Текст подпрограммы и версий
qsf9r_p.zip  qsf9e_p.zip  		Тексты тестовых примеров
tqsf9r_p.zip  tqsf9e_p.zip 

Подпрограмма:  QSF9R (модуль QSF9R_p)

Назначение

Вычисление определенных интегралов 

                  B
(1)             ∫  F (x) sin( ρ x + φ ) dx ,
                A
                  B
(2)             ∫  F (x) cos( ρ x + φ ) dx
                A

от векторной функции по формулам интерполяционного типа пятой степени точности.

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

Подпрограмма QSF9R вычисляет интегралы (1) и (2), или один из них от векторной функции F (x) = ( f1 (x), ..., fN (x) ) с покомпонентной погрешностью 

                     B
      E ( 1 + |  ∫  fi (x) sin( ρ x + φ) dx | )  или
                   A     
                     B                               
      E ( 1 + |  ∫  fi (x) cos( ρ x + φ) dx | )      
                   A

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

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

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

procedure QSF9R(var RINT :Array of Real; A :Real; B :Real;
                FUN :Proc_F3_A; RO :Real; FI :Real; N :Integer; EPS :Real;
                K :Integer; L1 :Integer; L2 :Integer;
                var IND :Array of Integer; var R :Array of Real;
                var IERR :Integer);


   Параметры
RINT - вещественный вектоp, содержащий вычисленные значения интегралов; его длина pавна 2N, если вычисляются оба интеграла (1) и (2), или N, если вычисляется один из этих интегралов;
A, B - заданные нижний и верхний пределы интегрирования (тип: вещественный);
F - имя подпрограммы, вычисляющей подинтегральные функции;
RO, FI - заданные значения параметров  ρ и  φ (тип: вещественный);
N - размерность вектора - функции F (x) (тип: целый);
E - заданная мера погрешности вычисления интеграла (тип: вещественный);
K - целая переменная, задающая начальное число частичных отрезков разбиения;
L1, L2 - задают режим работы подпрограммы (тип: целый); при этом, если:
 

L1 = 1 и L2 = 0, то вычисляется интеграл (1);

L1 = 0 и L2 = 1, то вычисляется интеграл (2);

L1 = 1 и L2 = 1, то вычисляются интегралы (1) и (2);
IND - целый вектоp длины N, каждая компонента котоpого pавна числу частичных отрезков разбиения, при котоpом достигается заданная точность по этой компоненте;
R - вещественный рабочий вектоp длины 3N, если вычисляются оба интеграла, и 2N, если вычисляется один из них;
IERR - целая переменная, служащая для диагностических сообщений:
IERR = 65 - если заданная точность не может быть достигнута при максимально возможном числе (1048576) частичных отрезков разбиения. 
   Версии
QSF9E - вычисление с расширенной (Extended) точностью определенных интегралов (1) и (2) от векторной функции по формулам интерполяционного типа пятой степени точности. 
   Вызываемые подпрограммы
UTQS10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QSF9R.
UTQS12 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QSF9E. 
   Замечания по использованию
  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 и 2M частичным отрезкам разбиения. Если заданная точность не достигнута, то значение M удваивается.

  6.  В подпрограмме QSF9E параметры RINT, A, B, F, RO, FI, E, R имеют тип Extended.

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

Unit TQSF9R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FQSF9R_p, QSF9R_p;

function TQSF9R: String;

implementation

function TQSF9R: String;
var
I,IERR :Integer;
PI0 :Real;
RINT :Array [0..9] of Real;
R :Array [0..14] of Real;
IND :Array [0..4] of Integer;
label
_1;
begin
Result := '';  { результат функции }
PI0 := 1.5707963268;

for I:=0 to 9 do
 begin
  RINT[I] := 0.e0;
 end;

QSF9R(RINT,-PI0,PI0,FQSF9R,1.0,0.0,5,1.E-4,2,1,1,IND,R,IERR);
for I:=1 to 5 do
 begin
  Result := Result + Format('%2d %20.16f %20.16f %5d ',
 [I,RINT[(I-1)+0],RINT[(I-1)+5],IND[I-1]]) + #$0D#$0A;
_1:
 end;
UtRes('TQSF9R',Result);  { вывод результатов в файл TQSF9R.res }
exit;
end;

end.

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

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

implementation

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

end.

Результаты:

              RINT (I, 1)                           RINT (I, 2)                 IND (I)

      1.9999999618 E+00             -8.661 E-12                            2
     -2.676 E-11                             9.3480217805 E-01              2
      2.8044065499 E+00             -6.428 E-12                            2
     -5.268 E-11                             9.5850994403 E-01              2
      4.7925452374 E+00             -2.098 E-11                            4

      IERR  =  0

      I  =   (1, 2, 3, 4, 5)