---
title: 9-1. Chebyshev approximations - derivative
tags: Spectral method
---

**Reference**
1. [Spectral method in Matlab](https://people.maths.ox.ac.uk/trefethen/spectral.html) by *Lloyd N. Trefethen*
2. Spectral methods by *J. Shen*, *T. Tang* and *L.-L. Wang*

---

% CHEB  compute D = differentiation matrix, x = Chebyshev grid

  function [D,x] = cheb(N)
  if N==0, D=0; x=1; return, end
  x = cos(pi*(0:N)/N)'; 
  c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
  X = repmat(x,1,N+1);
  dX = X-X';                  
  D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
  D  = D - diag(sum(D'));                 % diagonal entries

---

% CHEBFFT  Chebyshev differentiation via FFT. Simple, not optimal.  
%          If v is complex, delete "real" commands.

  function w = chebfft(v)
  N = length(v)-1; if N==0, w=0; return, end
  x = cos((0:N)'*pi/N);
  ii = 0:N-1;
  v = v(:); V = [v; flipud(v(2:N))];      % transform x -> theta          
  U = real(fft(V));
  W = real(ifft(1i*[ii 0 1-N:-1]'.*U));
  w = zeros(N+1,1);
  w(2:N) = -W(2:N)./sqrt(1-x(2:N).^2);    % transform theta -> x     
  w(1) = sum(ii'.^2.*U(ii+1))/N + .5*N*U(N+1);     
  w(N+1) = sum((-1).^(ii+1)'.*ii'.^2.*U(ii+1))/N + ...
              .5*(-1)^(N+1)*N*U(N+1);

---

% p11.m - Chebyshev differentation of a smooth function

  xx = -1:.01:1; uu = exp(xx).*sin(5*xx); clf
  for N = [10 20]
    [D,x] = cheb(N); u = exp(x).*sin(5*x);
      subplot('position',[.15 .66-.4*(N==20) .31 .28])
      plot(x,u,'.','markersize',14), grid on
      line(xx,uu)
      title(['u(x),  N=' int2str(N)])
    error = D*u - exp(x).*(sin(5*x)+5*cos(5*x));
      subplot('position',[.55 .66-.4*(N==20) .31 .28])
      plot(x,error,'.','markersize',14), grid on
      line(x,error)
      title(['  error in u''(x),  N=' int2str(N)])
  end

---

### C++ Code(亞衛)
#### Differentiation Matrix

#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<cmath>
#include<vector>
#include<conio.h>
#include<iomanip>

using namespace std;
#define N 10
#define pi 3.14159265358979323846

int main(void) {
	cout.precision(6); 
	vector<double> x(N+1,0);
	for (int i=0;i<=N;i++){
		x[i] = cos((pi*i) / N);
	}
	vector<double> c(N+1,0);
	c[0] = 2;
	c[N] = 2;
	for (int i=1;i<N;i++){
		c[i] = pow((-1.0),i);
	}
	double X[N+1][N+1];
	double Y[N+1][N+1];
	for(int i=0;i<=N;i++){
		for(int j=0;j<=N;j++){
			X[j][i] = x[j];
			Y[i][j] = x[j];    //Y = X transpose
		}
	}
	
	for(int i=0;i<=N;i++){
		for(int j=0;j<=N;j++){
			X[i][j] = X[i][j] - Y[i][j];   //dX = X-X'
		}
	}
	
	double C[N+1][N+1];
	double D[N+1][N+1];	
	for(int i=0;i<=N;i++){
		for(int j=0;j<=N;j++){
			C[i][j] = c[i]*(1.0/c[j]);
		}
	}
	for(int i=0;i<=N;i++)  X[i][i] = X[i][i] + 1.0;
	for(int i=0;i<=N;i++){
		for(int j=0;j<=N;j++){
			D[i][j] = C[i][j]/X[i][j];
		}
	}
	for(int i=0;i<=N;i++){
		for(int j=0;j<=N;j++){
			Y[i][j] = D[j][i];   //Y = D transpose
		}
	}
	vector<double> sum(N+1,0);
	for(int i=0;i<=N;i++){
		for(int j=0;j<=N;j++){
			sum[i] = sum[i] + Y[j][i];   //Y = D transpose
		}
	}
	for(int i=0;i<=N;i++){
		D[i][i] = D[i][i] - sum[i];
	}
	
	//Output Differentiation Matrix
	cout.setf(ios::fixed);
	for(int i=0;i<=N;i++){
		cout<<"[ ";
		for(int j=0;j<=N;j++){
			cout<<fixed<< setprecision(6) <<D[i][j]<<", ";
		}
		cout<<" ]"<<endl;
	}
}

Result:(with N=10)

#### Chebyshev differentiation via FFT

#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<cmath>
#include<vector>
#include<conio.h>
#include<complex>
#include<valarray>

using namespace std;

#define N 16
#define pi 3.14159265358979323846

typedef complex<double> Complex;
typedef valarray<Complex> CArray;

void Peek(vector<double> vector) {					//Output the vector
	for (int i = 0; i < vector.size(); i++) {
		cout <<"vec["<< i << "] = " << vector[i] << endl;
	}
}

void Peekcomplex(CArray& x){
		for(int i=0;i<=2*N-1;i++){
		cout <<"arr["<< i << "] = " << x[i] <<endl;
	}
}
void fft(CArray& x)
{
    const size_t K = x.size();
    if (K <= 1) return;
 
    // divide
    CArray even = x[std::slice(0, K/2, 2)];
    CArray  odd = x[std::slice(1, K/2, 2)];
 
    // conquer
    fft(even);
    fft(odd);
 
    // combine
    for (size_t k = 0; k < K/2; ++k)
    {
        Complex t = std::polar(1.0, -2 * pi * k / K) * odd[k];
        x[k    ] = even[k] + t;
        x[k+K/2] = even[k] - t;
    }
}
void ifft(CArray& x)
{
    // conjugate the complex numbers
    x = x.apply(std::conj);
 
    // forward fft
    fft( x );
 
    // conjugate the complex numbers again
    x = x.apply(std::conj);
 
    // scale the numbers
    x /= x.size();
}
int main(void) {
	vector<double> x(N+1,0);
	vector<double> v(N+1,0);
	for(int i=0;i<=N;i++) {  //create vector v
		x[i] = (-1) + (2.0/N)*i;
		v[i] = exp(x[i]) * sin(5*x[i]);
	}
	for(int i=0;i<=N;i++) x[i] = cos((pi*i) / N);
	Complex V[2*N] = {0.0};
	double U[2*N] = {0.0};
	for(int i=0;i<(2*N);i++) {  //create vector v
		if(i<=N)	V[i] = v[i];
		else {
			V[i] = v[2*N-i];
		}
	}
	CArray data(V, 2*N);
	fft(data);
	complex<double> temp;
	for(int i=0;i<(2*N);i++){
		temp = data[i];
		U[i] = temp.real();
	}
	vector<double> ii(2*N,0.0);
	for(int i=0;i<(2*N);i++) {
		if(i<N)		ii[i] = i;
		else{
			ii[i] = -ii[(2*N)-i];
		}
	}
	for(int i=0;i<(2*N);i++){
		temp.real() = 0.0;
		temp.imag() = ii[i] * U[i];
		data[i] = temp;
	}
	ifft(data);
	vector<double> W(2*N,0.0);
	for(int i=0;i<(2*N);i++){
		temp = data[i];
		W[i] = temp.real();
	}
	vector<double> w(N+1,0.0);
	for(int i=1;i<N;i++)	w[i] = -W[i] / sqrt(1-(pow(x[i],2))) ;
	double sum = 0.0;
	for(int i=0;i<N;i++)	sum = sum + (pow(i,2)*U[i]);
	w[0] = sum/N + 0.5*N*U[N];
	sum = 0.0;
	for(int i=0;i<N;i++)	sum = sum + (pow(i,2)*pow(-1,i+1)*U[i]);
	w[N] = sum/N + 0.5*pow(-1,N+1)*N*U[N];
	
	Peek(w);
}

Result:
(Using $N=16$ with $v_j = v(x_j) = \exp(x_j)\times\sin(5x_j)$ where $x_j=\cos(j\pi/N)$.)

#### FFT C++ Code References
[Fast Fourier transform](https://rosettacode.org/wiki/Fast_Fourier_transform#C.2B.2B)

---
### Python (尚文)

Some auxiliary codes:

import numpy as np
import matlotlib.pyplot as plt

def alt(n):
    alt = []
    for i in range(n):
        alt.append((-1)**i)
    return np.array(alt)

def _alt(n):
    alt = []
    for i in range(n):
        alt.append((-1)**(i+1))
    return np.array(alt)

def consec(n):
    return list(range(n))

def cesnoc(n):
    cesnoc = []
    for i in range(n):
        cesnoc.append(-n+i)
    return cesnoc

def arrconsec(n):
    return np.array(list(range(n)))

#### Differentiation Matrix

def cheb(N):
    if N==0:
        return 0, 1
    
    x = np.cos(np.pi*np.linspace(0,1,N+1))
    c = np.array([2] + [1]*(N-1)  + [2]) * alt(N+1)
    X = np.outer(x, np.ones(N+1))
    dX = X-X.T
    D = np.outer(c, np.array([1]*(N+1))/c) / (dX + np.identity(N+1))
    D = D - np.diag(np.sum(D,axis=1))
    return D, x

#### FFT differentiation

def chebfft(v):
    N, = np.shape(v)
    N = N-1
    if N==0:
        return 0
    
    vl = v.tolist()
    
    x = np.cos(np.pi*np.linspace(0,1,N+1)) # len = N+1 = len(v)
    ii = arrconsec(N) # len = N
    V = np.array(vl + vl[1:N][::-1]) # len = N+1+N-1 = 2N
    U = np.real(np.fft.fft(V)) # len = 2N
    W = np.real(np.fft.ifft( 1j*np.array(consec(N) + [0] + cesnoc(N-1)) * U)) # len = 2N
    w = np.zeros(N+1) # len = N+1 = len(v)
    
    # note - the last entry of (N+1,)-array is [N] in python
    w[0] = sum( np.square(ii)* U[:N] )/N + 0.5*N*U[N] # sum of matrix
    w[1:N] = -W[1:N]/np.sqrt(np.ones(N-1)-np.square(x[1:N])) # use np 1-array
    w[N] = sum( _alt(N)*np.square(ii)*U[:N] )/N + 0.5*(-1)**(N+1)*N*U[N]
    
    return w

#### Chebyshev recursion

def list2cheb(f):
    #given discrete evaluations on cosine grids f_i, return the Chebyshev coefficients
    N, = np.shape(f)
    N = N-1
    f_l = f.tolist()
    f_ext = np.array(f_l + f_l[1:N][::-1])
    f_hat = np.fft.fft(f_ext)/(2*N)

    a = np.zeros(N+1, dtype=complex)
    a[1:N] = 2*np.real(f_hat[1:N])
    a[0] = np.real(f_hat[0])
    a[N] = np.real(f_hat[N])
    return a

def chebdiff_3(f):
    N, = np.shape(f)
    N = N-1
    a = list2cheb(f)
    b = np.zeros(N+1)
    b[N-1] = 2*N*a[N]
    for run in range(2,N):
        i = N+1-run
        b[i-1] = 2*i*a[i] + b[i+1]
        print("i=",i,"b[{}]=".format(i+1), b[i+1])
    b[0] = a[1]+.5*b[2]
    
    ddx_f_hat = np.zeros(2*N)
    ddx_f_hat[0] = b[0]
    ddx_f_hat[N] = b[N]
    ddx_f_hat[1:N] = 0.5*b[1:N]
    ddx_f_hat[N+1:] = 0.5*np.flip(b[1:N])
    
    ddx_f = np.real(2*N*np.fft.ifft(ddx_f_hat))[0:N+1]
    return ddx_f

#### Results

import matplotlib.pyplot as plt

xx = np.linspace(-1,1,201)
uu = np.exp(xx)*np.sin(5*xx)

plt.figure(figsize=(20,8))
plt.subplots_adjust(hspace=0.5)
plot_count = 1

for N in [10,20]:
    D, x = cheb(N)
    u = np.exp(x)*np.sin(5*x)
    error1 = np.dot(D,u) - np.exp(x)*(np.sin(5*x)+5*np.cos(5*x)) # differentiation matrix
    error2 = chebfft(u) - np.exp(x)*(np.sin(5*x)+5*np.cos(5*x)) # FFT
    error3 = chebdiff_3(u) - np.exp(x)*(np.sin(5*x)+5*np.cos(5*x)) # recusion
    
    plt.subplot(3,2,plot_count)
    plot_count += 1
    plt.plot(xx,uu,'-',x,u,'o')
    plt.axis([-1,1,-4,2])
    plt.title('u(x), N={}'.format(N))
    
    plt.subplot(3,2,plot_count)
    plot_count += 1
    plt.plot(x,error1,'o',label="Diffmat")
    plt.plot(x,error2,'*',label="FFT")
    plt.plot(x,error3,'+',label="Recursion")
    plt.legend()
    plt.grid(True)
    plt.title('error in first derivative, N={}'.format(N))
    plt.plot(x,error1)
    plt.plot(x,error2)
    plt.plot(x,error3)


---
### Julia Code(浩恩)

#### Differentiation Matrix

using LinearAlgebra
using Plots
using FFTW

# CHEB  compute D = differentiation matrix, x = Chebyshev grid

function cheb(N)
    if N == 0
        D = 0; x = 1;
        return
    end
    x = cos.(pi*(0:N)/N);                  # generate grid x_i = cos(i*pi/N) i=0~N
    c = [2; ones(N-1,1); 2].*(-1).^(0:N);  # c_i*(-1)^i 
    X = repeat(x,1,N+1);                    
    dX = X-X';                             # compute x_i - x_j
    D = (c*(1 ./c)')./(dX+(I(N+1)));       # off-diagonal entries
    D = D - diagm(vec(sum(D,dims=2)));     # diagonal entries
    return D, x
end

#### Chebyshev differentiation via FFT

# CHEBFFT  Chebyshev differentiation via FFT. Simple, not optimal.  
#          If f is complex, delete "real" commands.
#          f = f(x_i) is a column vector

function chebfft(f)
    N = length(f)-1;
    if N==0
        f_prime = 0;
        return
    end
    x = cos.(pi*(0:N)/N);
    ii = 0:N-1;
    f_extend = [f; f[N:-1:2]];                     # transform f -> theta
    f_hat = real(fft(f_extend));
    w = real(ifft(1im*[0:N-1;0;1-N:-1;].*f_hat));  # transform theta -> w = df/(d theta)
    
    # transform w -> f_prime
    f_prime = zeros(N+1,1);
    f_prime[2:N] = -w[2:N]./sqrt.(ones(N-1,1)-x[2:N].^2);    
    f_prime[1] = sum(ii.^2 .*f_hat[1:N])/N + .5*N*f_hat[N+1];     
    f_prime[N+1] = sum((-1).^(1:N).*ii.^2 .*f_hat[1:N])/N + .5*(-1)^(N+1)*N*f_hat[N+1];
    return f_prime
end

#### Chebyshev differentiation via FFT in another way

# CHEBFFT2  Chebyshev differentiation via FFT in another way  
#           If v is complex, delete "real" commands.
#           f = f(x_i) is a column vector

function chebfft2(f)
    N = length(f)-1;
    if N==0
        f_prime = 0;
        return
    end
    x = -sin.((-N:2:N)*pi/(2*N));

    # extend the function and take scaled FFT
    f_extend = [f; f[N:-1:2]];
    f_hat = fft(f_extend)/(2*N);

    # obtain the Chebyshev coefficients of f
    a = zeros(N+1,1);
    a[1] = real(f_hat[1]);
    a[2:N] = 2*real(f_hat[2:N]);
    a[N+1] = real(f_hat[N+1]);
    
    # compute the Chebyshev coefficients of f'
    b = zeros(N+1,1);
    b[N] = 2*N*a[N+1];
    for i = N:-1:3
        b[i-1] = 2*(i-1)*a[i] + b[i+1];
    end
    b[1] = .5*(2*a[2]+b[3]);
    
    # compute the Fourier coefficient of f'
    f_prime_hat = zeros(2*N,1);
    f_prime_hat[1] = b[1];
    f_prime_hat[2:N] = 0.5*b[2:N];
    f_prime_hat[N+1] = b[N+1];
    f_prime_hat[N+2:end] = 0.5*b[N:-1:2];
    
    # take rescaled ifft to get extend f'
    f_prime_extend = real(ifft(f_prime_hat)*2*N);
    f_prime = f_prime_extend[1:N+1]; # get the unextended f'
    
    return f_prime
end

#### Example
$u(x) = e^x\sin5x$ and $u'(x) = e^x(\sin5x+5\cos5x)$

(1) $N=10$

xx=-1:1/100:1;
uu=(exp.(xx).*sin.(5*xx));

(D,x)=cheb(10);
u = (exp.(x).*sin.(5*x));
exact_u_prime = exp.(x).*(sin.(5*x)+5*cos.(5*x));
error1 = D*u-exact_u_prime;

u_prime = chebfft(u);
error2 = u_prime - exact_u_prime;

u_prime2 = chebfft2(u);
error3 = u_prime2 - exact_u_prime;

p1 = plot(xx,uu,label="",title="u(x), N=10",ylim=(-4,2),yticks=-4:2:2);
p1 = plot!(x,u,seriestype = :scatter,label="");
p1 = plot!(size=(800,210));

p2 = plot(x,error1,label="method 1",title="error in u'(x)",marker=:circle,markersize=6);
p2 = plot!(x,error2,label="method 2",seriestype = :scatter,marker=:utriangle,color=:red,markersize=5);
p2 = plot!(x,error3,label="method 3",seriestype = :scatter,marker=:hexagon,color=:yellow);
p2 = plot!(size=(800,210));

plot(p1,p2,layout = 2,legend=(0.8,0.3))

Result :

(2) $N=20$

xx=-1:1/100:1;
uu=(exp.(xx).*sin.(5*xx));

(D,x)=cheb(20);
u = (exp.(x).*sin.(5*x));
exact_u_prime = exp.(x).*(sin.(5*x)+5*cos.(5*x));
error1 = D*u-exact_u_prime;

u_prime = chebfft(u);
error2 = u_prime - exact_u_prime;

u_prime2 = chebfft2(u);
error3 = u_prime2 - exact_u_prime;

p3 = plot(xx,uu,label="",title="u(x), N=20",ylim=(-4,2),yticks=-4:2:2);
p3 = plot!(x,u,seriestype = :scatter,label="");
p3 = plot!(size=(900,240));

p4 = plot(x,error1,label="method 1",title="error in u'(x)",marker=:circle,markersize=6);
p4 = plot!(x,error2,label="method 2",seriestype = :scatter,marker=:utriangle,color=:red,markersize=5)
p4 = plot!(x,error3,label="method 3",seriestype = :scatter,marker=:hexagon,color=:yellow);
p4 = plot!(size=(900,240));
plot(p3,p4,layout = 2,legend=(0.4,0.88))

Result :