Текст подпрограммы и версий
qs21r_c.zip  qs21d_c.zip 
Тексты тестовых примеров
tqs21r_c.zip  tqs21d_c.zip 

Подпрограмма:  qs21r_c

Назначение

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

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

Пусть  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.

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

    int qs21r_c (real *u, real *v, real *a, real *b, S_fp f,
            integer *n, integer *ierr)

Параметры

u - вещественный вектор длины  n для вычисленных значений  Umn;
v - вещественный вектор длины  n для вычисленных значений  Vmn;
a, b - заданные нижний и верхний пределы интегрирования (тип: вещественный);
f - имя подпрограммы, вычисляющей подинтегральную функцию  f (x);
n - число узлов интегрирования (тип: вещественный);
ierr - целая переменная для диагностических сообщений:
ierr=65 - когда заданное  n превосходит максимально допустимое число узлов.

Версии

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

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

utqs10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qs21r_c;
utqs11_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы qs21d_c;
ftf1c_c - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины l = 2m (m - целое) методом быстрого преобразования Фурье при работе подпрограммы qs21r_c;
ftf1p_c - подпрограмма вычисления дискретного или обратного дискретного преобразования Фурье комплексного ряда длины l = 2m (m - целое) методом быстрого преобразования Фурье при работе подпрограммы qs21d_c.

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

  1. 

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

  2. 

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

  3. 

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

           int f(float *x, float *f1, float *f2)
       Здесь:
       x  - аргумент функции  f(x) (тип: вещественный);
       f1 - вещественная переменная,  f1 = re f(x);
       f2 - вещественная переменная,  f2 = im f(x) 
  4. 

Для подпрограммы qs21d_c параметры u, v, a, b, а также параметры x, f1 и f2 из подпрограммы  f имеют тип double;

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

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

int main(void)
{
    /* Local variables */
    static int ierr;
    extern int qs21r_c(float *, float *, float *, float *, U_fp,
                       int *, int *);
    extern int f_c();
    static int k;
    static float u[512], v[512];

    qs21r_c(u, v, &c_b1, &c_b2, (U_fp)f_c, &c__512, &ierr);

    printf("\n u[k] \n");
    for (k = 1; k <= 508; k += 3) {
         printf("\n %16.7f %16.7f %16.7f ",u[k - 1], u[k], u[k + 1]);
    }
    printf("\n %16.7f %16.7f \n ",u[510], u[511]);
    printf("\n v[k] \n");
    for (k = 1; k <= 508; k += 3) {
         printf("\n %16.7f %16.7f %16.7f ",v[k - 1], v[k], v[k + 1]);
    }
    printf("\n %16.7f %16.7f \n ",v[510], v[511]);
    return 0;
} /* main */

int f_c(float *x, float *f1, float *f2)
{
    /* Builtin functions */
    double exp(double);

    *f1 = (float)exp(*x);
    *f2 = 0.f;
    return 0;
} /* f_c */


Результаты:   

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

      u(1) = 1.718281828         v(1) = 0.000000000 
      u(2) = 0.042449333         v(2) = 0.266717025