Библиотека численных методов на языке Fortran 90. Вычисление сумм, произведений и коэффициентов. Преобразование Эйлера

Автор cайта:
Владимир
Потемкин
fortran-90@yandex.ru

Вычисление сумм, произведений и коэффициентов

SM10. Вычисление суммы бесконечного ряда с помощью преобразования Эйлера.

CNM. Вычисление биномиальных коэффициентов.

SM10

Программа SM10 вычисляет сумму бесконечного ряда с помощью преобразования Эйлера, которое заключаетсмя в следующем. Если ряд в левой части равенства
_∞_∞
∑ a_i = ∑ (∆ⁱa₀)/2ⁱ⁺¹
ⁱ⁼⁰ⁱ⁼⁰
сходится, то это равенство справедливо и ряд в правой части равенства сходится быстрее. Здесь ∆ⁱa₀ обозначает разность вперед i-го порядка для a₀
∆a₀ = a₁ - a₀
∆¹a₀ = ∆a₁ - ∆a₀
∆²a₀ = ∆¹a₁ - ∆¹a₀
. . . . . . . . . .
Программа особенно эффективна в случае медленно сходящихся, а также расходящихся знакопеременных рядов. Суммирование прекращается, как только Tim раз подряд абсолютные значения членов преобразованного ряда будут меньше, чем Eps.

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

SM10(Fct, Eps, Tim)

Параметры программы

Fct, Eps, Tim - входные параметры;
Real SM10 - возвращаемое значение;

Real Fct(Integer i) - функция, вычисляющая значения членов ряда;
Real Eps - погрешность вычисления суммы;
Integer Tim - параметр повторения.

Пример

!Вычисление суммы знакопеременного ряда
program TestSM10
use NML
implicit none
real, parameter:: pi=3.1415926, Eps=1.E-8
integer, parameter:: Tim=4
real:: S, E
!begin
  S=SM10(Pi4, Eps, Tim)
  E=0.25*pi
  print 10, S, E
10 format(/'  Summa =',E18.10/'   pi/4 =',E18.10)
!end

contains

real function Pi4(i)
integer, intent(in):: i
  Pi4=1.0/float(2*i+1)
  if(mod(i,2) /= 0) Pi4=-Pi4
  return
end function Pi4

end program TestSM10

  Summa =  0.7853981256E+00
   pi/4 =  0.7853981256E+00

Вернуться к оглавлению Скачать SM10

DFT

Программа DFT вычисляет дискретное (прямое и обратное) быстрое преобразование Фурье с прореживанием по времени. Алгоритм Кули-Тьюки.
Прямое и обратное дискретное преобразование Фурье задаётся формулами
DFT-defenithon
где s(n) - исходный сигнал, S(k) - его Фурье-образ, j - мнимая единица. Предполагается, что сигнал s(n) - периодический, в этом случае Фурье-образ будет тоже периодическим с тем же периодом.

Основная идея быстрого преобразования Фурье (БПФ) состоит в разбиении отсчётов сигнала на две части и вычислении БПФ уже для каждой из частей меньшей размерности N/2 (прореживание по времени). Во втором варианте БПФ (прореживание по частоте) на две части разбивается спектр сигнала. Подробное описание обоих алгоритмов можно найти в многочисленной литературе.

[Б14]

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

call DFT(N, S, Error, Dir)

Параметры программы

N, Dir - входные параметры;
S - входной и выходной параметр;
Error - выходной параметр;

Integer N - количество отсчётов сигнала или Фурье-образа, N должно быть степенью числа 2, большей нуля;
Complex S[0:N-1] - массив с отсчётами;
Integer Dir - направление преобразования, Dir>0 или отсутствует - прямое преобразование Фурье, Dir<0 - обратное преобразование Фурье;
Integer Error - индикатор ошибки, Error=0 в случае нормального завершения программы, Error>0 если размерность массива S не является степенью числа 2.

Пример

! Прямое и обратное преобразование Фурье треугольного сигнала
program DigitalFourierTransform
    use NML
    implicit none
! Variables
    integer, parameter:: n=8, n1=n-1
    complex S(0:n1)
    real Arg(0:n1), X(0:n1)
    real a, b, h, t
    integer k, err
    data a/0.0/, b/2.0/
! Body of DigitalFourierTransform
    h=(b-a)/n;
    do k=0, n1
        Arg(k)=a+h*k
        if(k<=n/2) X(k)=h*k
        if(k>n/2) X(k)=X(k-1)-h
        S(k)=cmplx(X(k), 0.0)
    enddo
    print *, 'Original signal'
    print 10, (k, Arg(k), X(k), k=0, n1)
    call DFT(n, S, err)
    print *, 'Error=',err
    print *, 'The direct transform signal'
    print 11, (k, Arg(k), Real(S(k)), Aimag(S(k)), k=0, n1)
    call DFT(n, S, err, -1)
    print *, 'The inverse transform signal'
    print 11, (k, Arg(k), Real(S(k)), Aimag(S(k)), k=0, n1)
10  format(I2,'  Arg=',F4.2,F10.5)
11  format(I2,'  Arg=',F4.2,'  Re=',F10.6,'  Im=',F10.6)
end program DigitalFourierTransform

 Original signal
 0  Arg=0.00   0.00000
 1  Arg=0.25   0.25000
 2  Arg=0.50   0.50000
 3  Arg=0.75   0.75000
 4  Arg=1.00   1.00000
 5  Arg=1.25   0.75000
 6  Arg=1.50   0.50000
 7  Arg=1.75   0.25000
 Error=           0
 The direct transform signal
 0  Arg=0.00  Re=  4.000000  Im=  0.000000
 1  Arg=0.25  Re= -1.707107  Im=  0.000000
 2  Arg=0.50  Re=  0.000000  Im=  0.000000
 3  Arg=0.75  Re= -0.292893  Im= -0.000000
 4  Arg=1.00  Re=  0.000000  Im=  0.000000
 5  Arg=1.25  Re= -0.292893  Im=  0.000000
 6  Arg=1.50  Re=  0.000000  Im=  0.000000
 7  Arg=1.75  Re= -1.707107  Im=  0.000000
 The inverse transform signal
 0  Arg=0.00  Re=  0.000000  Im=  0.000000
 1  Arg=0.25  Re=  0.250000  Im= -0.000000
 2  Arg=0.50  Re=  0.500000  Im=  0.000000
 3  Arg=0.75  Re=  0.750000  Im= -0.000000
 4  Arg=1.00  Re=  1.000000  Im=  0.000000
 5  Arg=1.25  Re=  0.750000  Im=  0.000000
 6  Arg=1.50  Re=  0.500000  Im=  0.000000
 7  Arg=1.75  Re=  0.250000  Im=  0.000000

Вернуться к оглавлению Скачать DFT

CNM

Программа CNM вычисляет биномиальные коэффициенты C_n^m = n!/(m!(n-m)!) по рекуррентной формуле C_nⁱ⁺¹ = (n-i)/(i+1)C_nⁱ для i = 0, 1, ..., m-1, причем C_n⁰ = 1.

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

CNM(n, m)

Параметры программы

Integer n, m - входные параметры;
Integer CNM - возвращаемое значение.

Вернуться к оглавлению Скачать CNM