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





_____________________________________
Do you know a company who employs DSP engineers?  
Is it already listed at http://dsprelated.com/employers.php ?