Текст подпрограммы и версий
qs21r_p.zip,  qs21e_p.zip  		Тексты тестовых примеров
tqs21r_p.zip,  tqs21e_p.zip 

Подпрограмма:  QS21R (модуль QS21R_p)

Назначение

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

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

Пусть  f (x) = f1 (x) + i f2 (x) - комплекснозначная функция вещественного аргумента, определенная на конечном отрезке [a, b]. Вычисляется серия из  N интегралов по формуле Симпсона по  N узлам: 

  b
 ∫ f (x) exp(-iwmx) dx  =
a
                                              1
                               =  (b-a)   ∫  f [ a+(b-a) y ] * exp { - iwm[ a+(b-a) y ]} dy  ≈
                                             0
                   N-1
 ≈ (b-a)/N    ∑   Ck f [ a+(b-a) k/N] * exp { - i wm[ a+(b-a) k/N ] }  =  ImN ,
                  k =0

 где  ImN  =  UmN + i VmN ,
        Ck - веса формулы Симпсона для отрезка [0,1] ,    i = √-1 ,
        wm = 2πm /(b-a) ,   m = 0, 1, 2, ..., N/2, N/2+1, N/2+2, ..., - 2, - 1 .

В.А.Морозов,Н.Н.Кирсанова,А.Ф.Сысоев. Комплекс алгоритмов быстрого преобразования Фурье дискретных рядов. Сб."Численный анализ на ФОРТРАНе", вып.15, Изд - во МГУ, 1976.

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

procedure QS21R(var U :Array of Real; var V :Array of Real; A :Real;
                B :Real; F :Proc_F3; N :Integer; var IERR :Integer);


   Параметры
U - вещественный вектор длины  N для вычисленных значений  UmN;
V - вещественный вектор длины  N для вычисленных значений  VmN;
A, B - заданные нижний и верхний пределы интегрирования (тип: вещественный);
F - имя подпрограммы, вычисляющей подинтегральную функцию  f (x);
N - число узлов интегрирования (тип: вещественный);
IERR - целая переменная для диагностических сообщений:
IERR=65 - когда заданное  N превосходит максимально допустимое число узлов. 
   Версии
QS21E - вычисление серии интегралов Фурье с помощью алгоритма быстрого преобразования Фурье (БПФ) от комплекснозначной функции вещественного аргумента по квадратурной формуле Симпсона в режиме расширенной (Extended) точности. 
   Вызываемые подпрограммы
UTQS10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QS21R.
UTQS11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QS21E.
FTF1C - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины L = 2M (M - целое) методом быстрого преобразования Фурье при работе подпрограммы QS21R.
FTF1Z - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины L = 2M (M - целое) методом быстрого преобразования Фурье при работе подпрограммы QS21E. 
   Замечания по использованию
  1. 

N должно быть степенью числа два.

  2. 

Максимальное значение  N не должно превосходить 220, т.е.  N ≤ 1048576 .

  3. 

Первый оператор подпрограммы  F должен иметь вид: 

       procedure F (X :Real; var F1 :Real; var F2 :Real);
       Здесь:
       X  - аргумент функции  f(x) (тип: вещественный);
       F1 - вещественная переменная,  F1 = Re f(x);
       F2 - вещественная переменная,  F2 = Im f(x).
  4. 

Для подпрограммы QS21E параметры U, V, A, B, а также параметры X, F1 и F2 из подпрограммы  F имеют тип Extended.

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

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

Unit TQS21R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FQS21R_p, QS21R_p;

function TQS21R: String;

implementation

function TQS21R: String;
var
IERR,K :Integer;
U :Array [0..511] of Real;
V :Array [0..511] of Real;
begin
Result := '';  { результат функции }
QS21R(U,V,0.0,1.0,FQS21R,512,IERR);
Result := Result + Format('%s',[' U(K)']) + #$0D#$0A;
RЕSUlt := Result + #$0D#$0A;
for K:=1 to 512 do
 begin
  Result := Result + Format('%16.7f ',[U[K-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' V(K)']) + #$0D#$0A;
Result := Result + #$0D#$0A;
for K:=1 to 512 do
 begin
  Result := Result + Format('%16.7f ',[V[K-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TQS21R',Result);  { вывод результатов в файл TQS21R.res }
exit;
end;

end.

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

procedure FQS21R(X :Real; var F1 :Real; var F2 :Real);

implementation

procedure FQS21R(X :Real; var F1 :Real; var F2 :Real);
begin
F1 := Exp(X);
F2 := 0.0;
end;

end.

Результаты:   

   Приводятся значения интегралов  UmN  и  VmN
   для  w0 = 0  и  w1 = 2π  при  N = 512:

      U(1) = 1.718281828         V(1) = 0.000000000 
      U(2) = 0.042449333         V(2) = 0.266717025