Reply by tchantaw May 8, 20072007-05-08
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