```Finally fixed it!

Here's a working Fortran subroutine for using the FFTW DCT to calculate
derivatives.  The error in the above two posted codes is an off-by-one on
the multiplier in the recurrence.  Also there is a lurking sign error that
I'm not quite able to trace from the FFTW definitions all the way through
to the recurrence.

subroutine dct_diff(dydx, y, n)
! use the discrete cosine transform from fftw to calculate the
derivative
integer, intent(in) :: n ! number of data points
double precision, intent(in), dimension(n) :: y ! function values
double precision, intent(out), dimension(n) :: dydx ! derivative
values
double precision, dimension(n) :: beta ! derivative coefficients
double precision, dimension(n) :: alpha ! function coefficients
integer :: plan1, plan2, i
integer :: n_logical ! the logical size of the transform, size of
! the corresponding real symmetric DFT

! the logical size depends on the type of transform, check the docs:
! http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
n_logical = 2*(n)

! forward DCT:
call dfftw_plan_r2r_1d(plan1, n, y, alpha, FFTW_REDFT10,
FFTW_ESTIMATE)
call dfftw_execute_r2r(plan1, y, alpha)
call dfftw_destroy_plan(plan1)

! recurrence for the derivative coefficients:
beta(n) = 0D0
beta(n-1) = 2D0 * dble(n-1) * alpha(n)
do i = n-1, 2, -1
beta(i-1) = beta(i+1) + 2D0 * dble(i-1) * alpha(i)
end do
beta = - beta ! this makes it work, but not sure why!

! inverse DCT:
call dfftw_plan_r2r_1d(plan2, n, beta, dydx, FFTW_REDFT01,
FFTW_ESTIMATE)
call dfftw_execute_r2r(plan2, beta, dydx)
call dfftw_destroy_plan(plan2)

! FFTW computes the un-normalized transforms, normalize by logical
size
dydx = dydx / dble(n_logical)

end subroutine dct_diff

```
```>
>>For one thing, you may be missing a scale factor.  FFTW computes the
>>unnormalized transform, so you need to multiply somewhere by a scale
>>factor yourself to get the inverse.  See:
>>
>>http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
>>
>>More generally, you need to pay close attention to FFTW's
>>normalization convention compared to the reference you are using for
>>the Chebyshev polynomials; compare the formulas in the manual above to
>>
>>Regards,
>>Steven G. Johnson
>>
>
>Hello Steven,
>
carried
>over with foactor of 2(n-1). I have taken care of that in the code right
>after doing the DCT transforms.
>
>      call dfftw_plan_r2r_1d ( plan, N, in, out,
>     ^                         FFTW_REDFT00, FFTW_ESTIMATE )
>      call dfftw_execute ( plan )
>      call dfftw_destroy_plan ( plan )
>      out = out/6
>
>I have already tested the inverse chebyshev transformation for the input
>function and it works fine. But if I perform the recurrence relation to
>calculate the differentiation coefficients and then use inverse
chebyshev,
>I am not able to obtain the derivative function.
>
>Regards,
>Sasidhar
>
>

Sasidhar,

Did you ever get your spectral derivative code to work?  I'm trying to do
the same thing, using the DCT from FFTW and following the recurrence in
Boyd's Spectral Methods book (pp. 298), but no luck for me either.

Here's my fortran (pretty much the same, but using different DCT
routines):
call dfftw_plan_r2r_1d(plan1, n, in, out, FFTW_REDFT10,
FFTW_ESTIMATE)
call dfftw_execute_r2r(plan1, in, out)
call dfftw_destroy_plan(plan1)
d(n) = 0.0
d(n-1) = 0.0
do i=n-2,2,-1
d(i) = d(i+2) + 2*(i+1)*out(i+1)
end do
d(1) = (d(3) + 2*2*out(2))/2.0
d = d/(2.0*dble(n))
call dfftw_plan_r2r_1d(plan2, n, d, out, FFTW_REDFT01, FFTW_ESTIMATE)
call dfftw_execute_r2r(plan2, d, out)
call dfftw_destroy_plan(plan2)

The test function I've been using is f(x) = x*exp(x) on [-1,1], and the
results from the above give me something that looks sort of like the
derivative flipped upside-down.

Anyone out there successfully accomplished calculating derivatives with
the discrete cosine transform that could lend a clue or two?

Cheers,
Josh

```
```On Nov 22, 6:31&#2013266080;pm, "sasidhar" <sasidha...@gmail.com> wrote:
> I have already tested the inverse chebyshev transformation for the input
> function and it works fine. But if I perform the recurrence relation to
> calculate the differentiation coefficients and then use inverse chebyshev,
> I am not able to obtain the derivative function.

I haven't looked at your recurrence relation, but you should realize
that there are several different normalization conventions for DCTs
besides just overall scale factors.  e.g. sometimes the x_0 element is
scaled differently to make the transform unitary, whereas FFTW's
normalization is chosen to enforce the correspondence with a DFT of
real-even data.  You should carefully compare FFTW's definition in the
manual with the definition in the reference you are using for the

Regards,
Steven G. Johnson
```
```>For one thing, you may be missing a scale factor.  FFTW computes the
>unnormalized transform, so you need to multiply somewhere by a scale
>factor yourself to get the inverse.  See:
>
>http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html
>
>More generally, you need to pay close attention to FFTW's
>normalization convention compared to the reference you are using for
>the Chebyshev polynomials; compare the formulas in the manual above to
>
>Regards,
>Steven G. Johnson
>

Hello Steven,

Thanks for the reply. I guess you are talking about normalization carried
over with foactor of 2(n-1). I have taken care of that in the code right
after doing the DCT transforms.

call dfftw_plan_r2r_1d ( plan, N, in, out,
^                         FFTW_REDFT00, FFTW_ESTIMATE )
call dfftw_execute ( plan )
call dfftw_destroy_plan ( plan )
out = out/6

I have already tested the inverse chebyshev transformation for the input
function and it works fine. But if I perform the recurrence relation to
calculate the differentiation coefficients and then use inverse chebyshev,
I am not able to obtain the derivative function.

Regards,
Sasidhar

```
```On Nov 21, 5:00&#2013266080;pm, "sasidhar" <sasidha...@gmail.com> wrote:
> Here dout array gives the manipulated coefficients and out array gives the
> chebyshev coefficients. However when I try to calculate the inverse, I am
> not getting the derivative function of the original.

For one thing, you may be missing a scale factor.  FFTW computes the
unnormalized transform, so you need to multiply somewhere by a scale
factor yourself to get the inverse.  See:

http://www.fftw.org/doc/1d-Real_002deven-DFTs-_0028DCTs_0029.html

More generally, you need to pay close attention to FFTW's
normalization convention compared to the reference you are using for
the Chebyshev polynomials; compare the formulas in the manual above to

Regards,
Steven G. Johnson
```
```>Basically it is a form of expansion of a function in terms of
>Chebyshev polys where they are sampled on a nonuniformly spaced set of
>nodes.  The nodes are found by cos(k*pi/N) where k goes from
>0,1,2,3,.., N. Using a set from -1 to +1. Scaling and translation may
>be used for other domains. The transform is also self inverting.
>On the "chebyshev nodes", the discrete chebyshev functions (formed by
>sampling the continuous ones) are orthogonal with respect to the
>chebyshev weighting function.
>
>Chebys have a bunch of neat properties such as a composition of chebys
>is also a chebys! But by doing a cosine transformation on the domain,
>the cheby transform can be found via a FFT. And since it is self
>inverting, the OP wants to take a sample function and go around the
>loop and see if all of his steps are working properly.
>
>He needs to do a discrete cheby transform by the direct method so as
>to compare that halfway step to what he gets from the FFT method. Then
>he can try the inverse to see if he gets back to where he started as a
>sanity check.

Hello Clay,

Once the chebyshev coefficients are calculated, they can be manipulated in
a way so that the inverse of the transformations gives a derivative of the
original function to which chebyshev transformation has been carried upon.
This is performed using the recurrence relation that I have shown above.

"dout(N) = 0
dout(N-1) = 2*(N)*out(N)

do j = N-2,1,-1
dout(j) = dout(j+2) + 2*(j+1)*out(j+1)
enddo"

Here dout array gives the manipulated coefficients and out array gives the
chebyshev coefficients. However when I try to calculate the inverse, I am
not getting the derivative function of the original.

Regards,
Sasidhar
```
```On Nov 21, 2:52&#2013266080;pm, Jerry Avins <j...@ieee.org> wrote:
> c...@claysturner.com wrote:
> > On Nov 21, 11:42 am, "sasidhar" <sasidha...@gmail.com> wrote:
> >> Hi all,
>
> >> I am trying to calculate derivatives of chebyshev coefficients using FFTW
> >> for cosine(x). However I am not getting back the right values when I
> >> perform inverse chebyshev transform. Can anyone check the code below and
> >> tell me the possible errors. Thank you.
>
> >> &#2013266080; &#2013266080; &#2013266080; program test
> >> &#2013266080; &#2013266080; &#2013266080; implicit none
> >> #include "fftw3.f"
>
> >> &#2013266080; &#2013266080; &#2013266080; INTEGER, PARAMETER &#2013266080;:: &#2013266080;N=4
> >> &#2013266080; &#2013266080; &#2013266080; INTEGER(8) :: &#2013266080;plan
> >> &#2013266080; &#2013266080; &#2013266080; REAL(8), DIMENSION(N) :: &#2013266080;in, out, out1, in1
> >> &#2013266080; &#2013266080; &#2013266080; REAL(8), DIMENSION(N) :: dout
> >> &#2013266080; &#2013266080; &#2013266080; REAL :: &#2013266080;pi, y
> >> &#2013266080; &#2013266080; &#2013266080; INTEGER :: j, p
>
> >> &#2013266080; &#2013266080; &#2013266080; pi = 3.14
> >> &#2013266080; &#2013266080; &#2013266080; write(*,*) 'Input array:'
> >> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; y =-1*cos(pi*(j)/N)
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; in(j) = cos(2*pi*y/N)
> >> &#2013266080; &#2013266080; &#2013266080; enddo
>
> >> &#2013266080; &#2013266080; &#2013266080; call dfftw_plan_r2r_1d ( plan, N, in, out,
> >> &#2013266080; &#2013266080; &#2013266080;^ &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; FFTW_REDFT00, FFTW_ESTIMATE )
> >> &#2013266080; &#2013266080; &#2013266080; call dfftw_execute ( plan )
> >> &#2013266080; &#2013266080; &#2013266080; call dfftw_destroy_plan ( plan )
> >> &#2013266080; &#2013266080; &#2013266080; out = out/6
> >> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; write(*,*) ' &#2013266080; &#2013266080;out(',j,') = ',out(j)
> >> &#2013266080; &#2013266080; &#2013266080; enddo
>
> >> &#2013266080; &#2013266080; &#2013266080; dout(N) = 0
> >> &#2013266080; &#2013266080; &#2013266080; dout(N-1) = 2*(N)*out(N)
>
> >> &#2013266080; &#2013266080; &#2013266080; do j = N-2,1,-1
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; dout(j) = dout(j+2) + 2*(j+1)*out(j+1)
> >> &#2013266080; &#2013266080; &#2013266080; enddo
>
> >> &#2013266080; &#2013266080; &#2013266080; call dfftw_plan_r2r_1d(plan, N+1, dout, out1,
> >> &#2013266080; &#2013266080; &#2013266080;^ &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; FFTW_REDFT00, FFTW_ESTIMATE)
> >> &#2013266080; &#2013266080; &#2013266080; call dfftw_execute ( plan )
> >> &#2013266080; &#2013266080; &#2013266080; call dfftw_destroy_plan ( plan )
>
> >> &#2013266080; &#2013266080; &#2013266080; write(*,*) 'Output array after inverse FFT:'
> >> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; write(*,*) ' &#2013266080; &#2013266080;out1(',j,') = ',out1(j)
> >> &#2013266080; &#2013266080; &#2013266080; enddo
>
> >> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; y =-1*cos(pi*(j)/N)
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; in1(j) = -1*sin(2*pi*y/N)
> >> &#2013266080; &#2013266080; &#2013266080; &#2013266080; write(*,*) ' &#2013266080; &#2013266080;in1(',j,') = ',in1(j)
> >> &#2013266080; &#2013266080; &#2013266080; enddo
> >> &#2013266080; &#2013266080; &#2013266080; end
>
> >> &#2013266080;Regards,
> >> &#2013266080;Sasidhar
>
> > Hello Sasidhar,
>
> > I don't know about your code, but derivatives of coefficients
> > (constants?) are zero. Or are you trying to find derivatives of
> > Chebyshev polynomials? These do have well defined forms in terms
> > Chebyshev polynomials (although both 1st and 2nd kind are needed
> > though). Or you can just stay with polynomials.
>
> > Maybe if you can give a clearer description of what you are really
> > trying to do, then we may be better able to help.
>
> > Are you trying to do differentiation via FFTs? You know aliasing can
> > be an issue here!
>
> Clay, What is a Chebychev transform?
>
> Jerry
> --
> Engineering is the art of making what you want from things you can get.
> &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;- Hide quoted text -
>
> - Show quoted text -

Hey Jerry,

Basically it is a form of expansion of a function in terms of
Chebyshev polys where they are sampled on a nonuniformly spaced set of
nodes.  The nodes are found by cos(k*pi/N) where k goes from
0,1,2,3,.., N. Using a set from -1 to +1. Scaling and translation may
be used for other domains. The transform is also self inverting.
On the "chebyshev nodes", the discrete chebyshev functions (formed by
sampling the continuous ones) are orthogonal with respect to the
chebyshev weighting function.

Chebys have a bunch of neat properties such as a composition of chebys
is also a chebys! But by doing a cosine transformation on the domain,
the cheby transform can be found via a FFT. And since it is self
inverting, the OP wants to take a sample function and go around the
loop and see if all of his steps are working properly.

He needs to do a discrete cheby transform by the direct method so as
to compare that halfway step to what he gets from the FFT method. Then
he can try the inverse to see if he gets back to where he started as a
sanity check.

IHTH,
Clay

Clay

```
```clay@claysturner.com wrote:
> On Nov 21, 11:42 am, "sasidhar" <sasidha...@gmail.com> wrote:
>> Hi all,
>>
>> I am trying to calculate derivatives of chebyshev coefficients using FFTW
>> for cosine(x). However I am not getting back the right values when I
>> perform inverse chebyshev transform. Can anyone check the code below and
>> tell me the possible errors. Thank you.
>>
>>       program test
>>       implicit none
>> #include "fftw3.f"
>>
>>       INTEGER, PARAMETER  ::  N=4
>>       INTEGER(8) ::  plan
>>       REAL(8), DIMENSION(N) ::  in, out, out1, in1
>>       REAL(8), DIMENSION(N) :: dout
>>       REAL ::  pi, y
>>       INTEGER :: j, p
>>
>>       pi = 3.14
>>       write(*,*) 'Input array:'
>>       do j = 1,N,1
>>         y =-1*cos(pi*(j)/N)
>>         in(j) = cos(2*pi*y/N)
>>       enddo
>>
>>       call dfftw_plan_r2r_1d ( plan, N, in, out,
>>      ^                         FFTW_REDFT00, FFTW_ESTIMATE )
>>       call dfftw_execute ( plan )
>>       call dfftw_destroy_plan ( plan )
>>       out = out/6
>>       do j = 1,N,1
>>         write(*,*) '    out(',j,') = ',out(j)
>>       enddo
>>
>>       dout(N) = 0
>>       dout(N-1) = 2*(N)*out(N)
>>
>>       do j = N-2,1,-1
>>         dout(j) = dout(j+2) + 2*(j+1)*out(j+1)
>>       enddo
>>
>>       call dfftw_plan_r2r_1d(plan, N+1, dout, out1,
>>      ^                         FFTW_REDFT00, FFTW_ESTIMATE)
>>       call dfftw_execute ( plan )
>>       call dfftw_destroy_plan ( plan )
>>
>>       write(*,*) 'Output array after inverse FFT:'
>>       do j = 1,N,1
>>         write(*,*) '    out1(',j,') = ',out1(j)
>>       enddo
>>
>>       do j = 1,N,1
>>         y =-1*cos(pi*(j)/N)
>>         in1(j) = -1*sin(2*pi*y/N)
>>         write(*,*) '    in1(',j,') = ',in1(j)
>>       enddo
>>       end
>>
>>  Regards,
>>  Sasidhar
>
> Hello Sasidhar,
>
> (constants?) are zero. Or are you trying to find derivatives of
> Chebyshev polynomials? These do have well defined forms in terms
> Chebyshev polynomials (although both 1st and 2nd kind are needed
> though). Or you can just stay with polynomials.
>
> Maybe if you can give a clearer description of what you are really
> trying to do, then we may be better able to help.
>
> Are you trying to do differentiation via FFTs? You know aliasing can
> be an issue here!

Clay, What is a Chebychev transform?

Jerry
--
Engineering is the art of making what you want from things you can get.
&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
```
```On Nov 21, 11:42&#2013266080;am, "sasidhar" <sasidha...@gmail.com> wrote:
> Hi all,
>
> I am trying to calculate derivatives of chebyshev coefficients using FFTW
> for cosine(x). However I am not getting back the right values when I
> perform inverse chebyshev transform. Can anyone check the code below and
> tell me the possible errors. Thank you.
>
> &#2013266080; &#2013266080; &#2013266080; program test
> &#2013266080; &#2013266080; &#2013266080; implicit none
> #include "fftw3.f"
>
> &#2013266080; &#2013266080; &#2013266080; INTEGER, PARAMETER &#2013266080;:: &#2013266080;N=4
> &#2013266080; &#2013266080; &#2013266080; INTEGER(8) :: &#2013266080;plan
> &#2013266080; &#2013266080; &#2013266080; REAL(8), DIMENSION(N) :: &#2013266080;in, out, out1, in1
> &#2013266080; &#2013266080; &#2013266080; REAL(8), DIMENSION(N) :: dout
> &#2013266080; &#2013266080; &#2013266080; REAL :: &#2013266080;pi, y
> &#2013266080; &#2013266080; &#2013266080; INTEGER :: j, p
>
> &#2013266080; &#2013266080; &#2013266080; pi = 3.14
> &#2013266080; &#2013266080; &#2013266080; write(*,*) 'Input array:'
> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; y =-1*cos(pi*(j)/N)
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; in(j) = cos(2*pi*y/N)
> &#2013266080; &#2013266080; &#2013266080; enddo
>
> &#2013266080; &#2013266080; &#2013266080; call dfftw_plan_r2r_1d ( plan, N, in, out,
> &#2013266080; &#2013266080; &#2013266080;^ &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; FFTW_REDFT00, FFTW_ESTIMATE )
> &#2013266080; &#2013266080; &#2013266080; call dfftw_execute ( plan )
> &#2013266080; &#2013266080; &#2013266080; call dfftw_destroy_plan ( plan )
> &#2013266080; &#2013266080; &#2013266080; out = out/6
> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; write(*,*) ' &#2013266080; &#2013266080;out(',j,') = ',out(j)
> &#2013266080; &#2013266080; &#2013266080; enddo
>
> &#2013266080; &#2013266080; &#2013266080; dout(N) = 0
> &#2013266080; &#2013266080; &#2013266080; dout(N-1) = 2*(N)*out(N)
>
> &#2013266080; &#2013266080; &#2013266080; do j = N-2,1,-1
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; dout(j) = dout(j+2) + 2*(j+1)*out(j+1)
> &#2013266080; &#2013266080; &#2013266080; enddo
>
> &#2013266080; &#2013266080; &#2013266080; call dfftw_plan_r2r_1d(plan, N+1, dout, out1,
> &#2013266080; &#2013266080; &#2013266080;^ &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; &#2013266080; FFTW_REDFT00, FFTW_ESTIMATE)
> &#2013266080; &#2013266080; &#2013266080; call dfftw_execute ( plan )
> &#2013266080; &#2013266080; &#2013266080; call dfftw_destroy_plan ( plan )
>
> &#2013266080; &#2013266080; &#2013266080; write(*,*) 'Output array after inverse FFT:'
> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; write(*,*) ' &#2013266080; &#2013266080;out1(',j,') = ',out1(j)
> &#2013266080; &#2013266080; &#2013266080; enddo
>
> &#2013266080; &#2013266080; &#2013266080; do j = 1,N,1
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; y =-1*cos(pi*(j)/N)
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; in1(j) = -1*sin(2*pi*y/N)
> &#2013266080; &#2013266080; &#2013266080; &#2013266080; write(*,*) ' &#2013266080; &#2013266080;in1(',j,') = ',in1(j)
> &#2013266080; &#2013266080; &#2013266080; enddo
> &#2013266080; &#2013266080; &#2013266080; end
>
> &#2013266080;Regards,
> &#2013266080;Sasidhar

Hello Sasidhar,

(constants?) are zero. Or are you trying to find derivatives of
Chebyshev polynomials? These do have well defined forms in terms
Chebyshev polynomials (although both 1st and 2nd kind are needed
though). Or you can just stay with polynomials.

Maybe if you can give a clearer description of what you are really
trying to do, then we may be better able to help.

Are you trying to do differentiation via FFTs? You know aliasing can
be an issue here!

Clay

```
```Hi all,

I am trying to calculate derivatives of chebyshev coefficients using FFTW
for cosine(x). However I am not getting back the right values when I
perform inverse chebyshev transform. Can anyone check the code below and
tell me the possible errors. Thank you.

program test
implicit none
#include "fftw3.f"

INTEGER, PARAMETER  ::  N=4
INTEGER(8) ::  plan
REAL(8), DIMENSION(N) ::  in, out, out1, in1
REAL(8), DIMENSION(N) :: dout
REAL ::  pi, y
INTEGER :: j, p

pi = 3.14
write(*,*) 'Input array:'
do j = 1,N,1
y =-1*cos(pi*(j)/N)
in(j) = cos(2*pi*y/N)
enddo

call dfftw_plan_r2r_1d ( plan, N, in, out,
^                         FFTW_REDFT00, FFTW_ESTIMATE )
call dfftw_execute ( plan )
call dfftw_destroy_plan ( plan )
out = out/6
do j = 1,N,1
write(*,*) '    out(',j,') = ',out(j)
enddo

dout(N) = 0
dout(N-1) = 2*(N)*out(N)

do j = N-2,1,-1
dout(j) = dout(j+2) + 2*(j+1)*out(j+1)
enddo

call dfftw_plan_r2r_1d(plan, N+1, dout, out1,
^                         FFTW_REDFT00, FFTW_ESTIMATE)
call dfftw_execute ( plan )
call dfftw_destroy_plan ( plan )

write(*,*) 'Output array after inverse FFT:'
do j = 1,N,1
write(*,*) '    out1(',j,') = ',out1(j)
enddo

do j = 1,N,1
y =-1*cos(pi*(j)/N)
in1(j) = -1*sin(2*pi*y/N)
write(*,*) '    in1(',j,') = ',in1(j)
enddo
end

Regards,
Sasidhar

```