2
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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].

$\endgroup$
  • $\begingroup$ now I understand, thanks a lot! $\endgroup$ – Joe Nov 9 at 15:59
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.