# Chebyshev differentiation via FFT with a domain [a,b]

I want to ask something about Chebyshev differentiation via FFT, which can be used to obtain with spectral accuracy the derivative of a smooth function. See for instance this code in python, which performs the derivative of a given function "u" in the domain [-1,1] and computes the L2 norm of the error.

from numpy.fft import fft,ifft
import numpy as np

def chebfft(v,x):
N = len(v)-1;
if N==0: return 0
ii = np.arange(0,N); iir = np.arange(1-N,0); iii = np.array(ii,dtype=int)
V = np.hstack((v,v[N-1:0:-1]))
U = np.real(fft(V))
W = np.real(ifft(1j*np.hstack((ii,[0.],iir))*U))
w = np.zeros(N+1)
w[1:N] = -W[1:N]/np.sqrt(1-x[1:N]**2)
w = sum(iii**2*U[iii])/N + .5*N*U[N]
w[N] = sum((-1)**(iii+1)*ii**2*U[iii])/N + .5*(-1)**(N+1)*N*U[N]
return w

N=50
L=10
aa = 0.34

x = np.cos(np.pi*np.arange(0,N+1)/N)
u=np.sin(2*np.pi*x/L)*np.exp(-aa*x)
uder = np.exp(-aa*x)*(2*np.pi*np.cos(2*np.pi*x/L)-aa*L*np.sin(2*np.pi*x/L))/L
err_ft = chebfft(u,x)-uder

print ("L2 norm:",np.linalg.norm(err_ft))


The problem that I have is that I don't really understand how to do the same with a general interval [a,b]. The new collocation points can be changed just with a linear mapping, but then the procedure inside chebfft fails. Provably the reason why is not woorking is very simple, but I don't see.

Thanks for the help.

Given a function $$u = u(X), \qquad X \in [X_{min},X_{max}]$$ Map to $$x \in [-1,1]$$ $$x = \frac{X - X_c}{\frac{1}{2}(X_{max} - X_{min})}, \qquad X_c = \frac{1}{2}(X_{min} + X_{max})$$ If $$D$$ is Chebyshev approximation of derivative $$D \approx \frac{d u}{d x}$$ then $$\frac{d u}{d X} = \frac{d x}{d X} \frac{d u}{d x} \approx \frac{2}{X_{max} - X_{min}} D$$ I modified your code for the function $$u(X) = X \sin(X), \qquad X \in [0,2\pi]$$ and it works fine.

from numpy.fft import fft,ifft
import numpy as np
import matplotlib.pyplot as plt

def chebfft(v,x):
N = len(v)-1;
if N==0: return 0
ii = np.arange(0,N); iir = np.arange(1-N,0); iii = np.array(ii,dtype=int)
V = np.hstack((v,v[N-1:0:-1]))
U = np.real(fft(V))
W = np.real(ifft(1j*np.hstack((ii,[0.],iir))*U))
w = np.zeros(N+1)
w[1:N] = -W[1:N]/np.sqrt(1-x[1:N]**2)
w = sum(iii**2*U[iii])/N + .5*N*U[N]
w[N] = sum((-1)**(iii+1)*ii**2*U[iii])/N + .5*(-1)**(N+1)*N*U[N]
return w

N = 50
Xmin, Xmax = 0.0, 2.0*np.pi
Xc = 0.5*(Xmin + Xmax)

x = np.cos(np.pi*np.arange(0,N+1)/N)
X = 0.5*(Xmax - Xmin)*x + Xc
u = X * np.sin(X)
uder = np.sin(X) + X * np.cos(X)
udercheb = chebfft(u,x)*2.0/(Xmax - Xmin)
err_ft = udercheb - uder
print ("L2 norm:",np.linalg.norm(err_ft))
plt.plot(X,uder,X,udercheb,'o')
plt.legend(('Exact','Chebyshev'))
plt.show()


So you just need a simple scaling to account for the domain being different from [-1,1].

• now I understand, thanks a lot!
– Joe
Nov 9, 2019 at 15:59

You need more than just simple linear mapping here. Let's say in general I have a complex variable $$z$$ as:

$$z = R e^{i\theta} + z_{C}$$

Your $$x$$ is defined as:

$$x = \Re\{z\} = R \cos(\theta) + x_{C}$$

The maximum and minimum values of $$x$$ could be calculated as:

$$\min(x) = x_{C} - R$$

$$\max(x) = x_{C} + R$$

Now you want to map into $$[a,b]$$ range, so:

$$x_{C} = \frac{b+a}{2}$$

$$R = \frac{b-a}{2}$$

If I take a derivative from $$x$$ with respect to $$\theta$$:

$$\frac{d x}{d \theta} = -R \sin(\theta) = -\sqrt{R^{2} - (x-x_{C})^{2}}$$

Finally Chebyshev polynomials are defined as:

$$T_{n}(x) = \Re \Bigg \{ \Big (\frac{z - z_{c}}{R} \Big )^{n} \Bigg \} = \cos(n\theta)$$

Furthermore, the derivative of Chebyshev polynomial is:

$$T^{'}_{n}(x) = -n \sin(n\theta) \frac{d \theta}{d x} = \frac{n \sin(n\theta)}{R \sin(\theta)}$$

Now, let's define discretized $$x_{j}$$ as:

$$x_{j} = R \cos(\theta_{j}) + x_{C}$$

You want to approximate your function $$f$$ by using Chebyshev polymonials as:

$$f(\theta) = \sum_{n=0}^{N} a_{n} \cos(n \theta)$$

The derivative of $$f$$, which you are mainly looking for is:

$$\frac{d f}{d x} = \sum_{n=0}^{N} \frac{n a_{n} \sin(n\theta)}{\sqrt{R^{2} - (x-x_{C})^{2}}}$$

So just find the $$a_{n}$$ from FFT and then at $$x = a$$ and $$x = b$$, from L'Hopital you have:

$$\frac{d f}{d x}|_{x = a} = \sum_{n = 0}^{N} (-1)^{n} \frac{n^{2}}{R} a_{n}$$

$$\frac{d f}{d x}|_{x = b} = \sum_{n = 0}^{N} \frac{n^{2}}{R} a_{n}$$

So, you need a little bit more effort to make it work for general $$[a,b]$$ range. I put the implementation for your homework.