Hello,
I am having a lot of trouble computing the Laplacian of a function
(with neumann boundary conditions) using fftw's DCT function. I do
this by taking the cosine transform using fftw, multiplying it by -k^2
(negative wave vectors squared), back transforming, and normalizing.
Here is my fortran code. Sorry that its a bit messy. It seems to work
when the periodic extention of the function lacks jump
discontinuities, but when these discontinuities exist the accuracy is
pretty bad. I expect this at the boundaries where the discontinuities
exist, but the accuracy is pretty bad for the whole range.
Program test_fftw
IMPLICIT NONE
INCLUDE '/usr/local/include/fftw3.f'
REAL (8), DIMENSION(0:399):: X, ft_X, X2
REAL (8), DIMENSION(0:399):: k2
REAL (8):: pi, L_x
INTEGER (8) :: planf, planb
INTEGER :: xgrid,i,k,j
DP = 8
xgrid = 400
L_x = 4
pi = 4*atan(1.0)
CALL dfftw_plan_r2r_1d( planf, xgrid, X, ft_X, FFTW_REDFT00,
FFTW_MEASURE )
CALL dfftw_plan_r2r_1d( planb, xgrid, ft_X, X2, FFTW_REDFT00,
FFTW_MEASURE )
! Compute k^2
DO i = 0, xgrid-1
IF ( ( i < xgrid/2 ) ) THEN
k2(i) = ( pi * REAL( i,8 ) / L_x )**2
ELSE IF (( i >= xgrid/2 )) THEN
k2(i) = ( pi * REAL(-xgrid + i ,8 ) / L_x )**2
END IF
END DO
!Define X
DO i = 0,xgrid-1
X(i) = sin( pi*REAL(i,8)*REAL(L_x,8)/(REAL(xgrid-1,8)))
END DO
call dfftw_execute(planf)
ft_X = -ft_X * k2
call dfftw_execute(planb)
X2 = X2 / REAL(2*(xgrid-1),8)
CALL dfftw_destroy_plan(planf)
CALL dfftw_destroy_plan(planb)
END Program