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

Question

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[0] = 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.

Answer 1

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[0] = 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].

Answer 2

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.