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 ?
DCT for computing the 2nd derivative of a function with Neumann Boundary Condition
Started by ●May 8, 2007