The following code calculates derivative of any sine or cosine trigonometric function using FFTW3 function.
program oned
implicit none
include "fftw3.f"
integer, parameter :: nx = 128
integer :: i
double complex, dimension(0:nx-1) :: input, solution
double precision, dimension(0:nx-1) :: x
double precision, parameter :: pi = acos(-1.0d0), dx = 2.0d0*pi/(nx)
double complex :: wrap_variable, num = (0, 1.0d0)
integer*8 :: plan_x_fwd, plan_x_bwd
call dfftw_plan_dft_1d(plan_x_fwd, nx, input, input, fftw_forward, fftw_measure)
call dfftw_plan_dft_1d(plan_x_bwd, nx, input, input, fftw_backward, fftw_measure)
do i = 0, nx-1
x(i) = i*dx
! function whose derivative is to be calculated.
input(i) = sin(8*x(i))
end do
! call fftw subroutine to calculate fourier coefficients of array
call dfftw_execute_dft(plan_x_fwd, input, input)
! normalization of array
input = input/nx
do i = 0, nx-1
if ( i .ge. 0 .and. i .le. nx/2-1 ) then
! calculating first derivatve. This involves multiplying by i*k
! where i**2 = -1 and k is the wavenumber
input(i) = input(i)*num*i
else if ( i .ge. nx/2 .and. i .le. nx-1) then
input(i) = input(i)*num*(i-nx)
end if
end do
do i = 0, nx-1
if ( i .gt. 0 .and. i .le. nx/2-1) then
! calculating second derivative.
input(i) = input(i)*num*i
else if ( i .ge. nx/2 .and. i .le. nx-1) then
input(i) = input(i)*num*(i-nx)
end if
end do
do i = 0, nx-1
if ( i .gt. 0 .and. i .le. nx/2-1) then
! calculating third derivative.
input(i) = input(i)*num*i .
else if ( i .ge. nx/2 .and. i .le. nx-1) then
input(i) = input(i)*num*(i-nx)
end if
end do
input(nx/2) = 0.0d0
call dfftw_execute_dft(plan_x_bwd, input, input)
! print the error in calculation.
print *, maxval(abs(real(solution - input)))
call dfftw_destroy_plan(plan_x_fwd)
call dfftw_destroy_plan(plan_x_bwd)
end program
The function I am using as "input" is sin(8x). When I calculate the derivative as above and compare against the answer which is 8*cos(8x), I get a maximum error of the order of E-14. When I calculate the second derivative and compare against the answer which is -64*sin(8x), then I get an error of the order E-12. Similarly for third derivative, I get E-10.
I find this very unsettling. I don't see why multiplying the number by i*k where i**2 = -1 and k is wavenumber causes precision loss.
Is this an issue with my code ? Or is it something to do with the nature of floating point numbers and if so why ?
