Текст подпрограммы и версий
qs20r_p.zip,  qs20e_p.zip  		Тексты тестовых примеров
tqs20r_p.zip,  tqs20e_p.zip 

Подпрограмма:  QS20R (модуль QS20R_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    ∑   f [ a+(b-a) k/N] * exp { - i wm[ a+(b-a) k/N ] }  =  ImN ,
                  k =0

 где  ImN  =  UmN + i VmN ,     i = √-1 ,
        wm = 2πm /(b-a) ,   m = 0, 1, 2, ..., N/2, N/2+1, N/2+2, ..., - 2, - 1

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

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

procedure QS20R(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 превосходит максимально допустимое число узлов. 
   Версии
QS20E - вычисление серии интегралов Фурье с помощью алгоритма быстрого преобразования Фурье (БПФ) от комплекснозначной функции вещественного аргумента по квадратурной формуле прямоугольников в режиме расширенной (Extended) точности. 
   Вызываемые подпрограммы
UTQS10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QS20R.
UTQS11 - подпрограмма выдачи диагностических сообщений при работе подпрограммы QS20E.
FTF1C - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины L = 2M (M - целое) методом быстрого преобразования Фурье при работе подпрограммы QS20R.
FTF1Z - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины L = 2M (M - целое) методом быстрого преобразования Фурье при работе подпрограммы QS20E. 
   Замечания по использованию
  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. 

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

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

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

Unit TQS20R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FQS20R_p, QS20R_p;

function TQS20R: String;

implementation

function TQS20R: String;
var
IERR,K :Integer;
U :Array [0..511] of Real;
V :Array [0..511] of Real;
begin
Result := '';  { результат функции }
QS20R(U,V,0.0,1.0,FQS20R,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('TQS20R',Result);  { вывод результатов в файл TQS20R.res }
exit;
end;

end.

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

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

implementation

procedure FQS20R(X :Real; var F1 :Real; var F2 :Real);
begin
F1 := (X*X-X+1.0/6.0)*2.0;
F2 := 0.0;
end;

end.

Результаты:   

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

      U(1) = 0.000001272         V(1) = 0.000000000 
      U(2) = 0.101322455         V(2) = 0.000000000